[UPDATE: Problem also happens in unitary time evolution, see below]
Hello everybody,
I've come around some weird error, with wich I've been struggling for quite some time now. I am trying to implement a stochastic wave function approach (arXiv:quant-ph/9806026, see sec II.A) for the time-evolution in open quantum systems. This method combines the non-unitary time evolution, with a non-Hermitian Hamiltonian @@H = H{closedsystem} -i \gamma C^{\dag} C @@, of wave functions with the stochastic application of jump operators.
For this I use IQTensors for the wave functions and gates for the bond operators of a second order even-odd time evolution. The model is a Heisenberg model for using the SpinHalf site-set . Most of the times this gives the right results (compared to ED), but after longer times I get an Index Mismatch error:
----------------------------------------
IQIndexSet 1 =
IQIndex(ul,4,Link,724) <Out>
(long,1,Link,853) QN(-2)
(long,2,Link,565) QN(0)
(long,1,Link,745) QN(2)
IQIndex(S=1/2 2,2,Site,416) <Out>
(Up 2,1,Site,957) QN(1)
(decltype(nullptr) 2,1,Site,622) QN(-1)
IQIndex(ul,2,Link,411) <Out>
(long,1,Link,282) QN(1)
(long,1,Link,761) QN(-1)
----------------------------------------
IQIndexSet 2 =
IQIndex(uc,4,Site,416) <Out>
(c0,1,Site,216) QN(2)
(c1,2,Site,631) QN(0)
(c2,1,Site,174) QN(-2)
IQIndex(S=1/2 2,2,Site,416) <In>
(Up 2,1,Site,957) QN(1)
(decltype(nullptr) 2,1,Site,622) QN(-1)
IQIndex(ul,2,Link,411) <In>
(long,1,Link,282) QN(1)
(long,1,Link,761) QN(-1)
----------------------------------------
Mismatched IQIndex IQIndex(S=1/2 2,2,Site,416) <Out>
(Up 2,1,Site,957) QN(1)
(decltype(nullptr) 2,1,Site,622) QN(-1)
A backtrace of the error shows, that the occurrence is in the position() function, during the SVD in the normalization process (see below). I also verified this by using couts right before and after the function call. Interstingly, the error is not deterministically reproducible. If I store the MPS just before entering the position() function where it crashes with writeToMPS(), reload it and perform position() everything works as expected. Also the point where it crashes is different for different runs of the same parameter set (incl. same seed). Therefore, I have thought the problem may originate in some dynamical memory errors. To circumvent this I tried to define the MPS explicitly on the heap
MPSt<IQTensor> *psi = new IQMPS(psi0);
where psi0 is a InitState, but that didn't solve the problem.
Backtrace of error:
[cmnode006:296505] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0xf030) [0x2b9079611030]
[cmnode006:296505] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x35) [0x2b9079850475]
[cmnode006:296505] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x180) [0x2b90798536f0]
[cmnode006:296505] [ 3] condor_exec.exe() [0x9f1dcb]
[cmnode006:296505] [ 4] condor_exec.exe(itensor::detail::checkArrows(itensor::IndexSetT<itensor::IQIndex> const&, itensor::IndexSetT<itensor::IQIndex> const&, bool)+0x3cc) [0xa7400c]
[cmnode006:296505] [ 5] condor_exec.exe(itensor::ITensorT<itensor::IQIndex>::operator*=(itensor::ITensorT<itensor::IQIndex> const&)+0x85d) [0xa8175d]
[cmnode006:296505] [ 6] condor_exec.exe(itensor::Spectrum itensor::svd<itensor::ITensorT<itensor::IQIndex> >(itensor::ITensorT<itensor::IQIndex>, itensor::ITensorT<itensor::IQIndex>&, itensor::ITensorT<itensor::IQIndex>&, itensor::ITensorT<itensor::IQIndex>&, itensor::Args)+0x466) [0xa079d6]
[cmnode006:296505] [ 7] condor_exec.exe(itensor::Spectrum itensor::orthMPS<itensor::ITensorT<itensor::IQIndex> >(itensor::ITensorT<itensor::IQIndex>&, itensor::ITensorT<itensor::IQIndex>&, itensor::Direction, itensor::Args const&)+0x24a) [0xb8ecea]
[cmnode006:296505] [ 8] condor_exec.exe(itensor::MPSt<itensor::ITensorT<itensor::IQIndex> >::position(int, itensor::Args const&)+0xc2) [0xb97eb2]
[cmnode006:296505] [ 9] condor_exec.exe(local_magnetization(itensor::MPSt<itensor::ITensorT<itensor::IQIndex> >, itensor::BasicSiteSet<itensor::SpinHalfSite>, int)+0x4f) [0xa051ef]
Also I once got the error message
Index::write: Index is default initialized
This might be an effect of the commonIndex() function used in orthMPS(). Generally I see the error appearing in every kind of parameter space, but it happens earlier (in terms of runtime) for smaller systems with lower cutoff, that means it might depend loosely on the number of gates applied. I tried this also for a Hamiltonian system, i.e. without Jump operators, but the program still crashes.
Maybe someone has any idea about that? I would appreciate any kind of help.
Best,
Stefan
P.S.: for completeness, I attach part of the code:
MPSt<IQTensor> *psi = new IQMPS(psi0);
*psi *= 1./norm(*psi);
compute_quantum_trajectory(sites,*psi,tensor_gates,
ofile_name,N, jump_operators,
uni,args, nstart, nfinal,Measurements);
void compute_quantum_trajectory(SpinHalf &sites,
MPSt<IQTensor> &psi,
vector< pair< int , IQTensor > > &tensor_gates,
string ofile_name,
int N,
vector<JumpOperator> &jump_operators,
my_uni_RNG uni,
Args args,
int nstart,
int nfinal,
vector<TimeSeries> &Measurements
)
{
int meas_counter =0;
double eta = uni(); // draws a random number from boost RNG
for(int step = nstart; step < nfinal; ++step)
{
measure(psi, sites, Measurements,step,ofile_name,args);
do_non_unitary_time_evolution(tensor_gates,psi,args);
if (norm(psi) * norm(psi) >= eta) continue;
else
{
JumpOperator Li
= choose_quantum_jump_operator(psi,sites,jump_operators,uni,args);
psi.Anc(Li.site) *= sites.op(Li.operator_string,Li.site);
psi.Anc(Li.site).noprime();
psi *= 1./ norm(psi);
}
eta = uni();
}
}
void
do_non_unitary_time_evolution(vector< pair< int , IQTensor > > &tensor_gates,
MPSt<IQTensor> &psi,
Args args)
{
for(auto& G : tensor_gates)
{
auto b = G.first;
psi.position(b,args);
IQTensor AA = psi.A(b)*psi.A(b+1);
AA *= G.second;
AA.noprime();
IQTensor D;
svd(AA,psi.Anc(b),D,psi.Anc(b+1), args);
psi.Anc(b+1) *= D;
}
}
[UPDATE]
As mentioned by G. Misguich, the problem also appears for the Hamiltonian time evolution. I tested it with a code very similar to the one provided itensor/tutorial/05_gates/gates.cc.