From 8a5bd406d53ba1a89334b1240cfa44d91b3ba810 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Fri, 29 Apr 2022 14:57:20 +0200
Subject: [PATCH 01/28] Handling complex numbers for PairwisePotential

---
 src/terms/pairwise.jl | 53 ++++++++++++++++++++++++++++++-------------
 1 file changed, 37 insertions(+), 16 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 2e20699026..b132055972 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -1,3 +1,8 @@
+function norm_cplx(x)
+    # TODO: ForwardDiff bug (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324)
+    sqrt(sum(x.*x))
+end
+
 struct PairwisePotential
     V
     params
@@ -10,10 +15,10 @@ Lennard—Jones terms.
 The potential is dependent on the distance between to atomic positions and the pairwise
 atomic types:
 For a distance `d` between to atoms `A` and `B`, the potential is `V(d, params[(A, B)])`.
-The parameters `max_radius` is of `100` by default, and gives the maximum (reduced) distance
-between nuclei for which we consider interactions.
+The parameters `max_radius` is of `1000` by default, and gives the maximum (full,
+non-reduced) distance between nuclei for which we consider interactions.
 """
-function PairwisePotential(V, params; max_radius=100)
+function PairwisePotential(V, params; max_radius=1000)
     params = Dict(minmax(key[1], key[2]) => value for (key, value) in params)
     PairwisePotential(V, params, max_radius)
 end
@@ -36,8 +41,10 @@ end
 @timing "forces: Pairwise" function compute_forces(term::TermPairwisePotential,
                                                    basis::PlaneWaveBasis{T}, ψ, occ;
                                                    kwargs...) where {T}
-    forces = zero(basis.model.positions)
-    energy_pairwise(basis.model, term.V, term.params; term.max_radius, forces)
+    TT = eltype(lattice)
+    forces = zero(TT, basis.model.positions)
+    energy_pairwise(basis.model, term.V, term.params; max_radius=term.max_radius,
+                    forces=forces, kwargs...)
     forces
 end
 
@@ -56,11 +63,17 @@ end
 
 # This could be factorised with Ewald, but the use of `symbols` would slow down the
 # computationally intensive Ewald sums. So we leave it as it for now.
+# TODO: *Beware* of using ForwardDiff to derive this function with complex numbers, use
+# multiplications and not powers (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324)
 function energy_pairwise(lattice, symbols, positions, V, params;
-                         max_radius=100, forces=nothing)
-    T = eltype(lattice)
+                         max_radius=1000, forces=nothing, ph_disp=nothing, q=0)
     @assert length(symbols) == length(positions)
 
+    T = complex(eltype(lattice))
+    if ph_disp !== nothing
+        T = promote_type(T, eltype(ph_disp[1]))
+    end
+
     if forces !== nothing
         @assert size(forces) == size(positions)
         forces_pairwise = copy(forces)
@@ -96,24 +109,32 @@ function energy_pairwise(lattice, symbols, positions, V, params;
 
                 ti = positions[i]
                 tj = positions[j]
+                # Phonons `q` points
+                if !isnothing(ph_disp)
+                    ti += ph_disp[i] * cis(2T(π)*dot(q, zeros(3)))
+                    tj += ph_disp[j] * cis(2T(π)*dot(q, R)) + R
+                end
                 ai, aj = minmax(symbols[i], symbols[j])
                 param_ij =params[(ai, aj)]
 
-                Δr = lattice * (ti .- tj .- R)
-                dist = norm(Δr)
+                Δr = lattice * (ti .- tj)
+                dist = norm_cplx(Δr)
 
                 # the potential decays very quickly, so cut off at some point
-                dist > max_radius && continue
+                abs(dist) > max_radius && continue
 
                 any_term_contributes = true
-                energy_contribution = V(dist, param_ij)
+                energy_contribution = real(V(dist, param_ij))
                 sum_pairwise += energy_contribution
                 if forces !== nothing
-                    # We use ForwardDiff for the forces
-                    dE_ddist = ForwardDiff.derivative(d -> V(d, param_ij), dist)
-                    dE_dti = lattice' * dE_ddist / dist * Δr
+                    dE_ddist = ForwardDiff.derivative(real(zero(eltype(dist)))) do ε
+                        res = V(dist + ε, param_ij)
+                        [real(res), imag(res)]
+                    end |> x -> complex(x...)
+                    dE_dti = lattice' * ((dE_ddist / dist) * Δr)
+                    # We need to "break" the symmetry for phonons; at equilibrium, expect
+                    # the forces to be zero at machine precision.
                     forces_pairwise[i] -= dE_dti
-                    forces_pairwise[j] += dE_dti
                 end
             end # i,j
         end # R
@@ -121,7 +142,7 @@ function energy_pairwise(lattice, symbols, positions, V, params;
     end
     energy = sum_pairwise / 2  # Divide by 2 (because of double counting)
     if forces !== nothing
-        forces .= forces_pairwise ./ 2
+        forces .= forces_pairwise
     end
     energy
 end

From 281e161ec3d7df802f6ab2ab1cf3f113b1a4a8fc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Fri, 29 Apr 2022 15:34:14 +0200
Subject: [PATCH 02/28] comments from Antoine

---
 src/terms/pairwise.jl | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index b132055972..9b57d4f86e 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -1,3 +1,5 @@
+# We cannot use `LinearAlgebra.norm` with complex numbers due to the need to use its
+# analytic continuation
 function norm_cplx(x)
     # TODO: ForwardDiff bug (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324)
     sqrt(sum(x.*x))
@@ -15,8 +17,8 @@ Lennard—Jones terms.
 The potential is dependent on the distance between to atomic positions and the pairwise
 atomic types:
 For a distance `d` between to atoms `A` and `B`, the potential is `V(d, params[(A, B)])`.
-The parameters `max_radius` is of `1000` by default, and gives the maximum (full,
-non-reduced) distance between nuclei for which we consider interactions.
+The parameters `max_radius` is of `1000` by default, and gives the maximum (Cartesian)
+distance between nuclei for which we consider interactions.
 """
 function PairwisePotential(V, params; max_radius=1000)
     params = Dict(minmax(key[1], key[2]) => value for (key, value) in params)
@@ -64,14 +66,17 @@ end
 # This could be factorised with Ewald, but the use of `symbols` would slow down the
 # computationally intensive Ewald sums. So we leave it as it for now.
 # TODO: *Beware* of using ForwardDiff to derive this function with complex numbers, use
-# multiplications and not powers (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324)
+# multiplications and not powers (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324).
+# `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
+# compute the Fourier transform of the force constant matrix.
 function energy_pairwise(lattice, symbols, positions, V, params;
                          max_radius=1000, forces=nothing, ph_disp=nothing, q=0)
     @assert length(symbols) == length(positions)
 
-    T = complex(eltype(lattice))
+    T = eltype(lattice)
     if ph_disp !== nothing
-        T = promote_type(T, eltype(ph_disp[1]))
+        T = promote_type(complex(T), eltype(ph_disp[1]))
+        @assert size(ph_disp) == size(positions)
     end
 
     if forces !== nothing
@@ -111,7 +116,8 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                 tj = positions[j]
                 # Phonons `q` points
                 if !isnothing(ph_disp)
-                    ti += ph_disp[i] * cis(2T(π)*dot(q, zeros(3)))
+                    ti += ph_disp[i] # * cis(2T(π)*dot(q, zeros(3))) === 1
+                                     #  as we use the forces at the nuclei in the unit cell
                     tj += ph_disp[j] * cis(2T(π)*dot(q, R)) + R
                 end
                 ai, aj = minmax(symbols[i], symbols[j])

From 42b82c0e1684a696bd22bee55abfda89c49ebdde Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Fri, 29 Apr 2022 15:41:23 +0200
Subject: [PATCH 03/28] assert

