[Q-e-developers] Stress Hooks in QE?

[Thomas Markovich]

Hi All,

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.

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:

E, -dE/dR, and dE/dh_ij

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
http://crm2.univ-lorraine.fr/pages_perso/Angyan/papers/jcp05v122p124508.pdf).
For TS, i've tried implementing this via:
    do i=1,3,1
      do j=1,3,1
        sigma_ts(i, j) = HtsvdW(1, i)*h(1, j) + HtsvdW(2, i)*h(2, j) +
HtsvdW(3, i)*h(3, j)
        end do
      end do
    sigma_ts = sigma_ts/omega
sigma = sigma + 2.0_DP*sigma_ts

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.

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.

Anyway, any ideas of what to do wrong? Perhaps my convention is incorrect?
Perhaps units?

Best,
Thomas
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20150325/ba026aa0/attachment.html>