search

TRIQS / cthyb

A fast and generic hybridization-expansion solver
https://triqs.github.io/cthyb
Other
21 stars 23 forks source link

cthyb crashes by saying Delta(infty) is not real #93

Closed asubedi closed 4 years ago

asubedi commented 6 years ago

I am trying to reproduce the results of https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.214119, which were obtained last year using git master of TRIQS v1.4 series.

In LiVO2, the V 3d Hamiltonian has large off-diagonal terms due to on-site e_g-e_g overlaps. However, the output of dmftproj shows that there are no imaginary terms in the Hamiltonian.

I modified the script that I used for TRIQS v1.4 to reflect the changes in the API for TRIQS v2.0. The script is attached. The h5 file after Wein2kConverter is run is also attached, in case you want to run the thing by yourselves. [The extension log is used because github does not seem to allow attaching python files.]

LiVO2-el1-py.log LiVO2-el1-h5.log

The script crashes with the following error:

Traceback (most recent call last):
  File "LiVO2-el1.py", line 124, in <module>
    S.solve(h_int=h_int, **p)
  File "/home/subedi/softwares/TRIQS-2.0-git/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve
    solve_status = SolverCore.solve(self, **params_kw)
Traceback (most recent call last):
  File "LiVO2-el1.py", line 124, in <module>
    S.solve(h_int=h_int, **p)
  File "/home/subedi/softwares/TRIQS-2.0-git/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve
    solve_status = SolverCore.solve(self, **params_kw)
RuntimeError: .. Error occurred at Wed Jun 27 00:44:14 2018

.. Error .. calling C++ overload 
.. solve() -> void 
.. in implementation of method SolverCore.solve
.. C++ error was : 
.. Triqs runtime error at /home/subedi/softwares/TRIQS-2.0-git/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 126

Delta(infty) is not real. Maximum imaginary part is 1.19357e-06
.. Error occurred on node 1

I am not sure what has changed between versions 1.4 and 2.0 in CTHYB that the same input crashes the program.

Best, Alaska

parcollet commented 6 years ago

The tail fit, the threshold is probably too strict. @Wentzell : assigned

parcollet commented 6 years ago

@Wentzell : can you also improve the error message ? Max Im part is given, but what was the threshold.

parcollet commented 6 years ago

Meanwhile, you can simply increase the number of Matsubara freqs, it should solve this normally ...

asubedi commented 6 years ago

Hi Olivier, I tried n_iw = 6*1025. However, I still get the error. In fact, the Max Im part is slowly increasing with n_iw.

Best, Alaska

HugoStrand commented 6 years ago

@parcollet this is the issue I mentioned last week. I have relaxed the imaginary threshold on cthyb/unstable and made it into a warning.

@asubedi, could you please try the unstable branch and let us know if that gives reasonable results?

Best regards, Hugo

asubedi commented 6 years ago

@HugoStrand I still get the same error:

Starting on 32 Nodes at : 2018-06-28 19:08:16.921489
WARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in ta
ilfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error 
in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big e
rror in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: 
Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitWARNING: Big error in tailfitTraceback (most recent call last):
  File "LiVO2-el1.py", line 124, in <module>
Traceback (most recent call last):
  File "LiVO2-el1.py", line 124, in <module>
    S.solve(h_int=h_int, **p)
  File "/home/subedi/softwares/TRIQS-2.0-git-unstable/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve
    S.solve(h_int=h_int, **p)
  File "/home/subedi/softwares/TRIQS-2.0-git-unstable/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve
    solve_status = SolverCore.solve(self, **params_kw)
RuntimeError: .. Error occurred at Thu Jun 28 19:08:56 2018

.. Error .. calling C++ overload 
.. solve() -> void 
.. in implementation of method SolverCore.solve
.. C++ error was : 
.. Triqs runtime error at /home/subedi/softwares/TRIQS-2.0-git-unstable/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 127

Delta(infty) is not real. Maximum imaginary part is 5.11958e-06
.. Error occurred on node 15

    solve_status = SolverCore.solve(self, **params_kw)
RuntimeError: .. Error occurred at Thu Jun 28 19:08:56 2018
asubedi commented 6 years ago

I tried couple more tests.

