diff --git a/src/structure.jl b/src/structure.jl index d5de01808..d9c584839 100644 --- a/src/structure.jl +++ b/src/structure.jl @@ -52,7 +52,7 @@ function estimate_integer_lattice_bounds(M::AbstractMatrix{T}, δ, shift=zeros(3 # then xi = = <= ||M^-T ei|| δ. inv_lattice_t = compute_inverse_lattice(M') xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i in 1:3] - + # Round up, unless exactly zero (in which case keep it zero in # order to just have one x vector for 1D or 2D systems) xlims = [xlim == 0 ? 0 : ceil(Int, xlim .- tol) for xlim in xlims] diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl index 39a21cbc4..b11b92ca2 100644 --- a/src/terms/pairwise.jl +++ b/src/terms/pairwise.jl @@ -1,3 +1,15 @@ +""" +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 + struct PairwisePotential V params @@ -10,8 +22,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 `100` by default, and gives the maximum (reduced) 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=100) params = Dict(minmax(key[1], key[2]) => value for (key, value) in params) @@ -56,11 +68,22 @@ 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. +# 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) - T = eltype(lattice) + 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 !isnothing(ph_disp) + @assert !isnothing(q) && !isnothing(forces) + T = promote_type(complex(T), eltype(ph_disp[1])) + @assert size(ph_disp) == size(positions) + end + if forces !== nothing @assert size(forces) == size(positions) forces_pairwise = copy(forces) @@ -89,22 +112,29 @@ function energy_pairwise(lattice, symbols, positions, V, params; 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) + ti = positions[i] + tj = positions[j] + R + if !isnothing(ph_disp) + ti += ph_disp[i] # * cis2pi(dot(q, zeros(3))) === 1 + # as we use the forces at the nuclei in the unit cell + tj += ph_disp[j] * cis2pi(dot(q, R)) + end + Δr = lattice * (ti .- tj) + dist = norm_cplx(Δ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_ddist = ForwardDiff.derivative(real(zero(eltype(dist)))) do ε + V(dist + ε, param_ij) + end 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 + forces .= forces_pairwise end energy end diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl index 34e7b5e2e..547eb4118 100644 --- a/src/workarounds/forwarddiff_rules.jl +++ b/src/workarounds/forwarddiff_rules.jl @@ -274,3 +274,26 @@ 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) && 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 diff --git a/test/phonons.jl b/test/phonons.jl new file mode 100644 index 000000000..31a24f867 --- /dev/null +++ b/test/phonons.jl @@ -0,0 +1,114 @@ +# TODO: Temporary, explanations too dry. To be changed with proper phonon computations. +using Test +using DFTK +using LinearAlgebra +using ForwardDiff +using StaticArrays + +# Convert back and forth between Vec3 and columnwise matrix +fold(x) = hcat(x...) +unfold(x) = Vec3.(eachcol(x)) + +function prepare_system(; n_scell=1) + positions = [[0.,0,0]] + 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) + 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) + + (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, max_radius=1e3) + blob = prepare_system(; 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 ε + 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, 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(; n_scell=1, max_radius=1e3, n_points=2) + blob = prepare_system(; n_scell) + 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, 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) + [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 + max_radius = 1e3 + tolerance = 1e-6 + n_points = 10 + + ph_bands = [] + 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 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, 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 diff --git a/test/runtests.jl b/test/runtests.jl index 3b14e2e8c..5b47483a5 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