---
 src/terms/pairwise.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 9b57d4f86e..f3008b1260 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -70,11 +70,12 @@ end
 # `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
 # compute the Fourier transform of the force constant matrix.
 function energy_pairwise(lattice, symbols, positions, V, params;
-                         max_radius=1000, forces=nothing, ph_disp=nothing, q=0)
+                         max_radius=1000, forces=nothing, ph_disp=nothing, q=nothing)
     @assert length(symbols) == length(positions)
 
     T = eltype(lattice)
     if ph_disp !== nothing
+        @assert q !== nothing
         T = promote_type(complex(T), eltype(ph_disp[1]))
         @assert size(ph_disp) == size(positions)
     end

From 4198c1c211ab5514e4c64838693970fe0ea0b805 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Mon, 2 May 2022 14:29:29 +0000
Subject: [PATCH 04/28] bug

---
 src/terms/pairwise.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index f3008b1260..8eedb380a0 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -114,12 +114,12 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                 rsh == 0 && i == j && continue
 
                 ti = positions[i]
-                tj = positions[j]
+                tj = positions[j] + R
                 # Phonons `q` points
                 if !isnothing(ph_disp)
                     ti += ph_disp[i] # * cis(2T(π)*dot(q, zeros(3))) === 1
                                      #  as we use the forces at the nuclei in the unit cell
-                    tj += ph_disp[j] * cis(2T(π)*dot(q, R)) + R
+                    tj += ph_disp[j] * cis(2T(π)*dot(q, R))
                 end
                 ai, aj = minmax(symbols[i], symbols[j])
                 param_ij =params[(ai, aj)]

From 1b9a3c7f1e49b4cf68eb408386bd3eb6bb3aaf1f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Tue, 3 May 2022 07:50:05 +0000
Subject: [PATCH 05/28] Update pairwise.jl

---
 src/terms/pairwise.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 8eedb380a0..e0f28f1d8d 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -43,7 +43,7 @@ end
 @timing "forces: Pairwise" function compute_forces(term::TermPairwisePotential,
                                                    basis::PlaneWaveBasis{T}, ψ, occ;
                                                    kwargs...) where {T}
-    TT = eltype(lattice)
+    TT = promote_type(T, eltype(basis.model.positions[1]))
     forces = zero(TT, basis.model.positions)
     energy_pairwise(basis.model, term.V, term.params; max_radius=term.max_radius,
                     forces=forces, kwargs...)

From 20c8a1764fd4c805ca0ddff6e6e3bc9159b65eeb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Tue, 3 May 2022 17:09:26 +0200
Subject: [PATCH 06/28] workarounds for complex exponentiation

---
 src/terms/pairwise.jl                | 23 +++++++++-------------
 src/workarounds/forwarddiff_rules.jl | 29 ++++++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 14 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index e0f28f1d8d..e5e1d1c768 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -1,8 +1,7 @@
 # We cannot use `LinearAlgebra.norm` with complex numbers due to the need to use its
 # analytic continuation
 function norm_cplx(x)
-    # TODO: ForwardDiff bug (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324)
-    sqrt(sum(x.*x))
+    sqrt(sum(x.^2))
 end
 
 struct PairwisePotential
@@ -17,10 +16,10 @@ Lennard—Jones terms.
 The potential is dependent on the distance between to atomic positions and the pairwise
 atomic types:
 For a distance `d` between to atoms `A` and `B`, the potential is `V(d, params[(A, B)])`.