First I tried to reproduce the calculations on LuNiO3 (https://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.075128). The is a two band problem (just the e_g electrons). CTHyb again crashes with the same error.

I thought maybe this has to do with the off-diagonal term in the Hamiltonian and/or the ensuing sign problem. [However, in both LiVO2 and LuNiO3, the average sign is greater than 0.5.]

So I tried the SrVO3 example given in the dft_tools documentation. But this time I artifically created a 3x3 block Green's function by doing SK.analyse_block_structure(threshold=1e-17). CTHyb crashes with the same error after two iterations. The average sign is 0.999967.

It looks like there is some issue with how the tail is handled when there are off-diagonal terms in the block Green's function.

Best, Alaska

asubedi commented 6 years ago

@mzingl and @aichhorn Looks like you have done calculations on LaVO3, which should have off-diagonal terms in the Hamiltonian. Can you please repeat them with CTHyb master?

parcollet commented 6 years ago

The issue here is that these tests should be in some C.I. (Continuous integration) for automatic testing.

We did not detect these issues before, simply because the tests are not run automatically (probably, we need so-called "hypothesis testing" as well here). This is something that Nils is working on, we expect to have this set up early Fall.

@Wentzell : can you add this on your list ? @mzingl, @aichhorn : if you run more tests, please coordinate with Nils to store them and have them in the CI.

HugoStrand commented 6 years ago

Dear all,

Here is a short update on what is happening with cthyb/master. There are, as you might know, issues with the rewritten transformation between imaginary time and Matsubara frequencies in triqs. This is being worked on.

Unfortunately this means that cthyb/master is suffering from these issues. When triqs have a production ready transform and we successfully can run the multiband cases discussed in this thread, I plan to bump up the version number of cthyb and make an official release.

Best regards, Hugo

sensudeshna commented 5 years ago

Dear TRIQS developers,

I am wondering if this issue is resolved. I bumped into this error when I am trying to perform DMFT self-consistency equations for a general density of states other than the Bethe lattice. For example constructing the S.G0_iw from the Dyson's equation involving the lattice Green's function and the solver provided self-energy. I write the code snippet below, in case I have made some mistake in feeding the S.G0_iw back to the solver:

for i in range(0,2*S.n_iw): g0.data[i,0,0] = 1.0/(sigma[i]+(1.0/G2.G_latt[i]))

S.G0_iw['up'].data[:,0,0] = g0.data[:,0,0]
S.G0_iw['down'].data[:,0,0] = S.G0_iw['up'].data[:,0,0]

where, G2.G_latt returns the Hilbert transform respective to the provided density of states for which I am currently using the semicircular density of states.

I have also carefully checked that the g0.data[:,0,0] matches perfectly with the g0 obtained on using the standard Bethe lattice relation: g0 << inverse(iOmega_n - e_f - g* t**2/4).

However, the subsequent iteration fail returning the following error:

WARNING: Big error in tailfitTraceback (most recent call last): File "python_class_obj.py", line 217, in main() File "python_class_obj.py", line 155, in main fit_max_n = 200) File "/home/sudeshna/triqs/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve solve_status = SolverCore.solve(self, **params_kw) RuntimeError: .. Error occurred at Thu Mar 21 17:10:35 2019

.. Error .. calling C++ overload .. solve() -> void .. in implementation of method SolverCore.solve .. C++ error was : .. Triqs runtime error at /home/sudeshna/codes/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 126

Delta(infty) is not real. Maximum imaginary part is 1.63529e-05 .. Error occurred on node 0

Your comments in this regard would be highly appreciated.

Best Regards, Sudeshna

HugoStrand commented 5 years ago

Dear @sensudeshna,

Thank you for reporting your problem! I am not sure that the origin of your problem is the same as the original posting of this issue.

If you can post your problem as a separate issue it would be great.

Regarding the error:

WARNING: Big error in tailfitTraceback (most recent call last):
File "python_class_obj.py", line 217, in 
main()
File "python_class_obj.py", line 155, in main
fit_max_n = 200)
File "/home/sudeshna/triqs/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve
solve_status = SolverCore.solve(self, **params_kw)
RuntimeError: .. Error occurred at Thu Mar 21 17:10:35 2019

.. Error .. calling C++ overload
.. solve() -> void
.. in implementation of method SolverCore.solve
.. C++ error was :
.. Triqs runtime error at /home/sudeshna/codes/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 126

Delta(infty) is not real. Maximum imaginary part is 1.63529e-05
.. Error occurred on node 0

In general this error is generated if the input G0_iw to the solver is not properly conditioned, wrt., noise level and number of frequencies. What the solver is doing under the hood is to split G0_iw in a static quadratic term Delta(infty) and a dynamic hybridization term Delta_iw that is then Fourier transformed to imaginary time Delta_tau.

The quadratic part Delta(infty) is not not purely real as expected by the real valued version of the solver. If this local crystal field of your model is supposed to be real then this is an numerical error, if it supposed to be complex you have to compile triqs/cthyb with complex support.

Best, Hugo

sensudeshna commented 5 years ago

Dear TRIQS/cthyb team,

Thanks for you response. You will notice that I have opened a new issue. It would be very helpful if you could elaborate on the above reply in that issue. Just to clarify, I am trying to do something very straightforward here. The non-interacting density of states is a semi-circular density of states. My final intention is to perform DMFT calculations for a model system that may not have a BL DOS, (but a 3D cubic DOS for e.g.) for which the exact relation g0 << inverse(iOmega_n -e_f - t*2g) is of course not correct. Looking forward to your reply.

Best Regards,

Sudeshna

On Thu, 21 Mar 2019 at 18:06, Hugo U.R. Strand notifications@github.com wrote:

Dear @sensudeshna https://github.com/sensudeshna,

Thank you for reporting your problem! I am not sure that the origin of your problem is the same as the original posting of this issue.

If you can post your problem as a separate issue it would be great.

Regarding the error:

WARNING: Big error in tailfitTraceback (most recent call last): File "python_class_obj.py", line 217, in main() File "python_class_obj.py", line 155, in main fit_max_n = 200) File "/home/sudeshna/triqs/install/lib/python2.7/site-packages/triqs_cthyb/solver.py", line 130, in solve solve_status = SolverCore.solve(self, **params_kw) RuntimeError: .. Error occurred at Thu Mar 21 17:10:35 2019

.. Error .. calling C++ overload .. solve() -> void .. in implementation of method SolverCore.solve .. C++ error was : .. Triqs runtime error at /home/sudeshna/codes/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 126

Delta(infty) is not real. Maximum imaginary part is 1.63529e-05 .. Error occurred on node 0

In general this error is generated if the input G0_iw to the solver is not properly conditioned, wrt., noise level and number of frequencies. What the solver is doing under the hood is to split G0_iw in a static quadratic term Delta(infty) and a dynamic hybridization term Delta_iw that is then Fourier transformed to imaginary time Delta_tau.

The quadratic part Delta(infty) is not not purely real as expected by the real valued version of the solver. If this local crystal field of your model is supposed to be real then this is an numerical error, if it supposed to be complex you have to compile triqs/cthyb with complex support.

Best, Hugo

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/TRIQS/cthyb/issues/93#issuecomment-475342499, or mute the thread https://github.com/notifications/unsubscribe-auth/AudusNoX4fjSkOhrLitc1ysbkQzuk-Udks5vY8o9gaJpZM4U40MC .

HugoStrand commented 5 years ago

The tail fit issues has been adressed in Triqs v.2.1.0 and the original tests should hopefully work now.

I am closing this issue due to inactivity, please reopen if needed.

asubedi commented 5 years ago

This is indeed fixed in version 2.1. Thanks!

the-hampel commented 4 years ago

I think I have to reopen this issue. I think the solver is not behaving as expected. When using DFT projectors to generate an input G0 it can happen, due to numerical convergence issues in Vasp, that the local Hamiltonian is indeed slightly complex. For example here for a t2g model:

Local Hamiltonian:
  Shell 1
    Site 1 (real | complex part)
    -0.9072065    -0.0003677    -0.0000003 |    -0.0000000    -0.0000054     0.0000014
    -0.0003677    -0.8659611     0.0000012 |     0.0000054     0.0000000    -0.0000065
    -0.0000003     0.0000012    -0.8986264 |    -0.0000014     0.0000065    -0.0000000

I encounter complex off-diag terms in Hloc that have a size of 10-5-10-6. We should of course also try to fix this in the Vasp interface, as these terms should not occur. However, from the solver docu I see that the solver_parameter imag_threshold should eliminate such elements if I set the threshold accordingly, e.g. 1e-4: "Threshold below which imaginary components of Delta and h_loc are set to zero" is written in the Doc. But, I obtain an error:

.. Error .. calling C++ overload 
.. solve() -> void 
.. in implementation of method SolverCore.solve
.. C++ error was : 
.. Triqs runtime error at /home/ahampel/git/triqs-dev/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 141

Delta(infty) is not real. Maximum imaginary part is 6.51749e-06
.. Error occurred on node 48

If I understand this correctly, in https://github.com/TRIQS/cthyb/blob/unstable/c++/triqs_cthyb/solver_core.cpp in line 141. The solver checks first for such off-diagonal elements with a fixed threshold of 1E-6 and later sets according elements to zero?

At least as a user this behavior seems not correct. Is there a reason for this, or should we set this threshold not fixed to 1E-6 but rather to solver_params.imag_threshold? Because later the elements are set to zero anyway. Maybe this should be rather a warning?

The input h5 file including the solver input parameters and the G0_iw can be found here: https://www.dropbox.com/s/8eutb8h78xphhwc/vasp_err.h5?dl=0 . The solver object is stored under DMFT_input/solver/it_1/S_0 (comment: I checked both in solver.py and solver_core.cpp just before the solver runs params.imag_threshold and it is properly set to the value I defined...)

Wentzell commented 4 years ago

Hey @the-hampel,

I agree! The imag_threshold should indeed be used also in the initial check of the imaginary part of Delta(infty). Maybe it should in fact be the only threshhold? My suggestion would be to

  1. Raise the default value for the imag_threshold from 1e-15 to e.g. 1e-14 or 1e-13 to allow some room for minor numerical inaccuracies.
  2. Set h_loc / Delta components smaller than imag_threshold to zero. For the non-complex version of cthyb throw an exception if either of the two is larger than imag_threshold (and not 1e-6 as it is now .. ).
  3. In the exception message, mention the imag_threshold parameter to the user so he will know what to adjust.

What do you think? Would this approach be too susceptible to error messages due to the delta infty extraction through fitting?

the-hampel commented 4 years ago

I created a PR that should provide all the mentioned improvements. I lowered the default threshold to 1e-13, added a more useful error message, and adjusted the documentation. I tested with the above given h5 archive the various options, and now the solver behaves as documented. I hope all looks good.

Wentzell commented 4 years ago

Fixed by 5a632cd