[QE-developers] Forces are not exact derivatives of the energy for U_projection_type = 'ortho-atomic'

[Timrov Iurii]

Dear Max,


Thanks for performing these tests!


I implemented Hubbard forces and stress with ortho-atomic orbitals, and I have thoroughly tested the implementation. I compared forces and stress computed using the exact analytical formulas versus finite differences, and the agreement was up to 10^{-6} Ry/Bohr which validates the implementation. I used the symmetric difference quotient (with the displacement parameter of 5 x 10^{-3} Bohr) to computed derivatives numerically:

https://en.wikipedia.org/wiki/Numerical_differentiation

For example, for NiO I obtained the following results:

Analytical:   F = -0.048891 Ry/Bohr

Numerical:  F = -0.048891 Ry/Bohr

As you can see forces agree up to 10^{-6} Ry/Bohr. In a similar way I tested also stress.


> It seems that there is a bug in the implementation of the forces for U_projection_type = 'ortho-atomic’, or that some terms have been neglected when computing the derivatives.


The analytical formulas are exact, we did not neglect anything. In fact, it took us a lot of time and effort to derive the exact formulas. Also I do not think that there is a bug, because my tests give excellent agreement between analytical and numerical derivatives for NiO and a couple of other systems.


The path integral method looks to me a bit overkilling. Why not simply to compute the total energy for two atomic configurations and then simply divide the difference of energies by the difference in atomic positions (I mean only one atom must be slightly displaced at a time), i.e. [ E(x + delta) - E(x - delta) ] / (2*delta)? But in any case, since your method works well for standard DFT forces and DFT+U atomic forces, then we should try to understand why it gives weird results for the ortho-atomic case.


Which version of Quantum ESPRESSO did you use? Was it 6.6 without any modifications, or did you apply the patch which I sent to you in the QE developers mailing list shortly after the release of QE 6.6? Are you running on your workstation or on Piz Daint (should not matter but good to know)? Which compilers and libraries do you use?

Perhaps, you can try to increase the kinetic-energy cutoff and see if the agreement for the ortho-atomic case is better. Forces and stress converge slowly with respect to the cutoff (and maybe ortho-atomic case is even slower, but in my tests I did not observe this though). Also, maybe it would useful to try other compilers and libraries.

If you want I can send you my input files for NiO which I used for the benchmarking and we can see what do you obtain using the path integral approach. Would you be so kind to perform these tests? On my side I can compute forces with finite differences for your system - could you send to me please the input and output files for reference?


> In the long run, it might be helpful to include such tests routinely as part of the QE development process.


The QE developers team makes big effort to benchmark various functionalities of QE before the release and during the development process. In particular, we have a test suite with many examples which are used to benchmark various codes and various functionalities. For the DFT+Hubbard ortho-atomic forces, I tested everything using finite differences (as described above), I was convinced that all is good, and only then I added tests to the test suite. So we do care a lot about testing and benchmarking, but of course sometimes there maybe be some bugs which were overlooked or which comes for some other reason.


Greetings,

Iurii


--
Dr. Iurii TIMROV
Postdoctoral Researcher
STI - IMX - THEOS and NCCR - MARVEL
Swiss Federal Institute of Technology Lausanne (EPFL)
CH-1015 Lausanne, Switzerland
+41 21 69 34 881
http://people.epfl.ch/265334
________________________________
From: developers <developers-bounces at lists.quantum-espresso.org> on behalf of Max Amsler <amsler.max at gmail.com>
Sent: Saturday, September 12, 2020 12:44:33 PM
To: General discussion list for Quantum ESPRESSO developers
Subject: [QE-developers] Forces are not exact derivatives of the energy for U_projection_type = 'ortho-atomic'

Dear all,
After extensively testing the recently implemented forces for U_projection_type = 'ortho-atomic’, I have come across several cases where the geometry optimizations would fail: the energy would fluctuate without convergence using BFGS, and even damped dynamics would frequently not converge. Such a behavior usually points to a problem with the forces, which are either not exact or noisy.

Based on this observation, I decided to perform a routine test: I compute the integral of the forces times displacement along a path in configurational space, and compare it with the energy difference between the initial and final point. If the forces are exact derivatives of the energy, the path integral must be identical to the energy difference (up to machine precision if the step size of the path integral goes to zero). The method is described in sec 3.3, https://doi.org/10.1016/j.cpc.2020.107415, see also the attachment of the relevant section.



I performed this test for a TiO2 Rutile structure with 6 atoms per cellI, along a path in a random direction with 0.001 Ang step size and with a total of 501 steps. I used 3 different settings: DFT without U, DFT+U with U_projection_type = ‘atomic’ (U=4.7147), and DFT+U with U_projection_type = 'ortho-atomic’ (U=4.7147). In the case of No-U and ‘atomic’, the relative error of the path integral and energy difference is less than1.d-4%, as expected when the forces are correct. But for 'ortho-atomic’ the relative error is a whopping 8%! This indicates that the forces are not the exact derivatives of the energy for U_projection_type = 'ortho-atomic’ (see also plot below, note that the y scale is logarithmic for the last two subplots).


It seems that there is a bug in the implementation of the forces for U_projection_type = 'ortho-atomic’, or that some terms have been neglected when computing the derivatives.

In the long run, it might be helpful to include such tests routinely as part of the QE development process.

Best
Max
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20200912/b9782bdf/attachment.html>