-The parameters `max_radius` is of `1000` by default, and gives the maximum (Cartesian)
-distance between nuclei for which we consider interactions.
+The parameters `max_radius` is of `100` by default, and gives the maximum distance (in
+Cartesian coordinates) between nuclei for which we consider interactions.
 """
-function PairwisePotential(V, params; max_radius=1000)
+function PairwisePotential(V, params; max_radius=100)
     params = Dict(minmax(key[1], key[2]) => value for (key, value) in params)
     PairwisePotential(V, params, max_radius)
 end
@@ -43,8 +42,7 @@ end
 @timing "forces: Pairwise" function compute_forces(term::TermPairwisePotential,
                                                    basis::PlaneWaveBasis{T}, ψ, occ;
                                                    kwargs...) where {T}
-    TT = promote_type(T, eltype(basis.model.positions[1]))
-    forces = zero(TT, basis.model.positions)
+    forces = zero(basis.model.positions)
     energy_pairwise(basis.model, term.V, term.params; max_radius=term.max_radius,
                     forces=forces, kwargs...)
     forces
@@ -65,15 +63,13 @@ end
 
 # This could be factorised with Ewald, but the use of `symbols` would slow down the
 # computationally intensive Ewald sums. So we leave it as it for now.
-# TODO: *Beware* of using ForwardDiff to derive this function with complex numbers, use
-# multiplications and not powers (https://github.com/JuliaDiff/ForwardDiff.jl/issues/324).
 # `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
 # compute the Fourier transform of the force constant matrix.
 function energy_pairwise(lattice, symbols, positions, V, params;
-                         max_radius=1000, forces=nothing, ph_disp=nothing, q=nothing)
+                         max_radius=100, forces=nothing, ph_disp=nothing, q=nothing)
     @assert length(symbols) == length(positions)
 
-    T = eltype(lattice)
+    T = eltype(positions[1])
     if ph_disp !== nothing
         @assert q !== nothing
         T = promote_type(complex(T), eltype(ph_disp[1]))
@@ -135,9 +131,8 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                 sum_pairwise += energy_contribution
                 if forces !== nothing
                     dE_ddist = ForwardDiff.derivative(real(zero(eltype(dist)))) do ε
-                        res = V(dist + ε, param_ij)
-                        [real(res), imag(res)]
-                    end |> x -> complex(x...)
+                        V(dist + ε, param_ij)
+                    end
                     dE_dti = lattice' * ((dE_ddist / dist) * Δr)
                     # We need to "break" the symmetry for phonons; at equilibrium, expect
                     # the forces to be zero at machine precision.
diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index c2ab835c1a..8e89286fef 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -272,3 +272,32 @@ function Smearing.occupation(S::Smearing.FermiDirac, d::ForwardDiff.Dual{T}) whe
     end
     ForwardDiff.Dual{T}(Smearing.occupation(S, x), ∂occ * ForwardDiff.partials(d))
 end
+
+# Workarounds for issue https://github.com/JuliaDiff/ForwardDiff.jl/issues/324
+ForwardDiff.derivative(f, x::Complex) = throw(DimensionMismatch("derivative(f, x) expects that x is a real number (does not support Wirtinger derivatives). Separate real and imaginary parts of the input."))
+@inline ForwardDiff.extract_derivative(::Type{T}, y::Complex) where {T}       = zero(y)
+@inline function ForwardDiff.extract_derivative(::Type{T}, y::Complex{TD}) where {T, TD <: ForwardDiff.Dual}
+    complex(ForwardDiff.partials(T, real(y), 1), ForwardDiff.partials(T, imag(y), 1))
+end
+function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Complex{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
+    xx = complex(ForwardDiff.value(real(x)), ForwardDiff.value(imag(x)))
+    yy = complex(ForwardDiff.value(real(y)), ForwardDiff.value(imag(y)))
+    dx = complex.(ForwardDiff.partials(real(x)), ForwardDiff.partials(imag(x)))
+    dy = complex.(ForwardDiff.partials(real(y)), ForwardDiff.partials(imag(y)))
+
+    expv = xx^yy
+    ∂expv∂x = yy * xx^(yy-1)
+    ∂expv∂y = log(xx) * expv
+    dxexpv = ∂expv∂x * dx
+    # TODO: Fishy and should be checked, but seems to catch most cases
+    if iszero(xx) && ForwardDiff.isconstant(real(y)) && ForwardDiff.isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
+        dexpv = zero(expv)
+    elseif iszero(xx)
+        throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))
+    else
+        dyexpv = ∂expv∂y * dy
+        dexpv = dxexpv + dyexpv
+    end
+    complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
+            ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
+end

From ca7f803127633ef4b260db5ea1e84b6370b95ddc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Tue, 3 May 2022 15:28:46 +0000
Subject: [PATCH 07/28] Update pairwise.jl

---
 src/terms/pairwise.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index e5e1d1c768..60c3313246 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -113,7 +113,7 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                 tj = positions[j] + R
                 # Phonons `q` points
                 if !isnothing(ph_disp)
-                    ti += ph_disp[i] # * cis(2T(π)*dot(q, zeros(3))) === 1
+                    ti += ph_disp[i] # * cis(2T(π)*dot(q, zeros(3))) ≡ 1
                                      #  as we use the forces at the nuclei in the unit cell
                     tj += ph_disp[j] * cis(2T(π)*dot(q, R))
                 end

From 98b65a0a4f1372b8f0742098df3d3d9b909f1040 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Thu, 5 May 2022 09:36:44 +0000
Subject: [PATCH 08/28] Update pairwise.jl

---
 src/terms/pairwise.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 60c3313246..aa0cd41d98 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -127,7 +127,7 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                 abs(dist) > max_radius && continue
 
                 any_term_contributes = true
-                energy_contribution = real(V(dist, param_ij))
+                energy_contribution = V(dist, param_ij)
                 sum_pairwise += energy_contribution
                 if forces !== nothing
                     dE_ddist = ForwardDiff.derivative(real(zero(eltype(dist)))) do ε

From 4185299eeb9a41f6f049f911dce741a7082ca96c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Thu, 5 May 2022 14:24:00 +0000
Subject: [PATCH 09/28] Update pairwise.jl

---
 src/terms/pairwise.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index aa0cd41d98..afaf487e19 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -134,8 +134,7 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                         V(dist + ε, param_ij)
                     end
                     dE_dti = lattice' * ((dE_ddist / dist) * Δr)
-                    # We need to "break" the symmetry for phonons; at equilibrium, expect
-                    # the forces to be zero at machine precision.
+                    # For the phonons, we compute the forces only in the unit cell.
                     forces_pairwise[i] -= dE_dti
                 end
             end # i,j

From 45da86cd51e38f0991ec54eb47a601fa1ca24b69 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Thu, 5 May 2022 16:08:53 +0000
Subject: [PATCH 10/28] Update pairwise.jl

remove `kwargs` for `compute_forces`
---
 src/terms/pairwise.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index afaf487e19..fd0589cb1f 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -44,7 +44,7 @@ end
                                                    kwargs...) where {T}
     forces = zero(basis.model.positions)
     energy_pairwise(basis.model, term.V, term.params; max_radius=term.max_radius,
-                    forces=forces, kwargs...)
+                    forces=forces)
     forces
 end
 

From d4d48e934a24e3d294c0d82ccc422fa18e527b7e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Tue, 7 Jun 2022 17:57:16 +0200
Subject: [PATCH 11/28] testing ph_disp

---
 test/phonons.jl  | 137 +++++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl |   5 ++
 2 files changed, 142 insertions(+)
 create mode 100644 test/phonons.jl

diff --git a/test/phonons.jl b/test/phonons.jl
new file mode 100644
index 0000000000..7035455c13
--- /dev/null
+++ b/test/phonons.jl
@@ -0,0 +1,137 @@
+using Test
+using DFTK
+using LinearAlgebra
+using ForwardDiff
+using StaticArrays
+
+# ## Helper functions
+# Some functions that will be helpful for this example.
+fold(x)   = hcat(x...)
+unfold(x) = Vec3.(eachcol(x))
+
+const MAX_RADIUS = 1e3
+const TOLERANCE = 1e-6
+const N_POINTS = 10
+
+function prepare_system(; case=1)
+    positions = [[0.,0,0]]
+    for i in 1:case-1
+        push!(positions, i*ones(3)/case)
+    end
+
+    a = 5. * length(positions)
+    lattice = a * [[1 0 0.]; [0 0 0.]; [0 0 0.]]
+
+    s = DFTK.compute_inverse_lattice(lattice)
+    if case === 1
+        directions = [[s * [1,0,0]]]
+    elseif case === 2
+        directions = [[s * [1,0,0], s * [0,0,0]],
+                      [s * [0,0,0], s * [1,0,0]]]
+    elseif case === 3
+        directions = [[s * [1,0,0], s * [0,0,0], s * [0,0,0]],
+                      [s * [0,0,0], s * [1,0,0], s * [0,0,0]],
+                      [s * [0,0,0], s * [0,0,0], s * [1,0,0]]]
+    end
+
+    params = Dict((:X, :X) => (; ε=1, σ=a / length(positions) /2^(1/6)))
+    V(x, p) = 4*p.ε * ((p.σ/x)^12 - (p.σ/x)^6)
+
+    (positions=positions, lattice=lattice, directions=directions, params=params, V=V)
+end
+
+# Compute phonons for a one-dimensional pairwise potential for a set of `q = 0` using
+# supercell method
+function test_supercell_q0(; N_scell=1)
+    blob = prepare_system(; case=N_scell)
+    positions = blob.positions
+    lattice = blob.lattice
+    directions = blob.directions
+    params = blob.params
+    V = blob.V
+
+    s = DFTK.compute_inverse_lattice(lattice)
+    n_atoms = length(positions)
+
+    directions = [reshape(vcat([[i==j, 0.0, 0.0] for i in 1:n_atoms]...), 3, :) for j in 1:n_atoms]
+
+    Φ = Array{eltype(positions[1])}(undef, length(directions), n_atoms)
+    for (i, direction) in enumerate(directions)
+        Φ[i, :] = - ForwardDiff.derivative(0.0) do ε
+            n_positions = unfold(fold(positions) .+ ε .* s * direction)
+            forces = zeros(Vec3{complex(eltype(ε))}, length(positions))
+            DFTK.energy_pairwise(lattice, [:X for _ in positions],
+                                 n_positions, V, params; forces=forces, max_radius=MAX_RADIUS)
+            [(s * f)[1] for f in forces]
+        end
+    end
+    sqrt.(abs.(eigvals(Φ)))
+end
+
+# Compute phonons for a one-dimensional pairwise potential for a set of `q`-points
+function test_ph_disp(; case=1)
+    blob = prepare_system(; case=case)
+    positions = blob.positions
+    lattice = blob.lattice
+    directions = blob.directions
+    params = blob.params
+    V = blob.V
+
+    pairwise_ph = (q, d; forces=nothing) ->
+                     DFTK.energy_pairwise(lattice, [:X for _ in positions],
+                                          positions, V, params; q=[q, 0, 0],
+                                          ph_disp=d, forces=forces,
+                                          max_radius=MAX_RADIUS)
+
+    ph_bands = []
+    qs = -1/2:1/N_POINTS:1/2
+    for q in qs
+        as = ComplexF64[]
+        for d in directions
+            res = -ForwardDiff.derivative(0.0) do ε
+                forces = zeros(Vec3{complex(eltype(ε))}, length(positions))
+                pairwise_ph(q, ε*d; forces=forces)
+                [DFTK.compute_inverse_lattice(lattice)' *  f for f in forces]
+            end
+            [push!(as, r[1]) for r in res]
+        end
+        M = reshape(as, length(positions), :)
+        @assert ≈(norm(imag.(eigvals(M))), 0.0, atol=1e-8)
+        push!(ph_bands, sqrt.(abs.(real(eigvals(M)))))
+    end
+    return ph_bands
+end
+
+@testset "Phonon consistency" begin
+    ph_bands_1 = test_ph_disp(; case=1)
+    ph_bands_2 = test_ph_disp(; case=2)
+    ph_bands_3 = test_ph_disp(; case=3)
+
+    min_1 = minimum(hcat(ph_bands_1...))
+    max_1 = maximum(hcat(ph_bands_1...))
+    min_2 = minimum(hcat(ph_bands_2...))
+    max_2 = maximum(hcat(ph_bands_2...))
+    min_3 = minimum(hcat(ph_bands_3...))
+    max_3 = maximum(hcat(ph_bands_3...))
+
+    # Recover the same extremum for the system whatever case we test
+    @test ≈(min_1, min_2, atol=TOLERANCE)
+    @test ≈(min_1, min_3, atol=TOLERANCE)
+    @test ≈(max_1, max_2, atol=TOLERANCE)
+    @test ≈(max_1, max_3, atol=TOLERANCE)
+
+    r1_q0 = test_supercell_q0(; N_scell=1)
+    @assert length(r1_q0) == 1
+    ph_bands_1_q0 = ph_bands_1[N_POINTS÷2+1]
+    @test norm(r1_q0 - ph_bands_1_q0) < TOLERANCE
+
+    r2_q0 = sort(test_supercell_q0(; N_scell=2))
+    @assert length(r2_q0) == 2
+    ph_bands_2_q0 = ph_bands_2[N_POINTS÷2+1]
+    @test norm(r2_q0 - ph_bands_2_q0) < TOLERANCE
+
+    r3_q0 = sort(test_supercell_q0(; N_scell=3))
+    @assert length(r3_q0) == 3
+    ph_bands_3_q0 = ph_bands_3[N_POINTS÷2+1]
+    @test norm(r3_q0 - ph_bands_3_q0) < TOLERANCE
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3b14e2e8ca..5b47483a5e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -125,5 +125,10 @@ Random.seed!(0)
         include("forwarddiff.jl")
     end
 
+    # Phonons
+    if "all" in TAGS
+        include("phonons.jl")
+    end
+
     ("example" in TAGS) && include("runexamples.jl")
 end

From ef145191079f3725c3ad51fb520591ce0b7ba9b1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Tue, 7 Jun 2022 18:01:11 +0200
Subject: [PATCH 12/28] comment

---
 test/phonons.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/phonons.jl b/test/phonons.jl
index 7035455c13..a667f59a37 100644
--- a/test/phonons.jl
+++ b/test/phonons.jl
@@ -120,6 +120,7 @@ end
     @test ≈(max_1, max_2, atol=TOLERANCE)
     @test ≈(max_1, max_3, atol=TOLERANCE)
 
+    # Test consistency between supercell method at `q = 0` and direct `q`-points computations
     r1_q0 = test_supercell_q0(; N_scell=1)
     @assert length(r1_q0) == 1
     ph_bands_1_q0 = ph_bands_1[N_POINTS÷2+1]

From 9135a169d7dd75d7acd822f4dccd657c4beb2ddd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Tue, 7 Jun 2022 18:14:29 +0200
Subject: [PATCH 13/28] renaming

---
 test/phonons.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/phonons.jl b/test/phonons.jl
index a667f59a37..946b00f992 100644
--- a/test/phonons.jl
+++ b/test/phonons.jl
@@ -58,10 +58,10 @@ function test_supercell_q0(; N_scell=1)
     Φ = Array{eltype(positions[1])}(undef, length(directions), n_atoms)
     for (i, direction) in enumerate(directions)
         Φ[i, :] = - ForwardDiff.derivative(0.0) do ε
-            n_positions = unfold(fold(positions) .+ ε .* s * direction)
+            new_positions = unfold(fold(positions) .+ ε .* s * direction)
             forces = zeros(Vec3{complex(eltype(ε))}, length(positions))
-            DFTK.energy_pairwise(lattice, [:X for _ in positions],
-                                 n_positions, V, params; forces=forces, max_radius=MAX_RADIUS)
+            DFTK.energy_pairwise(lattice, [:X for _ in positions], new_positions, V, params;
+                                 forces=forces, max_radius=MAX_RADIUS)
             [(s * f)[1] for f in forces]
         end
     end

From 3f343eaa262f4ec2a818d5fb051aac893726d17d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Tue, 7 Jun 2022 18:24:42 +0200
Subject: [PATCH 14/28] factorisation

---
 test/phonons.jl | 42 ++++++++++++++----------------------------
 1 file changed, 14 insertions(+), 28 deletions(-)

diff --git a/test/phonons.jl b/test/phonons.jl
index 946b00f992..3f06a85c95 100644
--- a/test/phonons.jl
+++ b/test/phonons.jl
@@ -103,36 +103,22 @@ function test_ph_disp(; case=1)
 end
 
 @testset "Phonon consistency" begin
-    ph_bands_1 = test_ph_disp(; case=1)
-    ph_bands_2 = test_ph_disp(; case=2)
-    ph_bands_3 = test_ph_disp(; case=3)
-
-    min_1 = minimum(hcat(ph_bands_1...))
-    max_1 = maximum(hcat(ph_bands_1...))
-    min_2 = minimum(hcat(ph_bands_2...))
-    max_2 = maximum(hcat(ph_bands_2...))
-    min_3 = minimum(hcat(ph_bands_3...))
-    max_3 = maximum(hcat(ph_bands_3...))
+    ph_bands = []
+    for case in [1, 2, 3]
+        push!(ph_bands, test_ph_disp(; case=case))
+    end
 
     # Recover the same extremum for the system whatever case we test
-    @test ≈(min_1, min_2, atol=TOLERANCE)
-    @test ≈(min_1, min_3, atol=TOLERANCE)
-    @test ≈(max_1, max_2, atol=TOLERANCE)
-    @test ≈(max_1, max_3, atol=TOLERANCE)
+    for case in [2, 3]
+        @test ≈(minimum(fold(ph_bands[1])), minimum(fold(ph_bands[case])), atol=TOLERANCE)
+        @test ≈(maximum(fold(ph_bands[1])), maximum(fold(ph_bands[case])), atol=TOLERANCE)
+    end
 
     # Test consistency between supercell method at `q = 0` and direct `q`-points computations
-    r1_q0 = test_supercell_q0(; N_scell=1)
-    @assert length(r1_q0) == 1
-    ph_bands_1_q0 = ph_bands_1[N_POINTS÷2+1]
-    @test norm(r1_q0 - ph_bands_1_q0) < TOLERANCE
-
-    r2_q0 = sort(test_supercell_q0(; N_scell=2))
-    @assert length(r2_q0) == 2
-    ph_bands_2_q0 = ph_bands_2[N_POINTS÷2+1]
-    @test norm(r2_q0 - ph_bands_2_q0) < TOLERANCE
-
-    r3_q0 = sort(test_supercell_q0(; N_scell=3))
-    @assert length(r3_q0) == 3
-    ph_bands_3_q0 = ph_bands_3[N_POINTS÷2+1]
-    @test norm(r3_q0 - ph_bands_3_q0) < TOLERANCE
+    for N_scell in [1, 2, 3]
+        r_q0 = test_supercell_q0(; N_scell=N_scell)
+        @assert length(r_q0) == N_scell
+        ph_band_q0 = ph_bands[N_scell][N_POINTS÷2+1]
+        @test norm(r_q0 - ph_band_q0) < TOLERANCE
+    end
 end

From 0157e8045efc962ae2cd168e2bf0467bc0155601 Mon Sep 17 00:00:00 2001
From: Niklas Schmitz <niklas.f.schmitz@gmail.com>
Date: Wed, 8 Jun 2022 00:08:25 +0200
Subject: [PATCH 15/28] Move estimate_integer_bounds to structure.jl and
 support 1D and 2D systems

---
 src/structure.jl   | 12 ++++++++++++
 src/terms/ewald.jl |  8 --------
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/src/structure.jl b/src/structure.jl
index 8e9f8f21ac..16c56190e5 100644
--- a/src/structure.jl
+++ b/src/structure.jl
@@ -42,3 +42,15 @@ function diameter(lattice::AbstractMatrix)
     end
     diam
 end
+
+"""
+Estimate integer bounds for dense space loops from a given inequality ||Mx|| ≤ δ.
+For 1D and 2D systems the limit will be zero in the auxiliary dimensions.
+"""
+function estimate_integer_lattice_bounds(M, δ, shift=zeros(3))
+    # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ, 
+    # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
+    inv_lattice_t = compute_inverse_lattice(M')
+    xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i in 1:3]
+    ceil.(Int, xlims)
+end
diff --git a/src/terms/ewald.jl b/src/terms/ewald.jl
index e04a7a87cb..91025a140e 100644
--- a/src/terms/ewald.jl
+++ b/src/terms/ewald.jl
@@ -60,14 +60,6 @@ function energy_ewald(lattice, charges, positions; η=nothing, forces=nothing)
     energy_ewald(lattice, compute_recip_lattice(lattice), charges, positions; η, forces)
 end
 
-function estimate_integer_lattice_bounds(M, δ, shift=zeros(3))
-    # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ, 
-    # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
-    # Below code does not support non-3D systems.
-    xlims = [norm(inv(M')[:, i]) * δ + shift[i] for i in 1:3]
-    ceil.(Int, xlims)
-end
-
 # This could be factorised with Pairwise, but its use of `atom_types` would slow down this
 # computationally intensive Ewald sums. So we leave it as it for now.
 function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, forces=nothing)

From b4beeee24cdaa2aa26b6a1531d106f522e28558c Mon Sep 17 00:00:00 2001
From: Niklas Schmitz <niklas.f.schmitz@gmail.com>
Date: Wed, 8 Jun 2022 00:09:12 +0200
Subject: [PATCH 16/28] Rewrite pairwise without shelll_indices

---
 src/terms/pairwise.jl | 73 +++++++++++++++----------------------------
 1 file changed, 25 insertions(+), 48 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 2e20699026..6cf40878fd 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -66,59 +66,36 @@ function energy_pairwise(lattice, symbols, positions, V, params;
         forces_pairwise = copy(forces)
     end
 
-    # Function to return the indices corresponding
-    # to a particular shell.
-    # Not performance critical, so we do not type the function
-    max_shell(n, trivial) = trivial ? 0 : n
-    # Check if some coordinates are not used.
-    is_dim_trivial = [norm(lattice[:,i]) == 0 for i=1:3]
-    function shell_indices(nsh)
-        ish, jsh, ksh = max_shell.(nsh, is_dim_trivial)
-        [[i,j,k] for i in -ish:ish for j in -jsh:jsh for k in -ksh:ksh
-         if maximum(abs.([i,j,k])) == nsh]
-    end
+    # We use the bound  ||A (ti - tj - R)|| ≤ max_radius
+    # where A is the real-space lattice, rj and rk are atomic positions.
+    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
+    Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)    
 
     #
     # Energy loop
     #
     sum_pairwise::T = zero(T)
-    # Loop over real-space shells
-    rsh = 0 # Include R = 0
-    any_term_contributes = true
-    while any_term_contributes || rsh <= 1
-        any_term_contributes = false
-
-        # Loop over R vectors for this shell patch
-        for R in shell_indices(rsh)
-            for i = 1:length(positions), j = 1:length(positions)
-                # Avoid self-interaction
-                rsh == 0 && i == j && continue
-
-                ti = positions[i]
-                tj = positions[j]
-                ai, aj = minmax(symbols[i], symbols[j])
-                param_ij =params[(ai, aj)]
-
-                Δr = lattice * (ti .- tj .- R)
-                dist = norm(Δr)
-
-                # the potential decays very quickly, so cut off at some point
-                dist > max_radius && continue
-
-                any_term_contributes = true
-                energy_contribution = V(dist, param_ij)
-                sum_pairwise += energy_contribution
-                if forces !== nothing
-                    # We use ForwardDiff for the forces
-                    dE_ddist = ForwardDiff.derivative(d -> V(d, param_ij), dist)
-                    dE_dti = lattice' * dE_ddist / dist * Δr
-                    forces_pairwise[i] -= dE_dti
-                    forces_pairwise[j] += dE_dti
-                end
-            end # i,j
-        end # R
-        rsh += 1
-    end
+    # Loop over real-space
+    for R1 in -Rlims[1]:Rlims[1], R2 in -Rlims[2]:Rlims[2], R3 in -Rlims[3]:Rlims[3]
+        R = Vec3(R1, R2, R3)
+        for i = 1:length(positions), j = 1:length(positions)
+            # Avoid self-interaction
+            R == zero(R) && i == j && continue
+            ai, aj = minmax(symbols[i], symbols[j])
+            param_ij = params[(ai, aj)]
+            Δr = lattice * (positions[i] - positions[j] - R)
+            dist = norm(Δr)
+            energy_contribution = V(dist, param_ij)
+            sum_pairwise += energy_contribution
+            if forces !== nothing
+                # We use ForwardDiff for the forces
+                dE_ddist = ForwardDiff.derivative(d -> V(d, param_ij), dist)
+                dE_dti = lattice' * dE_ddist / dist * Δr
+                forces_pairwise[i] -= dE_dti
+                forces_pairwise[j] += dE_dti
+            end
+        end # i,j
+    end # R
     energy = sum_pairwise / 2  # Divide by 2 (because of double counting)
     if forces !== nothing
         forces .= forces_pairwise ./ 2

From cbe29804c5efd43fb218691120d0f205d948dffe Mon Sep 17 00:00:00 2001
From: Niklas Schmitz <niklas.f.schmitz@gmail.com>
Date: Wed, 8 Jun 2022 00:13:09 +0200
Subject: [PATCH 17/28] trim whitespace

---
 src/structure.jl      | 2 +-
 src/terms/ewald.jl    | 8 ++++----
 src/terms/pairwise.jl | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/structure.jl b/src/structure.jl
index 16c56190e5..0fe31e3c7a 100644
--- a/src/structure.jl
+++ b/src/structure.jl
@@ -48,7 +48,7 @@ Estimate integer bounds for dense space loops from a given inequality ||Mx|| ≤
 For 1D and 2D systems the limit will be zero in the auxiliary dimensions.
 """
 function estimate_integer_lattice_bounds(M, δ, shift=zeros(3))
