New features for internal energy calculations

[ec147]

Greetings.

I have implemented several critical features in TRIQS/cthyb which were lacking to compute an internal energy.

The last time I tried to tell you about these features in the discussion thread, in the hope that you would at least agree to collaborate with me and answer my questions, (most of) you very kindly ignored me. Although you clearly show no interest at all in optimizing your code for energy calculations, you have to understand that some users would appreciate to be able to do that. So far, you are the only one to provide an automatic partition algorithm, so you are supposedly the fastest CT-QMC out-there for rotationnally invariant Slater parametrization, and this would really be a shame to be forced to use the slower embedded-DMFT in order to be able to compute an energy. Thus, I would really like my work to benefit everyone, so I am doing a pull request. I am not asking that you welcome these with open arms (since you obviously won't), but I have not changed the default behavior at all. A simple backdoor with some undocumented variables to activate my features would be enough. My custom interface with TRIQS/cthyb would do the rest.

Here's a non-exhaustive summary of the features I have implemented (this is by no means definitive, I am willing to discuss some implementation details and answer your questions, unlike you):

• Use of time invariance when sampling the density matrix -> in my system, the standard deviation of the density matrix was divided by 10, for simply twice the computation time, as I have claimed in the previous discussion I had with you.

• Possibility of starting the Monte-Carlo from a given configuration, so that thermalization can be skipped in a self-consistent calculation.

• Use of linear interpolation instead of 0th order interpolation to compute $\Delta(\tau)$. Right now, this is disabled so that I can pass the tests, but I put a flag to activate it. Ideally a cubic spline interpolation on a non-uniform grid would be better, but this is much harder to implement. I have tried DLR representation, which would theoretically remove the need for any interpolation scheme, but it simply takes too much time. Also, I do not see why you enforce the same grid for $\Delta(\tau)$ and $G(\tau)$, so I have added a parameter n_tau_delta specifically for this (the default value is n_tau).

• Possibility of providing a user-defined proposal distribution for the insert and remove moves. By default, they are not optimized at all in TRIQS/cthyb. You simply choose a random time interval on [0, $\beta$], while it has been known for a long time that only moves with small time intervals are accepted in general. In TRIQS/cthyb, you have the philosophy of rejecting moves quickly rather than proposing them cleverly, "quickly" does not mean instantaneously, and the latter philosophy is always much more elegant than the former. Thus, in a self-consistent calculation, I use the histogram of the performance analysis feature to estimate the acceptance probability for each $\tau$, and then use it as my proposal distribution for the next iteration. On my system, this has massively increased the acceptance rate of these two moves (from as low as 1% at low temperatures up to ~30%). More importantly, the acceptance rate is now constant with temperature, which is extremely important for low temperature calculations.

• Instead of forcing the user to manually tune a length_cycle parameter to decorrelate the measurement, I do not understand why you don't simply store the measure until the next update, then multiply it by the number of times the run has been stucked in this configuration. So far, I have only implemented this for the density matrix, since this is the real bottleneck for me, and this is much harder to do for the measurements of $G(\tau)$ with the way your code is structured. Also, I print the average update time at the end of the calculation, which is much more meaningful than an auto-correlation time measured from the perturbation order and the sign. (The shift move does not change the perturbation order, yet does it mean that it doesn't decorrelate the configurations ? In the density-density approximation, the sign is constant, yet does it mean that every single configuration in the density-density approximation is correlated ?)

• In qmc_data.hpp, you did not reserve the space for the dets variable before using emplace_back. This caused some issues, and the variables singular_threshold, precision_warning, etc. were kept at their default values instead of being attributed the user-provided values.

• In impurity_trace.hpp, the formula for the lazy trace evaluation is too tight. The correct upper bound for the trace is $dim(B) e^{-\Delta \tau E_{min}}$ and not $\sqrt{dim(B)} e^{-\Delta \tau E_{min}}$ as it is currently implemented in TRIQS/cthyb. This is eq 13 in this paper. This is even written in a comment a few lines before the formula.

• The threshold for the deviation of the trace of the density matrix from 1 at the end of the calculation is overly aggressive. On my 10-orbital system, at low temperature, the deviation was a bit over the tolerance, simply due to the size of the configurations. Thus it was preventing from doing my calculation, and I deleted it to simply leave a warning for now. Again, I am willing to discuss that.

• No idea why you do not allow the user to use the trace as weight when sampling the density matrix. You claim that using the norm of the matrix is extremely important in that case, but I have seen no difference in the results of my calculations between these two methods. Besides, to my knowledge, no other CT-QMC does that, and I do not understand why these two methods would yield different results. Thus I have deleted the check enforcing use_norm_as_weight=True when sampling the density matrix. I would have liked to discuss that with you, if only you had agreed to answer my questions about it.

