Miscellaneous questions/feature requests

[ec147]

Hi,

Just some quick questions about TRIQS/CT-HYB:

• I assume that by default a Monte-Carlo run starts from a configuration where all the orbitals are empty (?). Is there a way to start a Monte-Carlo run from a given configuration ? I am doing self-consistent run and thermalization is taking me a really long time, so this would help massively if I could start from a configuration of a previous iteration.
• When the density matrix is measured, is the time translation invariance taken into account ? Meaning, do you look only the configuration at time beta, or do you average over the whole segment [0,beta] ? The density matrix is very important for me to get a clean energy.
• Do you plan on releasing an improved estimator any time soon (ideally in the rotationally invariant picture and Legendre basis) ?

Again, all of this would help me massively. All of these issues are really impeding my work.

Thanks

[the-hampel]

Hi @ec147,

Do you plan on releasing an improved estimator any time soon (ideally in the rotationally invariant picture and Legendre basis) ?

Right now there is no concrete plan to implement improved estimators in triqs/cthyb. In case you are only interested in improved self-energies consider taking a look at: https://arxiv.org/abs/2310.01266 . Here, we showed how to leverage DLR + plus the measurement of high frequency moments: https://triqs.github.io/cthyb/latest/guide/high_freq_moments.html to achieve better noise suppression than standard improved estimators and the same noise level as symmetric improved estimators. I am currently working on a release version of the functionality within triqs for this, and I am happy to share a beta version of it. For the symmetric improved estimators we would require first a worm sampling algo in cthyb.

I assume that by default a Monte-Carlo run starts from a configuration where all the orbitals are empty (?). Is there a way to start a Monte-Carlo run from a given configuration ? I am doing self-consistent run and thermalization is taking me a really long time, so this would help massively if I could start from a configuration of a previous iteration.

As far as I know there is no way to initialize the startup configuration in MC right now and it is empty in the beginning

@Wentzell @parcollet please correct me if I am wrong. Maybe you can also answer the second question regarding the density matrix.

Best,
Alex

P.S. these questions are more qualified for a discussion. So far I do not see any obvious issue. So please consider next time to open this in a discussion instead.

[Wentzell]

Dear @ec147,

Thank you for reaching out and for your interest in the TRIQS software.

• I assume that by default a Monte-Carlo run starts from a configuration where all the orbitals are empty (?). Is there a way to start a Monte-Carlo run from a given configuration ? I am doing self-consistent run and thermalization is taking me a really long time, so this would help massively if I could start from a configuration of a previous iteration.

CTHyb samples the partition function.
The Monte-Carlo weight includes the atomic trace which is always calculated for the full Fock-Space of the atomic degrees of freedom, i.e. no particle-number constraint.
Regarding your point on using an existing configuration for the start of the Monte-Carlo run: This is currently not possible.

• When the density matrix is measured, is the time translation invariance taken into account ? Meaning, do you look only the configuration at time beta, or do you average over the whole segment [0,beta] ? The density matrix is very important for me to get a clean energy.

I am not sure I understand the question. The configurations always contain hybridization events (c, cdag operators) at random time points chosen from the full interval [0,beta]

P.S. these questions are more qualified for a discussion. So far I do not see any obvious issue. So please consider next time to open this in a discussion instead.

I agree on this point, I converted this to a discussion.

[HugoStrand]

Dear @ec147,

Let me try to reply on the (many-body) density matrix question.

When the density matrix is measured, is the time translation invariance taken into account ? Meaning, do you look only the configuration at time beta, or do you average over the whole segment [0,beta] ? The density matrix is very important for me to get a clean energy.

We currently do not use time-translational invariance when sampling the density matrix. We simply accumulate the entire product of atomic propagators and sandwiched operators (combined with the hybridization function determinant).

I think this amounts to the first case that you call "only looking at the configuration at $\beta$". If one thinks of the imaginary time interval $[0, \beta]$ as periodic (wrapping around at $\beta$), one can see this as only sampling the contribution to the many-body density matrix at the point $\beta$ (and $0$).

As you point out, one could "split" the matrix products of atomic propagators and operators at any other time $\tau$ on the interval $[0, \beta]$ and use this to get more samples per configuration. However, I am not convinced that this would improve the sampling since the samples from one configuration will be correlated.

@ec147 do you have any insights on this?

Cheers, Hugo

[ec147]

interaction_energy.pdf
nb_electrons.pdf

