Refactor and extend tests for forces (#1080)

Author: mfherbst

Project.toml

@@ -1,7 +1,7 @@
 name = "DFTK"
 uuid = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
 authors = ["Michael F. Herbst <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.7.10"
+version = "0.7.11"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/terms/magnetic.jl

@@ -10,6 +10,18 @@ struct Magnetic
 end
 (M::Magnetic)(basis) = TermMagnetic(basis, M.Afunction)
 
+struct MagneticFromValues
+    # Apotential[α] is an array of size fft_size for α=1:3
+    Apotential::Vector{<:AbstractArray}
+end
+function (M::MagneticFromValues)(Apotential)
+    @assert length(Apotential) == 3
+    @assert size(Apotential[1]) == basis.fft_size
+    @assert size(Apotential[2]) == basis.fft_size
+    @assert size(Apotential[3]) == basis.fft_size
+    TermMagnetic(basis, Apotential)
+end
+
 struct TermMagnetic <: Term
     # Apotential[α] is an array of size fft_size for α=1:3
     Apotential::Vector{<:AbstractArray}

test/forces.jl

@@ -86,7 +86,7 @@ end
     function compute_properties(lattice)
         model_strained = Model(model; lattice)
         basis = PlaneWaveBasis(model_strained; Ecut, kgrid)
-        scfres = self_consistent_field(basis; tol=1e-10, nbandsalg)
+        scfres = self_consistent_field(basis; tol=1e-11, nbandsalg)
         ComponentArray(
            eigenvalues=stack([ev[1:10] for ev in scfres.eigenvalues]),
            ρ=scfres.ρ,
@@ -97,9 +97,11 @@ end
     end
 
     strain_isotropic(ε) = (1 + ε) * model.lattice
-    strain_anisotropic(ε) = DFTK.voigt_strain_to_full([ε, 0., 0., 0., 0., 0.]) * model.lattice
+    function strain_anisotropic(ε)
+        DFTK.voigt_strain_to_full([ε, 0., 0., 0., 0., 0.]) * model.lattice
+    end
 
-    @testset "$strain_fn" for strain_fn in [strain_isotropic, strain_anisotropic]
+    @testset "$strain_fn" for strain_fn in (strain_isotropic, strain_anisotropic)
         f(ε) = compute_properties(strain_fn(ε))
         dx = ForwardDiff.derivative(f, 0.)
 
@@ -109,7 +111,7 @@ end
         dx_findiff = (x2 - x1) / 2h
         @test norm(dx.ρ - dx_findiff.ρ) * sqrt(basis.dvol) < 1e-6
         @test maximum(abs, dx.eigenvalues - dx_findiff.eigenvalues) < 1e-6
-        @test maximum(abs, dx.energies - dx_findiff.energies) < 1e-5
+        @test maximum(abs, dx.energies - dx_findiff.energies) < 3e-5
         @test dx.εF - dx_findiff.εF < 1e-6
         @test maximum(abs, dx.occupation - dx_findiff.occupation) < 1e-4
     end

test/forwarddiff.jl

@@ -0,0 +1,4 @@
+2
+pbc=[T, T, T] Lattice="2.87000000 0.00000000 0.00000000 0.00000000 2.87000000 0.00000000 0.00000000 0.00000000 2.87000000" Properties=species:S:1:pos:R:3:velocities:R:3:Z:I:1
+Fe         0.10000000      -0.03000000       0.05000000         0.00000000       0.00000000       0.00000000         26
+Fe         1.43500000       1.43500000       1.43500000         0.00000000       0.00000000       0.00000000         26

test/structures/Fe_rattled.extxyz

@@ -0,0 +1,8 @@
+6
+name=model_0 pbc=[T, T, T] Lattice="4.48609000 0.00000000 0.00000000 0.00000000 4.48609000 0.00000000 0.00000000 0.00000000 2.91186000" Properties=species:S:1:pos:R:3:velocities:R:3:Z:I:1
+Ge         0.00629788       0.00222026       0.01567598         0.00000000       0.00000000       0.00000000         32
+Ge         2.23696804       2.24090229       1.45515035         0.00000000       0.00000000       0.00000000         32
+O         1.38588856       1.40290802      -0.00302997         0.00000000       0.00000000       0.00000000          8
+O         3.09582248       3.09525318       0.00021341         0.00000000       0.00000000       0.00000000          8
+O         0.85224093       3.64776914       1.45805119         0.00000000       0.00000000       0.00000000          8
+O         3.62813145       0.85721112       1.45468113         0.00000000       0.00000000       0.00000000          8

test/structures/GeO2_rattled.extxyz

@@ -0,0 +1,8 @@
+6
+name=model_0 pbc=[T, T, T] Lattice="4.48609000 0.00000000 0.00000000 0.00000000 4.48609000 0.00000000 0.00000000 0.00000000 2.91186000" Properties=species:S:1:pos:R:3:charge:R:1:occupancy:R:1:velocities:R:3:vdw_radius:R:1:covalent_radius:R:1:Z:I:1
+Ge         0.10000000       0.00000000       0.00000000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         2.10000000         1.22000000         32
+Ge         2.24304500       2.22304500       1.45593000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         2.10000000         1.22000000         32
+O         1.39068790       1.39068790       0.00000000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         1.55000000         0.73000000          8
+O         3.09540210       3.08540210       0.00000000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         1.55000000         0.73000000          8
+O         0.75235710       3.53373290       1.45593000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         1.55000000         0.73000000          8
+O         3.73373290       0.85235710       1.45593000         0.00000000         1.00000000         0.00000000       0.00000000       0.00000000         1.55000000         0.73000000          8