-    # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ, 
+    # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ,
     # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
     inv_lattice_t = compute_inverse_lattice(M')
     xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i in 1:3]
diff --git a/src/terms/ewald.jl b/src/terms/ewald.jl
index 91025a140e..ea0022a8d7 100644
--- a/src/terms/ewald.jl
+++ b/src/terms/ewald.jl
@@ -83,17 +83,17 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, fo
     max_erfc_arg = sqrt(max_exp_arg)  # erfc(x) ~= exp(-x^2)/(sqrt(π)x) for large x
 
     # Precomputing summation bounds from cutoffs.
-    # In the reciprocal-space term we have exp(-||B G||^2 / 4η^2), 
-    # where B is the reciprocal-space lattice, and 
+    # In the reciprocal-space term we have exp(-||B G||^2 / 4η^2),
+    # where B is the reciprocal-space lattice, and
     # thus use the bound  ||B G|| / 2η ≤ sqrt(max_exp_arg)
     Glims = estimate_integer_lattice_bounds(recip_lattice, sqrt(max_exp_arg) * 2η)
 
-    # In the real-space term we have erfc(η ||A(rj - rk - R)||), 
+    # In the real-space term we have erfc(η ||A(rj - rk - R)||),
     # where A is the real-space lattice, rj and rk are atomic positions and
     # thus use the bound  ||A(rj - rk - R)|| * η ≤ max_erfc_arg
     poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
     Rlims = estimate_integer_lattice_bounds(lattice, max_erfc_arg / η, poslims)