Thank you all for your replies.

• About the density matrix:

Yes, this is exactly what I mean, @HugoStrand. In order to sample the density matrix, you could wrap around any arbitrary time $\tau_i$ between 0 and $\beta$ instead of using simply $\beta$. Formally, this amounts to shift every creation and annihilation operators by $\beta - \tau_i$ , and sample the density matrix just as you are currently doing for this new configuration. This is justified because all of these time-shifted configurations have equal probability due to time-invariance.

You could do this for several $\tau_i$ and thus increase the number of density matrix samples for each configuration. Even better, you could even do a continuous average over [0, $\beta$ ] instead of only looking at a set of discrete $\tau_i$.

Besides TRIQS, I am using a segment solver, and all our quantities (number of electrons, double occupancies, configurations) are sampled by taking a continuous average over [0, $\beta$ ]. I agree that in the case of a matrix solver like TRIQS, doing this would require much more implementation effort than for a segment solver. However, you can expect this to massively improve statistics. Indeed, if you have a simple system with one orbital and an average occupation number of 1/2, the sum of all the segments will be $\beta$/2 on average. So, taking the continuous time average of the number of electrons (ie the sum of the lengths of the segments divided by $\beta$) will yield the correct probability 1/2 by looking only at one configuration. Looking only at time $\beta$ will yield 0 half of the time and 1 electron half of the time, so you would need to sample at least two different configurations to get the correct probability.

However, you really piqued my curiosity about the benefits of this, since I have actually never done any formal test to compare both these techniques. So I just did a self-consistent DFT+DMFT calculation with my segment solver and plotted the interaction energy (computed from the double occupancies directly sampled from my solver) and the number of electrons (also directly sampled from my solver) along the DFT+DMFT iterations using both sampling techniques. Both calculations are done with the exact same number of sweeps and the exact same parameters of course. You can see the drastic improvement obtained by time-averaging. I have indicated the standard deviations, and you can see this is an order of magnitude smaller when time-averaging. The two plots are attached (I am sorry for the horrible format ; for technical reasons, I could not directly give you the original file and I had to print it and then scan it). Now I realize that to you, I am just a random Internet weirdo, but you can easily verify my claim within a few minutes using your unpublished segment solver. I hope this will motivate you to at least look that up, and (hopefully) implement this in TRIQS.

• Also, while I'm at it, this might be a very basic question, but why exactly do we have to use the norm of the matrix as weight in order to measure the density matrix in TRIQS ?

• About the noise reduction of the self-energy

I have read your article @the-hampel, and I have a few questions. In fig 3, you plot the self-energies computed with improved estimators, symmetric improved estimators and constrained residual minimization using 10^8 sweeps, and compare them to a reference self-energy computed from directly sampled Legendre coefficients, using 10^11 sweeps (if I understand correctly). You show that you get a similar noise level as a symmetric improved estimator, which is really impressive. But what about Legendre and 10^8 sweeps ? This would mean that the best techniques right now to reduce the noise levels are (by decreased order of efficiency): cRM ~ sIE > iE > Dyson equation ? How does direct sampling of Legendre compare to those ? Do you have any insight on that ? If cRM is truly the best technique out there, I would be really interested in getting a beta version of it.

• Conclusion

Basically, what I am trying to do with TRIQS/CTHYB is getting a 1 meV precision on free energy for a 3d system (thus with 10 orbitals, accounting for spin), in the rotationally invariant picture. I am really close on getting that, but I am still encountering a lot of issues. I was trying to compare my results with those that K. Haule obtained using his own CTQMC solver from embedded DMFT, and he told me that he routinely manage to get this accuracy with merely 500 CPUs, while I struggle to get this accuracy with 10x more CPUs using TRIQS. So I was trying to understand what is so magical about his CTQMC, and the three points I mentioned in my initial comment were the main differences I could spot. For instance, he clearly states in his article (https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.155113) that he samples the time-average of his operators to get improved accuracy. Just know that there is at least one guy out there who would be VERY interested in getting a solution to these three points.

PS: I am sorry, from now on I will post in the discussion thread in order to ask my questions.

[the-hampel]

Hi,

