Forces and stresses

[antoine-levitt]

Self-explanatory :-)

• [x] forces
• [ ] stresses

[cortner]

+1

[antoine-levitt]

Forces is https://github.com/JuliaMolSim/DFTK.jl/pull/106 (modulo https://github.com/JuliaMolSim/DFTK.jl/issues/105). Stresses I'd need to take a look at the formulas.

[cortner]

@hjchen1983 do you have a link or reference to compute the virial in DFT?

[antoine-levitt]

It's in R. Martin Electronic Structure appendix G. It's pretty straightforward but annoying because the lattice pops up in every energy expression through the reciprocal lattice vectors. OTOH this is typically the kind of expressions that ForwardDiff should like (it's only 9 perturbations, and not very costly), so we should think of doing that once the interface for the terms is stabilized.

[cortner]

You should be able to reduce it to 6 perturbations. If I remember correctly The virial is - up to rescaling - the Cauchy stress tensor, hence symmetric.

Would indeed be interesting to see if you can do this with ForwardDiff, but how will it exploit the the HF-theorem?

[antoine-levitt]

Would indeed be interesting to see if you can do this with ForwardDiff, but how will it exploit the the HF-theorem?

That you do yourself. It's partially automated differentiation: you tell it how to compute <psi, H(params) psi> for given psi, and you ask it to compute the derivative of that wrt the parameters. Figuring out if it's possible to AD through the whole solver (given some help), including the derivatives of the psi, is a lot more interesting, but harder.

We could also do that "HF AD" for the forces but to get the correct asymptotic scaling we'd need reverse diff which still sucks, so I just went ahead and coded it up explicitly.

[antoine-levitt]

You should be able to reduce it to 6 perturbations. If I remember correctly The virial is - up to rescaling - the Cauchy stress tensor, hence symmetric.

For future reference, the argument is: the energy E(A) where A is the matrix of lattice vectors is invariant wrt global rotations, so E(RA)=E(A) for all rotations R, so differentiating at R=I we get dE . (MA) = 0 when M is an antisymmetric matrix. It's enough to test the differential dE on perturbations of the form SA where S is symmetric - hence 6 dof.

[mfherbst]

If I get it right the idea is that R = e^M with M being the antisymmetric matrix, so that 0 = dE(e^M A) if only M varies, which yields dE . (MA) = 0 at M = 0 right?

[antoine-levitt]

Yup. Pedanticly, if your energy is invariant wrt a Lie group, the associated Lie algebra is in the kernel of the differential of the energy. I learned the other day that theoretical physicists call that with the funny name of Nambu-Goldstone modes (https://en.wikipedia.org/wiki/Goldstone_boson)

[mfherbst]

Stresses will come after 0.1.0.

[antoine-levitt]

This is still missing a nice API on top of it, but is basically done (via ForwardDiff) in https://github.com/JuliaMolSim/DFTK.jl/pull/476

[antoine-levitt]

@cortner I'm very confused by the expressions I see in the literature, can you perhaps shed some light on the proper definition of stresses? What we compute is E(lattice), the energy per unit cell. Is the stress sigma = - 1/det(lattice_0) dE/dlattice, taken at lattice_0? Or is the stress in "rescaled" coordinates, ie something like sigma = -1/det(lattice_0) dE/d(lattice_0 + M*lattice_0), taken at M=0 (which would make sigma symmetric)? I see these weird symmetrizations, is this just an optimization or will I get wrong results if I don't symmetrize explicitly? Any good source for this?

[cortner]

its been a while, but i'll remind myself and report back later. I\m fairly certain the virial is symmetric, but I need to remind myself why.

[antoine-levitt]

The definition sigma = -1/det(lattice_0) dE/d(lattice_0 + M*lattice_0) is symmetric because E(R * lattice) = E(lattice) for all rotation matrices R. The thing I'm not sure about is what is the actual definition of stresses, ie what coordinate system people use.

[cortner]

What we compute is E(lattice), the energy per unit cell. Is the stress sigma = - 1/det(lattice_0) dE/dlattice, taken at lattice_0? Or is the stress in "rescaled" coordinates, ie something like sigma = -1/det(lattice_0) dE/d(lattice_0 + M*lattice_0), taken at M=0 (which would make sigma symmetric)?

The virial / stress is defined for any configuration, it need not be a lattice (unless you equate periodic with lattice which I guess makes mathematical sense...) You are given a cell C and positions rj. As you deform the cell vectors c_a -> F c_a you deform the positions in the same way rj -> F rj. The virial is (minus) the derivative of E w.r.t. F at F = Id. If I may rewrite your expression as

   d / dM E(lattice_0 + M*lattice_0)

then you would evaluate this at M = 0 since F = I + M.

The rescaling with det is then called the stress, and is defined as

  stress = - virial / volume

So I think you have a sign error in your expression.

I double-checked this against my codes and this is exactly what I've implemented.

julia> virial(emt, at)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.0389347  -0.636515   -0.0149505
 -0.636515    0.0778579  -0.117448
 -0.0149505  -0.117448   -0.177202

This symmetry should come from exactly what you said above - i.e. dot with an infinitesimal roation gives 0.

I can't find a reference for all this right now, but I believe I tested the implementation against ASE so hopefully my expressions are correct.

[antoine-levitt]

Thank you, that's very helpful! That makes sense. Can you clarify the signs? You say virial is minus derivative of energy, then you write stress = - virial, meaning it's plus derivative of energy? I would have though virial = derivative, stress = -virial/vol = -derivative/vol?

We call it lattice in the code because it's the matrix of lattice vectors and defines the unit cell, but you're right it's not necessarily a lattice.

[cortner]

no virial is like a force i.e. - derivative of E

[cortner]

one of those things that I just don't question ...

[antoine-levitt]

OK, thanks! Martin says pressure = -1/3 tr(stress). So for a typical quadratic energy-volume curve, if you compress the material, it has dE/dvol negative, so negative stress, so positive pressure. That makes sense!

[mfherbst]

Also about the sign: I just read some comment in the ASE sources that there are both sign conventions used ...

[antoine-levitt]

Done in https://github.com/JuliaMolSim/DFTK.jl/pull/487