Forces

src/DFTK.jl

@@ -153,4 +153,7 @@ export load_density
 export load_atoms
 include("external/etsf_nanoquanta.jl")
 
+export forces
+include("forces.jl")
+
 end # module DFTK

src/energy_ewald.jl

@@ -1,4 +1,5 @@
 using SpecialFunctions: erfc
+using ForwardDiff
 
 """
     energy_ewald(lattice, [recip_lattice, ]charges, positions; η=nothing)
@@ -8,9 +9,10 @@ charges in a uniform background of compensating charge to yield net
 neutrality. the `lattice` and `recip_lattice` should contain the
 lattice and reciprocal lattice vectors as columns. `charges` and
 `positions` are the point charges and their positions (as an array of
-arrays) in fractional coordinates.
+arrays) in fractional coordinates. If `forces` is not nothing, minus the derivatives
+of the energy with respect to `positions` is computed.
 """
-function energy_ewald(lattice, charges, positions; η=nothing)
+function energy_ewald(lattice, charges, positions; η=nothing, forces=nothing)
     T = eltype(lattice)
 
     for i=1:3
@@ -21,20 +23,23 @@ function energy_ewald(lattice, charges, positions; η=nothing)
             return T(0)
         end
     end
-    energy_ewald(lattice, T(2π) * inv(lattice'), charges, positions, η=η)
+    energy_ewald(lattice, T(2π) * inv(lattice'), charges, positions; η=η, forces=forces)
 end
-function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
+function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, forces=nothing)
     T = eltype(lattice)
     @assert T == eltype(recip_lattice)
-    @assert(length(charges) == length(positions),
-            "Length of charges and positions does not match")
+    @assert length(charges) == length(positions)
     if η === nothing
         # Balance between reciprocal summation and real-space summation
         # with a slight bias towards reciprocal summation
         η = sqrt(sqrt(T(1.69) * norm(recip_lattice ./ 2T(π)) / norm(lattice))) / 2
     end
+    if forces !== nothing
+        @assert size(forces) == size(positions)
+        forces_real = copy(forces)
+        forces_recip = copy(forces)
+    end
 
-    #
     # Numerical cutoffs
     #
     # The largest argument to the exp(-x) function to obtain a numerically
@@ -89,11 +94,24 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
 
             any_term_contributes = true
             sum_recip += sum_strucfac * exp(-exponent) / Gsq
+
+            if forces !== nothing
+                for (ir, r) in enumerate(positions)
+                    Z = charges[ir]
+                    dc = -Z*2T(π)*G*sin(2T(π) * dot(r, G))
+                    ds = +Z*2T(π)*G*cos(2T(π) * dot(r, G))
+                    dsum = 2cos_strucfac*dc + 2sin_strucfac*ds
+                    forces_recip[ir] -= dsum * exp(-exponent)/Gsq
+                end
+            end
         end
         gsh += 1
     end
     # Amend sum_recip by proper scaling factors:
-    sum_recip = sum_recip * 4T(π) / abs(det(lattice))
+    sum_recip *= 4T(π) / abs(det(lattice))
+    if forces !== nothing
+        forces_recip .*= 4T(π) / abs(det(lattice))
+    end
 
     #
     # Real-space sum
@@ -109,7 +127,12 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
 
         # Loop over R vectors for this shell patch
         for R in shell_indices(rsh)
-            for (ti, Zi) in zip(positions, charges), (tj, Zj) in zip(positions, charges)
+            for i = 1:length(positions), j = 1:length(positions)
+                ti = positions[i]
+                Zi = charges[i]
+                tj = positions[j]
+                Zj = charges[j]
+
                 # Avoid self-interaction
                 if rsh == 0 && ti == tj
                     continue
@@ -124,9 +147,17 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
 
                 any_term_contributes = true
                 sum_real += Zi * Zj * erfc(η * dist) / dist
-            end # iat, jat
+                if forces !== nothing
+                    # Use ForwardDiff here because I'm lazy. TODO do it properly
+                    forces_real[i] -= ForwardDiff.gradient(r -> (dist=norm(lattice * (r - tj - R)); Zi * Zj * erfc(η * dist) / dist), ti)
+                    forces_real[j] -= ForwardDiff.gradient(r -> (dist=norm(lattice * (ti - r - R)); Zi * Zj * erfc(η * dist) / dist), tj)
+                end
+            end # i,j
         end # R
         rsh += 1
     end
+    if forces !== nothing
+        forces .= (forces_recip .+ forces_real) ./ 2
+    end
     (sum_recip + sum_real) / 2  # Divide by 1/2 (because of double counting)
 end

src/forces.jl

@@ -0,0 +1,77 @@
+"""
+Computes minus the derivatives of the energy with respect to atomic positions.
+"""
+function forces(scfres)
+    # By generalized Hellmann-Feynman, dE/dR = ∂E/∂R. The atomic
+    # positions come up explicitly only in the external and nonlocal
+    # part of the energy,
+
+    # TODO this assumes that the build_external and build_nonlocal are consistent with model.atoms
+    # Find a way to move this computation closer to each term
+
+    # minus signs here because f = -∇E
+    ham = scfres.ham
+    basis = ham.basis
+    model = basis.model
+    Psi = scfres.Psi
+    ρ = scfres.ρ
+    T = real(eltype(ham))
+
+    forces = []
+    for (type, pos) in model.atoms
+        @assert type.psp != nothing
+        f = zeros(Vec3{T}, length(pos))
+
+        # Local part
+        # energy = sum of form_factor(G) * struct_factor(G) * rho(G)
+        # where struct_factor(G) = cis(-2T(π) * dot(G, r))
+        form_factors = [Complex{T}(1/sqrt(model.unit_cell_volume) *
+                                   eval_psp_local_fourier(type.psp, model.recip_lattice * G))
+                        for G in G_vectors(basis)]
+        form_factors[1] = 0
+
+        for (ir, r) in enumerate(pos)
+            f[ir] -= real(sum(conj(ρ.fourier[iG]) .* form_factors[iG] .* cis(-2T(π) * dot(G, r)) .* (-2T(π)) .* G .* im
+                              for (iG, G) in enumerate(G_vectors(basis))))
+        end
+
+        # Nonlocal part
+        C = build_projection_coefficients_(type.psp)
+        for (ir, r) in enumerate(pos)
+            fr = zeros(T, 3)
+            for idir = 1:3
+                for (ik, kpt) in enumerate(basis.kpoints)
+                    form_factors = build_form_factors(type.psp, [model.recip_lattice * (kpt.coordinate + G) for G in kpt.basis])
+                    structure_factors = [cis(-2T(π)*dot(kpt.coordinate + G, r)) for G in kpt.basis]
+                    P = structure_factors .* form_factors ./ sqrt(model.unit_cell_volume)
+                    dPdR = [-2T(π)*im*(kpt.coordinate + G)[idir] * cis(-2T(π)*dot(kpt.coordinate + G, r)) for G in kpt.basis] .*
+                           form_factors ./
+                           sqrt(model.unit_cell_volume)
+                    # TODO BLASify this further
+                    for iband = 1:size(Psi[ik], 2)
+                        psi = Psi[ik][:, iband]
+                        fr[idir] -= basis.kweights[ik] * scfres.occupation[ik][iband] * real(dot(psi, P*C*dPdR'*psi) + dot(psi, dPdR*C*P'*psi))
+                    end
+                end
+            end
+            f[ir] += fr
+        end
+
+        push!(forces, f)
+    end
+
+    # Add the ewald term
+    charges_all = [charge_ionic(type) for (type, positions) in model.atoms for pos in positions]
+    positions_all = [pos for (elem, positions) in model.atoms for pos in positions]
+    forces_ewald = zeros(Vec3{T}, length(positions_all))
+    energy_ewald(model.lattice, charges_all, positions_all; forces=forces_ewald)
+
+    offset = 1
+    for i = 1:length(model.atoms)
+        for j = 1:length(model.atoms[i][2])
+            forces[i][j] += forces_ewald[offset]
+            offset += 1
+        end
+    end
+    forces
+end

test/energy_ewald.jl

@@ -46,3 +46,21 @@ end
     γ_E = energy_ewald(lattice, charges, positions)
     @test abs(γ_E - ref) < 1e-7
 end
+
+@testset "Forces" begin
+    lattice = [0.0  5.131570667152971 5.131570667152971;
+               5.131570667152971 0.0 5.131570667152971;
+               5.131570667152971 5.131570667152971  0.0]
+    # perturb positions away from equilibrium to get nonzero force
+    positions = [ones(3)/8+rand(3)/20, -ones(3)/8]
+    charges = [14, 14]
+
+    forces = zeros(Vec3{Float64}, 2)
+    γ1 = energy_ewald(lattice, charges, positions, forces=forces)
+
+    # Compare forces to finite differences
+    disp = [rand(3)/20, rand(3)/20]
+    ε = 1e-8
+    γ2 = energy_ewald(lattice, charges, positions .+ ε .* disp)
+    @test abs((γ2-γ1)/ε + dot(disp, forces))/γ1 < 1e-6
+end

test/forces.jl

@@ -0,0 +1,52 @@
+using DFTK
+using Test
+function energy(pos)
+    # Calculation parameters
+    kgrid = [2, 1, 2]        # k-Point grid
+    supercell = [1, 1, 1]    # Lattice supercell
+    Ecut = 5                # kinetic energy cutoff in Hartree
+
+    # Setup silicon lattice
+    a = 10.263141334305942  # Silicon lattice constant in Bohr
+    lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
+    Si = AtomType(14, psp=load_psp("hgh/lda/Si-q4"))
+    # Si = AtomType(14)
+    atoms = [Si => pos]
+
+    model = Model(lattice;
+                  atoms=atoms,
+                  external=term_external(atoms),
+                  nonlocal=term_nonlocal(atoms),
+                  hartree=term_hartree(),
+                  xc=term_xc(:lda_xc_teter93))
+
+    kcoords, ksymops = bzmesh_ir_wedge(kgrid, lattice, atoms)
+    basis = PlaneWaveBasis(model, Ecut, kcoords, ksymops)
+
+    # Run SCF. Note Silicon is a semiconductor, so we use an insulator
+    # occupation scheme. This will cause warnings in some models, because
+    # e.g. in the :reduced_hf model silicon is a metal
+    n_bands_scf = Int(model.n_electrons / 2)
+    ham = Hamiltonian(basis, guess_density(basis))
+    scfres = self_consistent_field(ham, n_bands_scf, tol=1e-12)
+    ham = scfres.ham
+
+    energies = scfres.energies
+    sum(values(energies)), forces(scfres)
+end
+
+using Random
+Random.seed!(0)
+
+pos1 = [(ones(3)+.1*randn(3))/8, -ones(3)/8]
+disp = randn(3)
+ε = 1e-7
+pos2 = [pos1[1]+ε*disp, pos1[2]]
+
+E1, F1 = energy(pos1)
+E2, F2 = energy(pos2)
+
+diff_findiff = -(E2-E1)/ε
+diff_forces = dot(F1[1][1], disp)
+
+@test abs(diff_findiff - diff_forces) < 1e-6