Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Handling complex numbers for PairwisePotential #655

Merged
merged 31 commits into from
Jun 28, 2022
Merged
Changes from 1 commit
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
8a5bd40
Handling complex numbers for PairwisePotential
epolack Apr 29, 2022
281e161
comments from Antoine
epolack Apr 29, 2022
42b82c0
assert
epolack Apr 29, 2022
4198c1c
bug
epolack May 2, 2022
1b9a3c7
Update pairwise.jl
epolack May 3, 2022
20c8a17
workarounds for complex exponentiation
epolack May 3, 2022
ca7f803
Update pairwise.jl
epolack May 3, 2022
98b65a0
Update pairwise.jl
epolack May 5, 2022
4185299
Update pairwise.jl
epolack May 5, 2022
45da86c
Update pairwise.jl
epolack May 5, 2022
d4d48e9
testing ph_disp
epolack Jun 7, 2022
ef14519
comment
epolack Jun 7, 2022
9135a16
renaming
epolack Jun 7, 2022
3f343ea
factorisation
epolack Jun 7, 2022
0157e80
Move estimate_integer_bounds to structure.jl and support 1D and 2D sy…
niklasschmitz Jun 7, 2022
b4beeee
Rewrite pairwise without shelll_indices
niklasschmitz Jun 7, 2022
cbe2980
trim whitespace
niklasschmitz Jun 7, 2022
de7405c
Fix pairwise bound comments
niklasschmitz Jun 7, 2022
86141e0
Fix comment
niklasschmitz Jun 7, 2022
aa2bf98
first batch of modifications
epolack Jun 8, 2022
afc7c14
some more
epolack Jun 8, 2022
a558373
Merge branch 'nfs/pairwise-bounds' into complex_pairwise
epolack Jun 8, 2022
f88dd86
workaround back with a vengeance
epolack Jun 8, 2022
1d1d790
bugfix
epolack Jun 8, 2022
610d404
Update forwarddiff_rules.jl
epolack Jun 9, 2022
284a915
Antoine's comment
epolack Jun 13, 2022
82c3bda
Merge remote-tracking branch 'origin/master' into complex_pairwise
epolack Jun 13, 2022
74908cf
bugfix
epolack Jun 13, 2022
0927546
Update phonons.jl
epolack Jun 28, 2022
1c3d395
Update pairwise.jl
epolack Jun 28, 2022
9200787
Merge remote-tracking branch 'origin/master' into complex_pairwise
epolack Jun 28, 2022
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
testing ph_disp
  • Loading branch information
epolack committed Jun 7, 2022
commit d4d48e934a24e3d294c0d82ccc422fa18e527b7e
137 changes: 137 additions & 0 deletions test/phonons.jl
Original file line number Diff line number Diff line change
@@ -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
5 changes: 5 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -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