Bug in complex cthyb: imaginary offdiagonal elements in hloc ignored

RobertTriebl commented 6 years ago

Dear all,

I think I just discovered a bug in the complex CTHYB solver together with @gkraberger.

In the file solver_core.cpp line 110 in the unstable branch the maximal imaginary part is calculated from the real part

 for (int b : range(gf_struct.size())) max_imag = std::max(max_imag, max_element(abs(real(_G0_iw[b].singularity()(2)))));

Consequences: When one has offdiagonal elements in H0_loc that are all purely imaginary, they are ignored.

Robert

mferrero commented 6 years ago

That indeed looks like a bug. @Wentzell , @HugoStrand , @parcollet do you agree?

parcollet commented 6 years ago

What is the consequence of this ?

HugoStrand commented 6 years ago

Dear Robert,

Thank you for reporting. I am trying to write a small test that reproduces the problem (before fixing the bug). If you have some simple example code that reproduces it would be helpful as reference.

Cheers, Hugo

gkraberger commented 6 years ago

Hi,

Olivier, the consequence is that the quadratic part of the local Hamiltonian (i.e., what comes from G0_iw) is ignored in the CTHYB calculation if it only consists of imaginary elements (real part = 0 for all elements).

Hugo, although I'm not Robert, I can give you an example:

import numpy as np
from pytriqs.gf import *
from pytriqs.operators import *
from pytriqs.applications.impurity_solvers.cthyb import Solver
Hloc = np.array([[0,0.1j],[-0.1j,0.0]])
S = Solver(beta=40,gf_struct={'ud':[0,1]})
S.G0_iw << inverse(Omega-Hloc)
Hint = Operator()
S.solve(h_int = Hint,n_cycles=1)

This writes out (using the complex version of the solver, of course), among other things,

The local Hamiltonian of the problem:
0

As soon as at least one real part of any element of Hloc is sufficiently large (10^-14 is enough), it works. I.e, replacing Hloc by

Hloc = np.array([[1.e-14,0.1j],[-0.1j,0.0]])

gives

The local Hamiltonian of the problem:
(0,0.1)*c_dag('ud',0)*c('ud',1) + (0,-0.1)*c_dag('ud',1)*c('ud',0)

Cheers, Gernot

parcollet commented 6 years ago

Ok, my question is then the impact on previous computations.

it is.a quite rare case
the error will be large and easily spotted. Correct ?

RobertTriebl commented 6 years ago

I think the impact is very low, one has to be careful in specific cases, though.

Rare case: Purely imaginary Hloc are indeed rare. In a real material the bug is, I think, never important, since one will have some (real) diagonal elements in Hloc as on-site energies. If one is doing models of degenerate orbitals, however, the bug might appear since then one usually shifts the energy so that the on-site energies are all zero. If one adds then a purely imaginary off-diagonal element (can happen with spin-orbit coupling or magnetic field), Hloc has the form as in Gernot's example above and the bug appears.
Easily spotted: One might overlook the error if one is sloppy. The self-energy has the correct structure since there is still a hybridisation at the positions in the matrix where Hloc is wrong. So just from the self-energy, the error is not obvious. There are, however, multiple ways to spot the error: First, simply look at Hloc in the output file, what people usually do at some point. Second, when plotting a correlated band structure after a calculation with small interactions, one should realize that the bands are wrong. To give a physical example: In some models, the bug leads to exactly the same results with and without spin-orbit coupling (if the matrix elements are purely imaginary). Therefore, something is obviously wrong.

HugoStrand commented 6 years ago

Dear all,

I have pushed a test and a fix on unstable Robert could you please test that the fixes works for you?

Best, Hugo

HugoStrand commented 6 years ago

I have done a backport of the bugfix for the master branch, see above.

Robert and Gernot can we close this issue?

gkraberger commented 6 years ago

In master, everything works fine now. But if the solver is compiled in the real version and one passes a G0_iw with a complex Hloc, the complex entries are just ignored and no error is thrown (I ran into this problem when testing because I forgot to turn on the complex compilation). This can lead to false results if one does not check whether the Hloc used in the solver is correct (the same is true for unstable).

In unstable, four tests fail