-        
+
     #
     # Reciprocal space sum
     #
diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 6cf40878fd..9ca71675b2 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -69,7 +69,7 @@ function energy_pairwise(lattice, symbols, positions, V, params;
     # We use the bound  ||A (ti - tj - R)|| ≤ max_radius
     # where A is the real-space lattice, rj and rk are atomic positions.
     poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
-    Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)    
+    Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)
 
     #
     # Energy loop

From de7405c7b035a352c67de468aa1ac3cee7948d78 Mon Sep 17 00:00:00 2001
From: Niklas Schmitz <niklas.f.schmitz@gmail.com>
Date: Wed, 8 Jun 2022 00:28:25 +0200
Subject: [PATCH 18/28] Fix pairwise bound comments

---
 src/terms/pairwise.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 9ca71675b2..d0aaf76547 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -66,7 +66,8 @@ function energy_pairwise(lattice, symbols, positions, V, params;
         forces_pairwise = copy(forces)
     end
 
-    # We use the bound  ||A (ti - tj - R)|| ≤ max_radius
+    # The potential V(dist) decays very quickly with dist = ||A (ri - rj - R)||,
+    # so we cut off at some point. We use the bound  ||A (ri - rj - R)|| ≤ max_radius
     # where A is the real-space lattice, rj and rk are atomic positions.
     poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
     Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)

