<div dir="ltr">Hi All,<div><br></div><div>I am nearing the point where I have a working many body dispersion code with gradients ready to be checked back into QE. I've followed the template of the TS module for how to hook into the rest of QE, but unfortunately the cell stresses haven't yet been implemented as far as I can tell. I'm going to ask about TS because once I have those hooks working, I should be able to replicate these for MBD without issue.</div><div><br></div><div>One thing that I'm struggling with is how to use the stress hooks. I have a working and verified energy expression as well as gradients with respect to lattice vectors and nuclear positions (just like TS). Just for clarity, I have:</div><div><br></div><div>E, -dE/dR, and dE/dh_ij</div><div><br></div><div>And now I want to take dE/dh_ij and compute stresses with it. There are known ways to do this in the literature (eqn A6 of <a href="http://crm2.univ-lorraine.fr/pages_perso/Angyan/papers/jcp05v122p124508.pdf">http://crm2.univ-lorraine.fr/pages_perso/Angyan/papers/jcp05v122p124508.pdf</a>). For TS, i've tried implementing this via:</div><div><div>    do i=1,3,1</div><div>      do j=1,3,1</div><div>        sigma_ts(i, j) = HtsvdW(1, i)*h(1, j) + HtsvdW(2, i)*h(2, j) + HtsvdW(3, i)*h(3, j)</div><div>        end do</div><div>      end do</div><div>    sigma_ts = sigma_ts/omega</div></div><div>sigma = sigma + 2.0_DP*sigma_ts</div><div><br></div><div>and by computing sigma_ts via MATMUL( HtsvdW, TRANSPOSE( h ) ) / omega (line 761 of vofrho.f90 in cpv) and neither worked. In fact, on a benzene crystal, previous work has shown that the TS optimized cell volume should be ~465ang^3, but with the stress implementation above, the TS optimization gives ~600ang^3. Something tells me that something is wrong here but I'm not sure what.</div><div><br></div><div>Note, the formula in A6 expects the lattice vectors to be in the columns of h, but QE keeps them in the rows. Also, I've multiplied by the factor of 2 to convert from hartree to ry. This is why I have h(1, j) instead of h(j, 1). I expect that I should get the same cell parameters from a vc-relax job as are in the literature with numerical gradients, but i don't unfortunately. Also, swapping the sign on sigma_ts didn't help.</div><div><br></div><div>Anyway, any ideas of what to do wrong? Perhaps my convention is incorrect? Perhaps units?</div><div><br></div><div>Best,</div><div>Thomas</div></div>