The following tests FAILED:
         15 - atomic_observables (Failed)
         18 - measure_static (Failed)
         23 - complex_bug81 (Failed)
         24 - complex_Gtau_ED (Failed)

Interestingly, in the LastTest.log file, in the part for the complex_bug81, the reported Hloc at the beginning is still the one with the missing off-diagonal elements, but the value of S.h_loc is totally different. To me, this looks like a data type issue.

23/24 Testing: complex_bug81
23/24 Test: complex_bug81
Command: "/usr/bin/python" "/temp/AG-CMS/TRIQS/complex/new_debian/unstable/src/cthyb/test/python/complex_bug81.py"
Directory: /temp/AG-CMS/TRIQS/complex/new_debian/unstable/build/cthyb/test/python
"complex_bug81" start time: Jan 26 14:48 CET
Output:
----------------------------------------------------------
Starting on 1 Nodes at : 2018-01-26 14:48:56.351035

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌─┐┌┬┐┬ ┬┬ ┬┌┐ 
 ║ ╠╦╝║║═╬╗╚═╗  │   │ ├─┤└┬┘├┴┐
 ╩ ╩╚═╩╚═╝╚╚═╝  └─┘ ┴ ┴ ┴ ┴ └─┘

The local Hamiltonian of the problem:
1*c_dag('ud',0)*c('ud',0) + 3*c_dag('ud',1)*c('ud',1)
Using autopartition algorithm to partition the local Hilbert space
Found 4 subspaces.

Warming up ...

Accumulating ...

[Rank 0] Timings for all measures:
Measure            | seconds   
Average sign       | 1.4332e-05
G_tau measure      | 0.000162534
Total measure time | 0.000176866
[Rank 0] Acceptance rate for all moves:
Move set Insert two operators: 0
  Move  Insert Delta_ud: 0
Move set Remove two operators: 0
  Move  Remove Delta_ud: 0
Move  Shift one operator: 0
[Rank 0] Warmup lasted: 0.0769272 seconds [00:00:00]
[Rank 0] Simulation lasted: 1.3214e-05 seconds [00:00:00]
[Rank 0] Number of measures: 1
Total number of measures: 1
Average sign: 1
h_loc_ref = (0+0.1j)*c_dag('ud',0)*c('ud',1) + 1*c_dag('ud',0)*c('ud',0) + 3*c_dag('ud',1)*c('ud',1) + (0-0.1j)*c_dag('ud',1)*c('ud',0)
S.h_loc = (10j)*c_dag('ud',0)*c('ud',0) + (30j)*c_dag('ud',1)*c('ud',1)
diff = (0+0.1j)*c_dag('ud',0)*c('ud',1) + (0-0.1j)*c_dag('ud',1)*c('ud',0)
Traceback (most recent call last):
  File "/temp/AG-CMS/TRIQS/complex/new_debian/unstable/src/cthyb/test/python/complex_bug81.py", line 62, in <module>
    check_h_loc(h_loc_ref, h_loc_mat, orb=orb)
  File "/temp/AG-CMS/TRIQS/complex/new_debian/unstable/src/cthyb/test/python/complex_bug81.py", line 53, in check_h_loc
    assert( diff.is_zero() )
AssertionError
<end of output>
Test time =   1.64 sec
----------------------------------------------------------
Test Failed.
"complex_bug81" end time: Jan 26 14:48 CET
"complex_bug81" time elapsed: 00:00:01
----------------------------------------------------------

Unfortunately, I don't have the time to dig into this right now, but I'll be happy to provide quick answers if necessary or to look more closely in two weeks.

gkraberger commented 6 years ago

I just couldn't stop looking at it :-).

The problem is, when the the real version of CTHYB is installed and one builds the complex version, for testing purposes, still the real version is used. I think I ran into that issue before, but I never think of it. After installing the complex version, all tests pass.

So, now just one possible improvement remains: warnings in the real CTHYB if complex quantities are passed.

RobertTriebl commented 6 years ago

I also tried now my script where the bug appeared first in addition to the tests. Everything works fine now. From my point, we can close the issue.

Thanks for fixing!

TRIQS / cthyb

Bug in complex cthyb: imaginary offdiagonal elements in hloc ignored #81