From 86141e03d263ce6c6e7ff07c7d9afa8720257ae4 Mon Sep 17 00:00:00 2001
From: Niklas Schmitz <niklas.f.schmitz@gmail.com>
Date: Wed, 8 Jun 2022 00:29:16 +0200
Subject: [PATCH 19/28] Fix comment

---
 src/terms/pairwise.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index d0aaf76547..e94cd0d236 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -66,8 +66,8 @@ function energy_pairwise(lattice, symbols, positions, V, params;
         forces_pairwise = copy(forces)
     end
 
-    # The potential V(dist) decays very quickly with dist = ||A (ri - rj - R)||,
-    # so we cut off at some point. We use the bound  ||A (ri - rj - R)|| ≤ max_radius
+    # The potential V(dist) decays very quickly with dist = ||A (rj - rk - R)||,
+    # so we cut off at some point. We use the bound  ||A (rj - rk - R)|| ≤ max_radius
     # where A is the real-space lattice, rj and rk are atomic positions.
     poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
     Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)

From aa2bf98307612ee3530a378d3ad019b106d97f72 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Wed, 8 Jun 2022 11:39:40 +0200
Subject: [PATCH 20/28] first batch of modifications

---
 src/terms/pairwise.jl |  6 ++--
 test/phonons.jl       | 65 ++++++++++++++++++-------------------------
 2 files changed, 29 insertions(+), 42 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index fd0589cb1f..8eb1a7dde7 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -43,8 +43,7 @@ end
                                                    basis::PlaneWaveBasis{T}, ψ, occ;
                                                    kwargs...) where {T}
     forces = zero(basis.model.positions)
-    energy_pairwise(basis.model, term.V, term.params; max_radius=term.max_radius,
-                    forces=forces)
+    energy_pairwise(basis.model, term.V, term.params; term.max_radius, forces)
     forces
 end
 
@@ -67,6 +66,7 @@ end
 # compute the Fourier transform of the force constant matrix.
 function energy_pairwise(lattice, symbols, positions, V, params;
                          max_radius=100, forces=nothing, ph_disp=nothing, q=nothing)
+    isnothing(ph_disp) && @assert isnothing(q)
     @assert length(symbols) == length(positions)
 
     T = eltype(positions[1])
@@ -111,7 +111,6 @@ function energy_pairwise(lattice, symbols, positions, V, params;
 
                 ti = positions[i]
                 tj = positions[j] + R
-                # Phonons `q` points
                 if !isnothing(ph_disp)
                     ti += ph_disp[i] # * cis(2T(π)*dot(q, zeros(3))) ≡ 1
                                      #  as we use the forces at the nuclei in the unit cell
@@ -134,7 +133,6 @@ function energy_pairwise(lattice, symbols, positions, V, params;
                         V(dist + ε, param_ij)
                     end
                     dE_dti = lattice' * ((dE_ddist / dist) * Δr)
-                    # For the phonons, we compute the forces only in the unit cell.
                     forces_pairwise[i] -= dE_dti
                 end
             end # i,j
diff --git a/test/phonons.jl b/test/phonons.jl
index 3f06a85c95..c530fcfb39 100644
--- a/test/phonons.jl
+++ b/test/phonons.jl
@@ -4,35 +4,21 @@ using LinearAlgebra
 using ForwardDiff
 using StaticArrays
 
-# ## Helper functions
-# Some functions that will be helpful for this example.
+# Convert back and forth between Vec3 and columnwise matrix
 fold(x)   = hcat(x...)
 unfold(x) = Vec3.(eachcol(x))
 
-const MAX_RADIUS = 1e3
-const TOLERANCE = 1e-6
-const N_POINTS = 10
-
-function prepare_system(; case=1)
+function prepare_system(; n_scell=1)
     positions = [[0.,0,0]]
-    for i in 1:case-1
-        push!(positions, i*ones(3)/case)
+    for i in 1:n_scell-1
+        push!(positions, i*ones(3)/n_scell)
     end
 
     a = 5. * length(positions)
     lattice = a * [[1 0 0.]; [0 0 0.]; [0 0 0.]]
 
     s = DFTK.compute_inverse_lattice(lattice)
-    if case === 1
-        directions = [[s * [1,0,0]]]
-    elseif case === 2
-        directions = [[s * [1,0,0], s * [0,0,0]],
-                      [s * [0,0,0], s * [1,0,0]]]
-    elseif case === 3
-        directions = [[s * [1,0,0], s * [0,0,0], s * [0,0,0]],
-                      [s * [0,0,0], s * [1,0,0], s * [0,0,0]],
-                      [s * [0,0,0], s * [0,0,0], s * [1,0,0]]]
-    end
+    directions = [[s * [i==j,0,0] for i in 1:n_scell] for j in 1:n_scell]
 
     params = Dict((:X, :X) => (; ε=1, σ=a / length(positions) /2^(1/6)))
     V(x, p) = 4*p.ε * ((p.σ/x)^12 - (p.σ/x)^6)
@@ -42,8 +28,8 @@ end
 
 # Compute phonons for a one-dimensional pairwise potential for a set of `q = 0` using
 # supercell method
-function test_supercell_q0(; N_scell=1)
-    blob = prepare_system(; case=N_scell)
+function test_supercell_q0(; n_scell=1, max_radius=1e3)
+    blob = prepare_system(; n_scell)
     positions = blob.positions
     lattice = blob.lattice
     directions = blob.directions
@@ -61,7 +47,7 @@ function test_supercell_q0(; N_scell=1)
             new_positions = unfold(fold(positions) .+ ε .* s * direction)
             forces = zeros(Vec3{complex(eltype(ε))}, length(positions))
             DFTK.energy_pairwise(lattice, [:X for _ in positions], new_positions, V, params;
-                                 forces=forces, max_radius=MAX_RADIUS)
+                                 forces, max_radius)
             [(s * f)[1] for f in forces]
         end
     end
@@ -69,8 +55,8 @@ function test_supercell_q0(; N_scell=1)
 end
 
 # Compute phonons for a one-dimensional pairwise potential for a set of `q`-points
-function test_ph_disp(; case=1)
-    blob = prepare_system(; case=case)
+function test_ph_disp(; n_scell=1, max_radius=1e3, n_points=2)
+    blob = prepare_system(; n_scell)
     positions = blob.positions
     lattice = blob.lattice
     directions = blob.directions
@@ -80,17 +66,16 @@ function test_ph_disp(; case=1)
     pairwise_ph = (q, d; forces=nothing) ->
                      DFTK.energy_pairwise(lattice, [:X for _ in positions],
                                           positions, V, params; q=[q, 0, 0],
-                                          ph_disp=d, forces=forces,
-                                          max_radius=MAX_RADIUS)
+                                          ph_disp=d, forces, max_radius)
 
     ph_bands = []