We debated quite a bit about this comparison with the IE against the legendre sampled MC results. As you may noticed, we have rigorous comparison against analytic results as well, but to showcase IE etc we decided to use an example that is not exactly solvable. We just tried manually to create the "best" possible result by sampling using legendre coeff and carefully choosing a cut-off in coefficients to remove remaining MC noise. However, ideally we should have sampled G(tau) infinitely long to create the possible reference data to remove any noise and compare against this bare sampling, but to speed things up we used legendre as an noise filter (which may be inherently unstable if not done carefully). So it is not necessarily about legendre sampling. This still has in general the problem of an arbitrary cuf-off as long as MC noise is present. We also think one could implement the constrained fitting using legendre coefficients instead of DLR, but DLR is just superior when it comes to temperature scaling (similar to IR)

So I was trying to understand what is so magical about his CTQMC

Let me also comment on K. Haule's CTQMC implementation. In general it would very valuable to extend the comparison we did in the paper mentioned above to K. Haule's CTQMC code. Right now we compared in the paper implicitly how well triqs/cthyb compares to w2dynamics. We found very good agreement with respect to noise levels for given number of sweeps. Also run times where on the same magnitude. So if you or someone else has time to spare such rigorous comparison for a controlled example between QMC impurity solvers would be amazing to check progress and consistency in the community! Our paper is just a small starting point for that.

When it comes to speed I would like to point out that the CTQMC code of Haule often does not not use the full impurity fock space. You have to manually specify for (as far as I know) which particle sectors should be included etc. For larger impurity problems this greatly improves speed and also maybe statistics. Especially in f-systems this can be a real deal changer, when you have only one electron in the f-shell etc. So for a fair comparison of speed please check that CTQMC used also all particle sectors, or remove those in cthyb using loc_n_min and loc_n_max (https://triqs.github.io/cthyb/latest/_ref/triqs_cthyb.solver.Solver.solve_parameters.html) . But this feature should be used with great care, since you could throw away contributions to G(tau).

P.S.: "Now I realize that to you, I am just a random Internet weirdo, but you can easily verify my claim ..." Feel also free to use real names and or affiliation. ;-)

Best,
Alex

[ec147]

Hi,

So it is not necessarily about legendre sampling. This still has in general the problem of an arbitrary cuf-off as long as MC noise is present.

I have never claimed that Legendre sampling was perfect. However, on top of filtering some noise, it also provides an analytical Fourier transform and thus removes any aliasing due to the binning error. In practice, even by taking a large number of $\tau$, doing a Fourier transform on some binned data (either by fitting Legendre coefficients or using spline interpolation for instance) ends up getting me an accuracy which is an order of magnitude lower than directly sampled Legendre. So I was merely how well it compares to IE and cRM.

Right now we compared in the paper implicitly how well triqs/cthyb compares to w2dynamics. We found very good agreement with respect to noise levels for given number of sweeps. Also run times where on the same magnitude. So if you or someone else has time to spare such rigorous comparison for a controlled example between QMC impurity solvers would be amazing to check progress and consistency in the community!