[HugoStrand]

Dear @ec147,

Thank you for sharing your work on TRIQS/cthyb. I have had a brief look at your changes and I think I have a vague understanding of some of the features you have implemented.

Could you please help me understand what your goal is with this pull request? and what expectations you have on the small community of old developers?

Best, Hugo

[ec147]

Sure.

First of, my original motivation for this work was to be able to compute an internal energy on big systems, using TRIQS/chyb.

Since DFT Tools has not been designed for that, I first started to build my own interface to TRIQS/cthyb (which I have not published yet, but I intend to). However, I still had a lot of issues due to high autocorrelation times, thermalization and a noisy density matrix. So I implemented the features I showed you in order to fix these.

Now, although some people will be able to maintain my interface with TRIQS/cthyb so that it stays up-to-date, no one will be able to maintain my version of TRIQS/cthyb, since I am the only one to understand it in my lab. I would really like my work not to go to waste, so here's the reason for my pull request.

Basically, I am showing you my work to see if you would be interested in collaborating with me. The way I see it, there are three possibilities:

• You are not interested at all in adding this in TRIQS/cthyb, and I publish my own version of TRIQS/cthyb anyway because I really want to provide to the users a way of computing an energy. In that case, this will not be maintained, so this will quickly die out.
• You are not interested, but still agree to add it to TRIQS/cthyb since this does not change the default behavior. In that case, no need to document these features nor add them to the Python API. They would be accessible only through my custom interface with some undocumented keywords.
• You are interested in my work, and want to document them and fully add them to TRIQS/cthyb. In that case, we would need to work through this pull request together.
  It's up to you.

Oh, and apologies for the tests I do not pass. I only use gcc, so I had no idea my code would not compile with different compilers. I'll fix that, in the case where you intend to merge with me.

Best

[ec147]

Usually, this is basic decency to at least acknowledge the person who is trying to improve your code.
This is a simple question: do you want to merge, yes or no ?

[HugoStrand]

Dear @ec147,

I am interested in the features you have implemented.

Do I understand you correctly that your expectations on the small community of old developers is that we should document and test your pull request?

Would you be willing to contribute to the testing and documentation?

Best, Hugo

[ec147]

Dear @HugoStrand,

I'm glad to see you're interested in my features. I am willing to contribute to the documentation and addition of new tests.

Unless you have specific guidelines, I'll start working on 3 new tests (time invariance, restart from a previous configuration
and user-provided proposal distribution).

There is also the issue of the linear interpolation for $\Delta(\tau)$. Either you want to keep it as an option, and this would also
require a new test; either you want to completely replace the 0th order interpolation by a linear interpolation,
and this would require modifying the reference of almost every current test.

I will also work on adding the new variables to the Python API, changing a few details, and on making sure my version passes all the tests.
In the meantime, I am happy to hear any suggestions.

Best

PS: If you are interested, we can also discuss by email.

[Wentzell]

Dear @ec147,

Thank you for your various efforts on improving the TRIQS/cthyb solver.

I agree with @HugoStrand that the features are very much of interest
to us and the broader TRIQS community.

It would be great if you could break your changes into several individual
Pull Requests, for each of the bullets you have raised above.
This would help a lot in making the review process more viable and
to disentagle the discussions on the various points.

Could you possibly look into that?

It would also be useful if you joined our next Zoom Developers Meeting
this Thursday July 25th at 10am ET. Please provide me with your Email
address so I can send you an invite.

Best,

Nils

[ec147]

Dear @Wentzell,

I'm glad you are interested in the features I have implemented.
No problem, I will break my changes into several smaller pull requests; but please give me a bit of time to do that.
Sure, I'll join your meeting. My email address is e.castiel@hotmail.fr.

Best

[ec147]

Closing this PR to split it into smaller ones.

[Wentzell]

Thank you @ec147 for starting to break your contributions down into managable chunks.
I noticed that you opened four new PR's. Are you still working on creating additional PR's
to cover the 9 bullet points you have listed above?

[ec147]

So far, my PR cover all my original modifications except bullets 2 and 6. Since the last three bullets are only one-line modifications, I have combined them with the first bullet, but I can do a separate PR for them if you prefer. Likewise, since bullet 6 is only one line and is related to a part of the code that is also modified by bullet 2, I think it's easier to merge them into a single PR, but it's up to you. I haven't done a PR for bullet 2 yet since adding it to the Python API is a bit tricky, so please let me a bit of additional time. And once we agree on the way these new features will be included and called in TRIQS/CTHYB, I will start working on the documentation and testing.

I have also squashed all my commits. Hope it's more to your liking.