-    qs = -1/2:1/N_POINTS:1/2
+    qs = -1/2:1/n_points:1/2
     for q in qs
         as = ComplexF64[]
         for d in directions
             res = -ForwardDiff.derivative(0.0) do ε
                 forces = zeros(Vec3{complex(eltype(ε))}, length(positions))
-                pairwise_ph(q, ε*d; forces=forces)
+                pairwise_ph(q, ε*d; forces)
                 [DFTK.compute_inverse_lattice(lattice)' *  f for f in forces]
             end
             [push!(as, r[1]) for r in res]
@@ -103,22 +88,26 @@ function test_ph_disp(; case=1)
 end
 
 @testset "Phonon consistency" begin
+    max_radius = 1e3
+    tolerance = 1e-6
+    n_points = 10
+
     ph_bands = []
-    for case in [1, 2, 3]
-        push!(ph_bands, test_ph_disp(; case=case))
+    for n_scell in [1, 2, 3]
+        push!(ph_bands, test_ph_disp(; n_scell, max_radius, n_points))
     end
 
     # Recover the same extremum for the system whatever case we test
-    for case in [2, 3]
-        @test ≈(minimum(fold(ph_bands[1])), minimum(fold(ph_bands[case])), atol=TOLERANCE)
-        @test ≈(maximum(fold(ph_bands[1])), maximum(fold(ph_bands[case])), atol=TOLERANCE)
+    for n_scell in [2, 3]
+        @test ≈(minimum(fold(ph_bands[1])), minimum(fold(ph_bands[n_scell])), atol=tolerance)
+        @test ≈(maximum(fold(ph_bands[1])), maximum(fold(ph_bands[n_scell])), atol=tolerance)
     end
 
     # Test consistency between supercell method at `q = 0` and direct `q`-points computations
-    for N_scell in [1, 2, 3]
-        r_q0 = test_supercell_q0(; N_scell=N_scell)
-        @assert length(r_q0) == N_scell
-        ph_band_q0 = ph_bands[N_scell][N_POINTS÷2+1]
-        @test norm(r_q0 - ph_band_q0) < TOLERANCE
+    for n_scell in [1, 2, 3]
+        r_q0 = test_supercell_q0(; n_scell, max_radius)
+        @assert length(r_q0) == n_scell
+        ph_band_q0 = ph_bands[n_scell][n_points÷2+1]
+        @test norm(r_q0 - ph_band_q0) < tolerance
     end
 end

From afc7c143e07924ff76610ada1ec9cdbcbc4161db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Wed, 8 Jun 2022 13:47:08 +0200
Subject: [PATCH 21/28] some more

---
 src/terms/pairwise.jl                | 27 ++++++++++++++++++--------
 src/workarounds/forwarddiff_rules.jl | 29 ----------------------------
 2 files changed, 19 insertions(+), 37 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index 8eb1a7dde7..6bb16bc997 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -1,5 +1,11 @@
-# We cannot use `LinearAlgebra.norm` with complex numbers due to the need to use its
-# analytic continuation
+"""
+Complex-analytic extension of `LinearAlgebra.norm(x)` from real to complex inputs.
+Needed for phonons as we want to perform a matrix-vector product `f'(x)·h`, where `f` is
+a real-to-real function and `h` a complex vector. To do this using automatic
+differentiation, we can extend analytically f to accept complex inputs, then differentiate
+`t -> f(x+t·h)`. This will fail if non-analytic functions like norm are used for complex
+inputs, and therefore we have to redefine it.
+"""
 function norm_cplx(x)
     sqrt(sum(x.^2))
 end
@@ -60,18 +66,23 @@ function energy_pairwise(model::Model{T}, V, params; kwargs...) where {T}
 end
 
 
-# This could be factorised with Ewald, but the use of `symbols` would slow down the
-# computationally intensive Ewald sums. So we leave it as it for now.
-# `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
-# compute the Fourier transform of the force constant matrix.
+"""
+This could be factorised with Ewald, but the use of `symbols` would slow down the
+computationally intensive Ewald sums. So we leave it as it for now.
+`q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
+compute the Fourier transform of the force constant matrix. Only the computations of the
+forces make sense.
+For phonons computations, this gives the forces of particles `ti` in the unit cell w.r.t. to
+a displacement of the particles `tj` of the form `ph_disp·e^{i q·R}`.
+"""
 function energy_pairwise(lattice, symbols, positions, V, params;
                          max_radius=100, forces=nothing, ph_disp=nothing, q=nothing)
     isnothing(ph_disp) && @assert isnothing(q)
     @assert length(symbols) == length(positions)
 
     T = eltype(positions[1])
-    if ph_disp !== nothing
-        @assert q !== nothing
+    if !isnothing(ph_disp)
+        @assert !isnothing(q) && !isnothing(forces)
         T = promote_type(complex(T), eltype(ph_disp[1]))
         @assert size(ph_disp) == size(positions)
     end
diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index 8e89286fef..c2ab835c1a 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -272,32 +272,3 @@ function Smearing.occupation(S::Smearing.FermiDirac, d::ForwardDiff.Dual{T}) whe
     end
     ForwardDiff.Dual{T}(Smearing.occupation(S, x), ∂occ * ForwardDiff.partials(d))
 end
-
-# Workarounds for issue https://github.com/JuliaDiff/ForwardDiff.jl/issues/324
-ForwardDiff.derivative(f, x::Complex) = throw(DimensionMismatch("derivative(f, x) expects that x is a real number (does not support Wirtinger derivatives). Separate real and imaginary parts of the input."))
-@inline ForwardDiff.extract_derivative(::Type{T}, y::Complex) where {T}       = zero(y)
-@inline function ForwardDiff.extract_derivative(::Type{T}, y::Complex{TD}) where {T, TD <: ForwardDiff.Dual}
-    complex(ForwardDiff.partials(T, real(y), 1), ForwardDiff.partials(T, imag(y), 1))
-end
-function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Complex{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
-    xx = complex(ForwardDiff.value(real(x)), ForwardDiff.value(imag(x)))
-    yy = complex(ForwardDiff.value(real(y)), ForwardDiff.value(imag(y)))
-    dx = complex.(ForwardDiff.partials(real(x)), ForwardDiff.partials(imag(x)))
-    dy = complex.(ForwardDiff.partials(real(y)), ForwardDiff.partials(imag(y)))
-
-    expv = xx^yy
-    ∂expv∂x = yy * xx^(yy-1)
-    ∂expv∂y = log(xx) * expv
-    dxexpv = ∂expv∂x * dx
-    # TODO: Fishy and should be checked, but seems to catch most cases
-    if iszero(xx) && ForwardDiff.isconstant(real(y)) && ForwardDiff.isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
-        dexpv = zero(expv)
-    elseif iszero(xx)
-        throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))
-    else
-        dyexpv = ∂expv∂y * dy
-        dexpv = dxexpv + dyexpv
-    end
-    complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
-            ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
-end

From f88dd86fdcd630cc17dc40c792e9df34221a859d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Wed, 8 Jun 2022 15:28:38 +0200
Subject: [PATCH 22/28] workaround back with a vengeance

---
 src/workarounds/forwarddiff_rules.jl | 41 ++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index c2ab835c1a..b50766232f 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -272,3 +272,44 @@ function Smearing.occupation(S::Smearing.FermiDirac, d::ForwardDiff.Dual{T}) whe
     end
     ForwardDiff.Dual{T}(Smearing.occupation(S, x), ∂occ * ForwardDiff.partials(d))
 end