Yes, but you compare self-energies. I was talking about the sampling of the density matrix, which is crucial to get a clean energy (which is what I'm interested in). Besides, I do not know whether time invariance is taken into account in w2dynamics to sample the observables.

When it comes to speed I would like to point out that the CTQMC code of Haule often does not not use the full impurity fock space.

In the system K. Haule studied and that I am trying to compare with, restricting the Fock space is out of question. All the states have non-negligible probabilities. Even though, I have tried removing the least likely states with loc_n_min, but not only this degraded a lot my Green's function (since, again, these states had non-negligible contributions), but the computation time did not change much. I want to point out that the fact that I need 10x more CPUs than him to get the same accuracy on the energy corresponds to the accuracy loss you get when not taking into account time invariance (at least in my system). I have just shown that in my previous comment. If you do not believe me, again, this claim can be verified within a few minutes by using a segment solver.

By proposing that you implement time invariance and the possibility of starting a Monte-Carlo run from a previous configuration, I was merely trying to give you a feedback of what one of your TRIQS users need. I thought this was the purpose of this github thread. Personally, I do not understand why you are so reluctant about implementing features which would only require a few days work and would drastically improve the accuracy and calculation time, but I cannot force you.

Now, since, again, I critically need these features for my work, I was thinking about implementing them by myself. I do not get the feeling that you would be willing to cooperate with me for that, but in the case you are, I would appreciate a few answers. So let me reiterate my questions:

• Why exactly do we have to use the norm of the matrix as weight when sampling the density matrix ?
• From my understanding of your code, it would thus seem to me that the sampling of the different configurations and thus the result of the Monte-Carlo run would be different when sampling the density matrix (since we use the norm of the matrix as weight instead of the trace). Isn't it an issue ?
• When deriving the equations of the CT-HYB solver, the weights of the MC configurations involve a trace. How do you justify the use of the norm as weight then ? Isn't the principle of a CT-HYB to sample the configurations in such a way that they reproduce the probability distribution given by the partition function of the Anderson Hamiltonian ?

Thanks

[aichhorn]

Thanks for the feedback ec147. You need to understand that there is no such thing as implementations on demand. If you do not get the feeling that somebody wants to cooperate with you on this, the reason might be that you are not even using your real name for this thread. This makes interactions really complicated.
Best,
Markus

[the-hampel]

I am not sure at which point you felt that we did not properly replied to your questions. Just let me try to explain that the people replied here work on very different parts of TRIQS, and and most of the contributions do come from people that do not work full time on TRIQS. In fact there are only 3-4 full time TRIQS developers. If you had the feeling I try to circumvent your questions, that is because I am not an expert on triqs/cthyb. I just tried to answer your questions to the best of my knowledge. When it comes to answering whether we can implement such feature you have to wait for the corresponding developers to respond. So please do not give up that quickly. Implementations in such big code take time and most importantly we first need to find someone who is willing to implement this and do the testing. You are more than welcome to give it a shot and make a PR.

And maybe there was some misunderstanding about my statement regarding speed. But you wrote:

I have just shown that in my previous comment. If you do not believe me, again, this claim can be verified within a few minutes by using a segment solver.

I do believe you. My statement was more a motivation. If you can proof that triqs/cthyb is 10x slower when sampling a certain quantity, this will very much motivate us to understand more why that is, i.e. in case sampling G(tau) is actually 10x faster in K. Haule's code we would be very interested in this. I was merely pointing out that we only compared so far against w2dynamics and found similar performance.

I guess if we want to improve accuracy in cthyb we first have to understand what K. Haule actually implements (this was not clear to me, and probably is not to the other developers here). Are you referring to Eq. 35 in https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.155113 ? It would have helped to point that out I think. @Wentzell @krivenko maybe you can take shot at answering briefly the question of @ec147 regarding using the norm of the matrix as weight. Maybe there is even a paper for this, and how this compares to K. Haule's approach described in his paper.

[ec147]

You need to understand that there is no such thing as implementations on demand.

This is not at all what I am asking, and apologies if this is the feeling you got from my comment. I know that you are all working on different projects. I was simply expecting either "Thanks for the feedback, this is has been noted and we will think about it for future releases", or "Thanks for the feedback, but this cannot be done at the moment because we are working on other projects (or other reasons)".

If you can proof that triqs/cthyb is 10x slower when sampling a certain quantity, this will very much motivate us to understand more why that is, i.e. in case sampling G(tau) is actually 10x faster in K. Haule's code we would be very interested in this. I was merely pointing out that we only compared so far against w2dynamics and found similar performance.

Indeed, I misunderstood you. Regarding the sampling of G($\tau$), it seems to me that you already do take into account time invariance. Not taking into account time invariance would mean that you remove only part of the hybridization lines when sampling. It seems to me that you remove all of them (?). Anyway, even though the comparison between a segment solver and a matrix solver is not clear, I did my best to have two (more or less) comparable runs between my segment solver and TRIQS, and I seem to have more or less the same accuracy regarding G($\tau$). Thus, to me, the only issue is the sampling of the density matrix. Now, I realize that the density matrix is not your main concern, and that your main focus is about the sampling of the Green's function. But I still think this could interest you, yes ?
Also, just to clarify myself, I do not claim that this will increase tenfold the accuracy in all cases. This is the case for my system, but the accuracy gain will obviously strongly depends on the number of orbitals, perturbation order, etc. Still, I am convinced that all in all, this will drastically improve the accuracy of the density matrix no matter what.

Are you referring to Eq. 35 in https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.155113 ?

Yes, precisely. This is the continuous time average I was talking about.
However, this whole time invariance thing relies on the fact that time shifting every creation/annihilation operator by the same amount does not change the probability of the configuration. But it seems to me that this is no longer true when sampling with the norm of the matrix instead of the trace. Hence my questions about it.

I am sure you have a good reason for doing all this, especially since you state in an old TRIQS documentation that this allows to take into account very important contributions for the density matrix, but I am just trying to understand.

Best