Symmetrize forces in DFTKCalculator

[Technici4n]

This fix makes for example BaTiO3 geometry optimization converge! Without it, the symmetries get lost at some point in the optimization step.

[mfherbst]

I actually prefer this. I think having a calculator interface that behaves fundamentally different compared to when using DFTK's API directly is unexpected behaviour and bound to lead to subtle inconsistencies down the road. If we find this needs to be done explicitly to ensure proper behaviour than we should always do it.

• Does it need to be the model symmetries ? Are the basis symmetries not sufficient ? Because right now we make sure the forces of each term are symmetric with the basis symmetries only. This is indeed only explicitly done in the nonlocal term, but the others should be symmetric out of the box. but maybe there are numerical issues and we need an explicit symmetrisation wrt. basis symmetries as well. Did you test that explicit symmetrisation wrt. basis symmetries is not yet sufficient ?
• How expensive is this symmetrisation anyways ? Do we even care if this is done every time or not ?
• In contrast to using the basis symmetries, symmetrising with the Model symmetries changes behaviour and I am not convinced that this is the right thing in all cases, so if we switch the default behaviour to using the model symmetries, there should be a kwarg (e.g. symmetries ?) to enable an easy change to using the basis symmetries instead. This subtlety than has to be documented with an appropriate sentence in the docstring of compute_forces and compute_forces_cart as well as clear comments in the code:
  • That his was needed to get geometry optimisation to work
  • how one can switch from one kind of symmetries to the other and what this implies (once the symmetries of the physical problem encoded in Model, once only those symmetries that the discretisation can represent)

[Technici4n]

I put it here because the DFTKCalculator reuses the previous model's symmetries. That's also why I used the model's symmetries specifically. IMO for geometry optimization you don't want to lose any model symmetry.

[Technici4n]

Quick check with silicon:

# Very basic setup, useful for testing
using DFTK
using PseudoPotentialData

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf"))
atoms     = [Si, Si]
positions = [ones(3)/8 + [0.001, 0.001, 0.001], -ones(3)/8]

model = model_DFT(lattice, atoms, positions; functionals=LDA())
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])
scfres = self_consistent_field(basis, tol=1e-8)

Fred = compute_forces(scfres)
Fred_sym_basis = DFTK.symmetrize_forces(basis, Fred)
Fred_sym_model = DFTK.symmetrize_forces(model, Fred; model.symmetries)

let
    println("Raw reduced forces: ", Fred)
    println("  Total: ", sum(Fred))
    println("Basis-symmetrized forces: ", Fred_sym_basis)
    println("  Total: ", sum(Fred_sym_basis))
    println("Model-symmetrized forces: ", Fred_sym_model)
    println("  Total: ", sum(Fred_sym_model))
end

This gives

Raw reduced forces: StaticArraysCore.SVector{3, Float64}[[-0.014528303367595304, -0.014528303367597068, -0.014528303367596604], [0.014528291138446968, 0.014528291138447471, 0.014528291138445375]]
  Total: [-1.2229148336345697e-8, -1.2229149597489664e-8, -1.2229151228129731e-8]
Basis-symmetrized forces: StaticArraysCore.SVector{3, Float64}[[-0.014528303367596326, -0.014528303367596326, -0.014528303367596326], [0.014528291138446605, 0.014528291138446605, 0.014528291138446607]]
  Total: [-1.222914972065503e-8, -1.222914972065503e-8, -1.2229149718920307e-8]
Model-symmetrized forces: StaticArraysCore.SVector{3, Float64}[[-0.014528297253021466, -0.014528297253021466, -0.014528297253021466], [0.014528297253021466, 0.014528297253021466, 0.014528297253021466]]
  Total: [0.0, 0.0, 0.0]

So it looks like the forces coming out of compute_forces more or less already satisfy the basis' symmetries, but they do not satisfy the model's symmetries. Confirmed by printing norm(Fred_sym_basis - Fred) / norm(Fred_sym_basis) which is 5e-14.

[Technici4n]

So my conclusion is:

• compute_forces seems already basis-symmetrized
• geometry optimization needs model-symmetrization

[antoine-levitt]

Can you explain a bit more the behavior of the optimizer? IMO it's reasonable that the forces are not model symmetrized, no?

[Technici4n]

Can you explain a bit more the behavior of the optimizer? IMO it's reasonable that the forces are not model symmetrized, no?

Without model symmetrization, the optimizer step (gradient based) can break the symmetries. Combined with how DFTKCalculator reuses the symmetries from the previous model this can lead to gradients that cannot be zeroed by the optimizer. Disabling symmetries entirely is also not satisfactory (ignoring the slowness aspect) since the optimizer can then converge to a non-symmetric equilibrium geometry.

Generally speaking if you have a symmetric structure you have two choices:

• Accept that the symmetry might be lost.
• Force the optimizer to remain on the symmetric subspace.

For my particular case I really do not want to lose the symmetries. IMO this is what you generally want when relaxing a structure with some known symmetries.

[antoine-levitt]

That's not the right fix I think. If you actually want to preserve the symmetries, do it in the optimizer (ie symmetrize the positions). Otherwise the force is not the derivative of the energy, which I think we really don't want.

[Technici4n]

Indeed this might be better fixed at the geometry optimization level. However, it currently does not know the concept of symmetries and I don't see it knowing about it in the near future either. Fixing this is essential for any serious geometry optimization run, and this PR is IMO a clear improvement in that regard. I have this fix locally, it works, and it allowed me to move on with my actual research; clearly good enough for now IMO. There is a TODO in case someone wants to investigate further.

[mfherbst]

Here is my version of this. Essentially geometry optimisation automatically enables this symmetrisation, else its off.

[Technici4n]

LGTM