+
+# Fix for https://github.com/JuliaDiff/ForwardDiff.jl/issues/514
+function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Complex{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
+    xx = complex(ForwardDiff.value(real(x)), ForwardDiff.value(imag(x)))
+    yy = complex(ForwardDiff.value(real(y)), ForwardDiff.value(imag(y)))
+    dx = complex.(ForwardDiff.partials(real(x)), ForwardDiff.partials(imag(x)))
+    dy = complex.(ForwardDiff.partials(real(y)), ForwardDiff.partials(imag(y)))
+
+    expv = xx^yy
+    ∂expv∂x = yy * xx^(yy-1)
+    ∂expv∂y = log(xx) * expv
+    dxexpv = ∂expv∂x * dx
+    if iszero(xx) && isconstant(real(y)) && isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
+        dexpv = zero(expv)
+    elseif iszero(xx)
+        throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))
+    else
+        dyexpv = ∂expv∂y * dy
+        dexpv = dxexpv + dyexpv
+    end
+    complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
+            ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
+end
+function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Int64) where {T,V,N}
+    xx = complex(ForwardDiff.value(real(x)), ForwardDiff.value(imag(x)))
+    dx = complex.(ForwardDiff.partials(real(x)), ForwardDiff.partials(imag(x)))
+
+    expv = xx^y
+    ∂expv∂x = y * xx^(y-1)
+    dxexpv = ∂expv∂x * dx
+    if iszero(xx) && isconstant(real(y)) && isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
+        dexpv = zero(expv)
+    elseif iszero(xx)
+        throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))
+    else
+        dexpv = dxexpv
+    end
+    complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
+            ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
+end
+

From 1d1d7906e73335677c9ffd6563d7f3ae36679a52 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Wed, 8 Jun 2022 15:49:14 +0200
Subject: [PATCH 23/28] bugfix

---
 src/workarounds/forwarddiff_rules.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index b50766232f..be68fc4abc 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -302,7 +302,7 @@ function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Int64) where {T,V,N}
     expv = xx^y
     ∂expv∂x = y * xx^(y-1)
     dxexpv = ∂expv∂x * dx
-    if iszero(xx) && isconstant(real(y)) && isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
+    if iszero(xx) && imag(y) === zero(imag(y)) && real(y) > 0
         dexpv = zero(expv)
     elseif iszero(xx)
         throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))

From 610d40495b1b58c0fd0b4919f1a41a221347d521 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Thu, 9 Jun 2022 14:55:41 +0000
Subject: [PATCH 24/28] Update forwarddiff_rules.jl

---
 src/workarounds/forwarddiff_rules.jl | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index be68fc4abc..5b024450f7 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -295,21 +295,3 @@ function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Complex{ForwardDiff.Dua
     complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
             ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
 end
-function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Int64) where {T,V,N}
-    xx = complex(ForwardDiff.value(real(x)), ForwardDiff.value(imag(x)))
-    dx = complex.(ForwardDiff.partials(real(x)), ForwardDiff.partials(imag(x)))
-
-    expv = xx^y
-    ∂expv∂x = y * xx^(y-1)
-    dxexpv = ∂expv∂x * dx
-    if iszero(xx) && imag(y) === zero(imag(y)) && real(y) > 0
-        dexpv = zero(expv)
-    elseif iszero(xx)
-        throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))
-    else
-        dexpv = dxexpv
-    end
-    complex(ForwardDiff.Dual{T,V,N}(real(expv), ForwardDiff.Partials{N,V}(tuple(real(dexpv)...))),
-            ForwardDiff.Dual{T,V,N}(imag(expv), ForwardDiff.Partials{N,V}(tuple(imag(dexpv)...))))
-end
-

From 284a9157a587a61a0cf1420159abda1912855413 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Mon, 13 Jun 2022 10:20:41 +0200
Subject: [PATCH 25/28] Antoine's comment

---
 src/terms/pairwise.jl | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index c709ddf242..bbe208ec79 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -66,15 +66,13 @@ function energy_pairwise(model::Model{T}, V, params; kwargs...) where {T}
 end
 
 
-"""
-This could be factorised with Ewald, but the use of `symbols` would slow down the
-computationally intensive Ewald sums. So we leave it as it for now.
-`q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
-compute the Fourier transform of the force constant matrix. Only the computations of the
-forces make sense.
-For phonons computations, this gives the forces of particles `ti` in the unit cell w.r.t. to
-a displacement of the particles `tj` of the form `ph_disp·e^{i q·R}`.
-"""
+# This could be factorised with Ewald, but the use of `symbols` would slow down the
+# computationally intensive Ewald sums. So we leave it as it for now.
+# `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
+# compute the Fourier transform of the force constant matrix. Only the computations of the
+# forces make sense.
+# Computes the local energy and forces on the atoms of the reference unit cell 0, for an
+# infinite array of atoms at positions r_{iR} = positions[i] + R + ph_disp[i]*e^{iq·R}.
 function energy_pairwise(lattice, symbols, positions, V, params;
                          max_radius=100, forces=nothing, ph_disp=nothing, q=nothing)
     isnothing(ph_disp) && @assert isnothing(q)

From 74908cfaef78a43cee2ea2c46e5aae85106ca5e3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?= <polack@hollved.net>
Date: Mon, 13 Jun 2022 16:47:35 +0200
Subject: [PATCH 26/28] bugfix

---
 src/workarounds/forwarddiff_rules.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl
index 5b024450f7..d41b195b70 100644
--- a/src/workarounds/forwarddiff_rules.jl
+++ b/src/workarounds/forwarddiff_rules.jl
@@ -284,7 +284,7 @@ function Base.:^(x::Complex{ForwardDiff.Dual{T,V,N}}, y::Complex{ForwardDiff.Dua
     ∂expv∂x = yy * xx^(yy-1)
     ∂expv∂y = log(xx) * expv
     dxexpv = ∂expv∂x * dx
-    if iszero(xx) && isconstant(real(y)) && isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
+    if iszero(xx) && ForwardDiff.isconstant(real(y)) && ForwardDiff.isconstant(imag(y)) && imag(y) === zero(imag(y)) && real(y) > 0
         dexpv = zero(expv)
     elseif iszero(xx)
         throw(DomainError(x, "mantissa cannot be zero for complex exponentiation"))

From 0927546569bc821edf415b109abded0abb40a729 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Tue, 28 Jun 2022 14:49:32 +0000
Subject: [PATCH 27/28] Update phonons.jl

---
 test/phonons.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/phonons.jl b/test/phonons.jl
index c530fcfb39..31a24f867d 100644
--- a/test/phonons.jl
+++ b/test/phonons.jl
@@ -1,3 +1,4 @@
+# TODO: Temporary, explanations too dry. To be changed with proper phonon computations.
 using Test
 using DFTK
 using LinearAlgebra

From 1c3d395e7c021beb3442022623d0135681712cdf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89tienne=20Polack?=
 <87805034+epolack@users.noreply.github.com>
Date: Tue, 28 Jun 2022 14:49:58 +0000
Subject: [PATCH 28/28] Update pairwise.jl

---
 src/terms/pairwise.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index bbe208ec79..b11b92ca2d 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -69,8 +69,7 @@ end
 # This could be factorised with Ewald, but the use of `symbols` would slow down the
 # computationally intensive Ewald sums. So we leave it as it for now.
 # `q` is the phonon `q`-point (`Vec3`), and `ph_disp` a list of `Vec3` displacements to
-# compute the Fourier transform of the force constant matrix. Only the computations of the
-# forces make sense.
+# compute the Fourier transform of the force constant matrix.
 # Computes the local energy and forces on the atoms of the reference unit cell 0, for an
 # infinite array of atoms at positions r_{iR} = positions[i] + R + ph_disp[i]*e^{iq·R}.
 function energy_pairwise(lattice, symbols, positions, V, params;