<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">Dear Iurii,<div class="">Thanks for your response. Actually, I made a mistake in my previous mail. The energies are in eV, not Ry! But the relative errors are correct.<br class=""><div><br class=""><blockquote type="cite" class=""><div class="">On 12 Sep 2020, at 15:32, Timrov Iurii <<a href="mailto:iurii.timrov@epfl.ch" class="">iurii.timrov@epfl.ch</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class="">Dear Max,</div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><br class=""></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">Thanks for performing these tests!</div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><br class=""></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">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:</div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><a href="https://en.wikipedia.org/wiki/Numerical_differentiation" class="OWAAutoLink" id="LPlnk565393" previewremoved="true">https://en.wikipedia.org/wiki/Numerical_differentiation</a><br class=""></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">For example, for NiO I obtained the following results:</div><div style="margin-top: 0px; margin-bottom: 0px;" class="">Analytical: F = -0.048891 Ry/Bohr</div><div style="margin-top: 0px; margin-bottom: 0px;" class="">Numerical: <span style="font-size: 12pt;" class="">F = -0.048891 Ry/Bohr</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">As you can see forces agree up to <span style="font-size: 12pt;" class="">10^{-6} Ry/Bohr. In a similar way I tested also stress.</span></div></div></div></blockquote>Is the relative error of the same magnitude as the ones form “atomic” or no U at all?<br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="caret-color: rgb(33, 33, 33); color: rgb(33, 33, 33); font-family: wf_segoe-ui_normal, "Segoe UI", "Segoe WP", Tahoma, Arial, sans-serif, serif, EmojiFont; font-size: 15px;" class="">> 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.</span><span style="font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">The analytical formulas are<span class="Apple-converted-space"> </span><b class="">exact</b>, 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 </font>agreement<font size="3" class=""> between analytical and numerical derivatives for NiO and a couple of other systems.</font></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class=""><br class=""></font></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">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. </font></div></div></div></blockquote>What you are suggesting is the central difference. For a particular structure this might give the correct force, but if you do the integral for a completely arbitrary path (and many of them) you can uncover particular configurations where the forces are wrong (e.g., when atoms cross the periodic boundaries), or you accumulate small errors in the forces that show up as significant deviations from the energy difference. On several occasions, I have discovered in this way several bugs in my own implementations of forces that were not easily detectable, and I found that this kind of test is rather robust. <br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class=""><br class=""></font></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">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? </font></div></div></div></blockquote>I am using the new version with your patch. I am running it on daint, compiled with Intel version 19.1.1.217 and 2020.1.217 MKL libraries.<br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">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).</font></div></div></div></blockquote>I can increase the cutoff and rerun the test. But I am already using very tight parameters (100.00 Ry). I could increase the k-point mesh (which is not fully converged), but then again I would expect that I would see larger errors for No-U and “atomic”.<br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class=""> Also, maybe it would useful to try other compilers and libraries.</font></div></div></div></blockquote>I can recompile with PGI and see if I see a change in behavior.<br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">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? </font></div></div></div></blockquote>Yes, I will perform the test with NiO. Please send me your input files.<br class=""><blockquote type="cite" class=""><div class=""><div id="divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">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?</font></div></div></div></blockquote>I have attached the atomic positions for every displacement along the path in ascii format. The PW parameters are given below:</div><div><br class=""></div><div> &CONTROL<br class=""> restart_mode = 'from_scratch' ,<br class=""> outdir='./' ,<br class=""> pseudo_dir = './' ,<br class="">calculation = "scf" ,<br class="">tstress = .true. ,<br class="">tprnfor = .true. ,<br class="">/<br class=""> &SYSTEM<br class="">U_projection_type = 'ortho-atomic'<br class=""> ibrav= 0, <br class=""> nat= 6, <br class=""> ntyp = 2<br class=""> nosym=.false.<br class=""> ecutwfc = 100.0 ,<br class=""> ecutrho = 1200.0 ,<br class=""> occupations = 'smearing' ,<br class=""> smearing = 'gaussian' ,<br class=""> degauss = 0.01 ,<br class=""> lda_plus_u=.true.<br class=""> lda_plus_u_kind = 0,<br class=""> Hubbard_U(1) = 4.7147<br class="">nat = 6 ,<br class="">ntyp = 2 ,<br class="">nosym = .true. ,<br class="">/<br class=""> &ELECTRONS<br class=""> conv_thr = 1.0d-10 ,<br class="">/<br class="">ATOMIC_SPECIES<br class="">Ti 2.00000 Ti.PSP<br class=""> O 2.00000 O.PSP<br class="">K_POINTS automatic<br class=""> 4 4 6 0 0 0</div><div></div></div></body></html>