Fix symmetrization of stresses

[jaemolihm]

As @antoine-levitt spotted in #511 (comment), there were several issues in my implementation of symmetrize_stresses.

1. since both the stress tensor and the symmetry operation are in reduced coordinates, why does the lattice even appear?
2. why do you use S (which is in reciprocal space) rather than W (which is in real space)? Did you copy the formula from somewhere or did you derive it yourself?

The problems are:

• The implementation was wrong in that a transpose is incorrectly applied to lattice.
• Test did not fail because lattice was symmetric. I updated the test to use non-symmetric lattice. Test fails with the old implementation but passes with the new one.
• S_reduced in the previous version was actually W_cart (Regarding point 2, I was using W indeed, but incorrectly referred it to S.)

One thing that still needs to be discussed is the definition of stress tensor and point 1 ("why does the lattice even appear").
In

DFTK.jl/src/postprocess/stresses.jl

Line 33 in 70fc3c8

stresses_not_symmetrized = ForwardDiff.gradient(M -> HF_energy((I+M) * L), zero(L)) / Ω

, the lattice is transformed as L=(I+M) * L, so M acts on the columns of L which are Cartesian axis. In this case, the first lattice vector transform as L'[i, 1] = L[i, 1] + sum_j M[i, j] * L[j, 1].
But the formula I found in the literature (Appendix G of Richard Martin's book) is that the lattice vector transform as L'[i, 1] = L[i, 1] + sum_b M_Martin[a, b] * L[i, b], i.e. L'_Martin = L * (I + M_Martin^T).
So the relation between M and M_Martin is M = L * M_Martin^T * L^{-1}. I think this is the reason why lattice appears in the symmetrization routine.

[antoine-levitt]

Thank you for looking at this again! I knew I shouldn't have been lazy at the time and checked it carefully :) I think the main difference is that our stress is in Cartesian (our M is a transformation of R3 of each lattice vector expressed in the canonical basis) and Martin's is in reduced a transformation of R3 expressed in the lattice basis. The transpose is weird though. I'll take a close look and add some comments to have it right.

[jaemolihm]

Right, I think there was a comment that stresses was in reduced coordinates which confused me. If you figure out a solution (either converting to reduced, or using Cartesian), I can fix those in this PR.

[antoine-levitt]

Where do you find the deformation in appendix G of Martin?

[jaemolihm]

Sorry, I checked again and it was just my confusion... 😓 The definition in Martin is consistent with what is done in DFTK.

[antoine-levitt]

OK I think I get it. Our definition of stresses is in cartesian coordinates. (A.28) of the QE paper is given in cartesian coordinates (they use the fact that inv(R) = R', which is only valid in cartesian). In the PR you pass the stress to reduced coordinates, symmetrize and pass back. Would it not be clearer to pass the symmetry to cartesian coordinates instead, and directly apply (A.28)?

So there's still an API issue. Our current compute_stresses is very much in cartesian coordinates. Should we mimic the forces interface, and have a compute_stresses_cart (which does what is currently done) and a compute_stresses (which computes them in reduced coordinates)? Or is it too confusing to make the distinction and when people talk about stresses they always mean in cartesian coordinates? I could see stresses in reduced coordinates being pretty useful (eg to identify the shears and whatnot), but I don't use these in practice so I don't know.

(also minor style point, A / B is slightly better numerically than A * inv(B))

[antoine-levitt]

OK, let's just rename compute_stresses to compute_stresses_cart for uniformity with the forces, and leave compute_stresses undefined (perhaps with a @deprecate compute_stresses(args...) compute_stresses_cart(args...)). @jaemolihm will you have time to finish this? It's making me a bit nervous to have such a bug in the code right now. Happy to take it over if you're busy.

[jaemolihm]

@antoine-levitt sorry for the delay. I think I wont have the time for a few days, so It would be great if you could take this over.

[antoine-levitt]

OK! I opened a new PR with some added functions to hopefully get this right in the future, would be great if you could review - if not no worries, the code is functionally identical to this PR and so should be OK.

[mfherbst]

Thanks for your efforts on this, @jaemolihm and @antoine-levitt. I'm very happy that we are close to (hopefully) closing this chapter again for some time. Once #603 is in, I'll tag a new version with all polished up symmetry routines!

[antoine-levitt]

I'm very happy that we are close to (hopefully) closing this chapter again for some time

Where have I heard this before... :-) But yeah, I understand this bit much better now, thanks @jaemolihm for all your help!