Investigate forces test failure

[antoine-levitt]

Eg https://github.com/JuliaMolSim/DFTK.jl/runs/5439547530?check_suite_focus=true with println(abs(diff_findiff - diff_forces)) on the first forces test I can't get it to go down far below 1e-5 with symmetry. Without symmetry everything is fine and I get 1e-10. Possibly linked with SYMMETRY_TOLERANCE. We should investigate and set the tolerance far below 1e-5.

[antoine-levitt]

Decreasing SYMMETRY_TOLERANCE does not change anything

[antoine-levitt]

Reproducer

# Very basic setup, useful for testing
using DFTK

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [[1.01, 1.02, 1.03]/8, -ones(3)/8]]

Ecut = 7
kgrid = [2, 2, 2]
kshift = [0, 0, 0]
is_converged = DFTK.ScfConvergenceDensity(1e-12)

model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model; Ecut, kgrid, kshift)
scfres = self_consistent_field(basis; is_converged)
F1 = compute_forces(scfres)

model = model_LDA(lattice, atoms, symmetries=false)
basis = PlaneWaveBasis(model; Ecut, kgrid, kshift)
scfres = self_consistent_field(basis; is_converged)
F2 = compute_forces(scfres)

println(norm(F1-F2))

The error goes down when using larger supersampling factors, and the error is zero when using a model without xc

[antoine-levitt]

It's not the fault of lowpass_for_symmetry! (which actually is probably completely useless and should be removed...)

[antoine-levitt]

A simpler test: run a computation without symmetry, and then reload it in a symmetric computation and compare the density to the symmetrized one. There's a similar difference. This is because the xc term depends on the actual fft_size, which itself is not preserved by symmetries. The Hartree is not affected because it is exact (does not depend on fft_size as long as supersampling >= 2), and the variational basis sets for the psi are symmetry-invariant. Still not completely clear to me if there's something we can do or if we should just write it up in the docs and live with it.

[jaemolihm]

Do you mean that there are symmetry operations (fractional translations) that are not compatible with the real space grid ?

[mfherbst]

Depends which grid one is talking about. In DFTK we have the Ecut and the fft_size.

The Ecut determines the "sphere" of plane waves used for each k-Point to represent the orbitals, the fft_size the grid on which FFTs are done and quantities like the potential and density are represented. The fft_size is by default chosen such that |ψ|^2 and the convolution of the density with the Hartree potential can be represented exactly. Therefore if the k-Points and the respective Ecut-grids respect the symmetry, so does the Hartree term and its action. This is what we enforce in our symmetry setup.

However, the XC potential is a more complicated function of the density (even LDA is ρ^{1/3}), which for a ψ represented on an Ecut-basis cannot be represented exactly on any truncated plane-wave grid. Now the error one makes in ρ^{1/3} is generally not evenly distributed amongst the symmetry operations (e.g. simple inversion is usually fine, but rotatations could rotate a G-vector out of the considered grid). Therefore even if ρ is perfectly symmetric ρ^{1/3} and thus the XC potential may not be.

[antoine-levitt]

The operation that doesn't preserve the symmetries is to compute the xc energy (and therefore potential) from a rho given by its Fourier coefficients. I think translations and rotations both contribute to the symmetry breaking, but that should be easy to check (can't do it now though) : on the mwe above, the translations should be preserved for some values of fft grid.

[antoine-levitt]

But I'm pretty sure we cannt do anything about it. @mfherbst is it easy for you to run the above case in eg QE with and without symmetries and see if the forces differ?

[jaemolihm]

Sorry for interrupting again, but I was curious about this, as it seemed a fundamental problem that should have been known.

atoms = [Si => [[1.01, 1.02, 1.03]/8, -ones(3)/8]] has symmetry "inversion + fractional translation by [1/800, 2/800, 3/800]". This nontrivial fractional translation does not preserve the real-space grid (unless there are 800^3 grid points...). In QE, such a symmetry is not used (maybe not detected at all - QE does not use spglib).

I expect that atoms = [Si => [[1.005, 1.01, 1.015]/8, -[1.005, 1.01, 1.015]/8] would give the same forces with and without symmetry in DFTK, because its symmetry is "inversion", without any fractional translation. In this case, QE detects the inversion symmetry and uses it.

My expectation as a user (may not be universal) is that symmetries that break fft_grid should be discarded at the basis level. (This is what QE does.)

 &control
    calculation  = 'scf'
    prefix       = 'si'
    restart_mode = 'from_scratch'
    pseudo_dir   = '.'
    outdir       = './temp/'
    tprnfor      = .true.
    tstress      = .true.
    verbosity    = "high"
 /
 &system
    ibrav           = 0
    nat             = 2
    ntyp            = 1
    ecutwfc         = 14
    !nosym = .true.
    !noinv = .true.
 /
 &electrons
    diagonalization = 'david'
    mixing_beta     = 0.7
    conv_thr        = 1.0d-15
 /
ATOMIC_SPECIES
  Si  1.  Si.dojo.sr.ONCVv0.4.lda.standard.upf
CELL_PARAMETERS bohr
  0.00  5.13  5.13
  5.13  0.00  5.13
  5.13  5.13  0.00
K_POINTS automatic
 2 2 2 0 0 0
ATOMIC_POSITIONS crystal
  Si   -0.125              -0.125              -0.125
  Si    0.125+0.000625*2    0.125+0.000625*4    0.125+0.000625*6
!  Si   -0.125-0.000625    -0.125-0.000625*2   -0.125-0.000625*3
!  Si    0.125+0.000625     0.125+0.000625*2    0.125+0.000625*3

[antoine-levitt]

Thanks! Hm, discarding in the basis the symmetries that break the discrete grid indeed seems like a reasonable thing to do: that way, hopefully, symmetry handling becomes a pure optimization (at least at SCF convergence), which would be nice. The rotations can also be checked easily, I guess. I'll try to implement it and see how it goes, but my bandwidth on this is about to be severely limited so it'll have to wait!

(I hate symmetries)

[antoine-levitt]

Btw this prediction is correct:

I expect that atoms = [Si => [[1.005, 1.01, 1.015]/8, -[1.005, 1.01, 1.015]/8] would give the same forces with and without symmetry in DFTK, because its symmetry is "inversion", without any fractional translation. In this case, QE detects the inversion symmetry and uses it.

So OK, filtering out these symmetries sounds very reasonable.

[antoine-levitt]

Turns out this is super annoying because it means bzmesh_ir_wedge needs to know fft_size, so we need to overhaul the pwbasis constructors before doing this.