Skip to content

Commit

Permalink
stuff
Browse files Browse the repository at this point in the history
  • Loading branch information
@ejmeitz
ejmeitz committed Oct 17, 2023
1 parent 8e5d2db commit 4c05d88
Show file tree
Hide file tree
Showing 9 changed files with 407 additions and 293 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
cuTENSOR = "011b41b2-24ef-40a8-b3eb-fa098493e9e1"

[compat]
SimpleCrystals = "0.1.6"
SimpleCrystals = "0.1.7"
TensorOperations = "4.0.0"
julia = "1"

Expand Down
60 changes: 44 additions & 16 deletions src/autodiff_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,24 +97,16 @@ function third_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol) w
end
end

#Acoustic Sum Rule #&TODO
# Threads.@threads for i in range(1, N_atoms) # index of block matrix
# for α in range(1,D)
# for β in range(1,D)
# ii = D*(i-1) + α
# jj = D*(i-1) + β # i == j because we're on diagonal
# IFC2[ii,jj] = -1*sum(IFC2[ii, β:D:end])
# end
# end
# end
#Acoustic Sum Rule
ASR!(IFC3, N_atoms, D)

IFC3 = apply_tols!(IFC3,tol)

return ThirdOrderMatrix(IFC3, unit(pot.ϵ / pot.σ^3), tol)

end

function fourth_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol; float_type = Float32) where D
function fourth_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol) where D
vars = make_variables(:r, D)
r_norm = sqrt(sum(x -> x^2, vars))
pot_symbolic = potential_nounits(pot, r_norm)
Expand All @@ -124,7 +116,7 @@ function fourth_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol;
FO_exec = make_function(fourth_order_symbolic, vars)

N_atoms = n_atoms(sys)
IFC4 = FC_val{float_type,4}[]
IFC4 = FC_val{Float64,4}[]

for i in range(1, N_atoms)
for j in range(i + 1, N_atoms)
Expand All @@ -143,10 +135,10 @@ function fourth_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol;
#iiij
ii = D*(i-1) + α; jj = D*(i-1) + β; kk = D*(i-1) + γ; ll = D*(j-1) + δ
#& THERE WILL BE DUPLICATES SINCE THIS USES PUSH INSTEAD OF JUST OVERWRITTING DATA!!!
set_terms_fourth_order!(χ, ii, jj, kk, ll, float_type(iiij_block[α,β,γ,δ]))
set_terms_fourth_order!(χ, ii, jj, kk, ll, iiij_block[α,β,γ,δ])
#iijj
ii = D*(i-1) + α; jj = D*(i-1) + β; kk = D*(j-1) + γ; ll = D*(j-1) + δ
set_terms_fourth_order!(χ, ii, jj, kk, ll, float_type(iiij_block[α,β,γ,δ]))
set_terms_fourth_order!(χ, ii, jj, kk, ll, iiij_block[α,β,γ,δ])
end
end
end
Expand All @@ -167,7 +159,7 @@ function fourth_order_AD_test(sys::SuperCellSystem{D}, pot::PairPotential, tol;
if abs(iiij_block[α,β,γ,δ]) > tol
#ijjj
ii = D*(i-1) + α; jj = D*(j-1) + β; kk = D*(j-1) + γ; ll = D*(j-1) + δ
set_terms_fourth_order!(χ, ii, jj, kk, ll, float_type(ijjj_block[α,β,γ,δ]))
set_terms_fourth_order!(χ, ii, jj, kk, ll, ijjj_block[α,β,γ,δ])
end
end
end
Expand Down Expand Up @@ -219,4 +211,40 @@ function energy_loop_mcc(sys::SuperCellSystem{D}, pot::PairPotential, q, phi, ma

return U_total

end
end

#q_vars = make_variable(:q, N_modes)
#r_vars = zeros(eltype(q_vars), length(q_vars))
#for i in 1:N_modes
# r_vars[i] = dot(q_vars, phi[i,:])
#end

# divide r_vars by mass_sqrt

# pass r_vars into energy_loop2

# function energy_loop2(pot::PairPotential, posns, r_cut, box_sizes, N_atoms)

# U_total = 0.0
# box_sizes = ustrip.(box_sizes)

# for i in range(1,N_atoms)
# for j in range(i+1, N_atoms)
# r = posns[3*(i-1) + 1: 3*i] .- posns[3*(j-1) + 1: 3*j ]
# nearest_mirror_AD!(r, box_sizes)
# dist = sqrt(sum(x-> x*x, r))
# # if dist < r_cut
# U_total += potential_nounits(pot, dist)
# # end
# end
# end

# return U_total

# end

# #& not quite right but I just wanna test the AD
# function nearest_mirror_AD!(r_ij, box_sizes)
# r_ij .+= sign.(r_ij .- box_sizes./2) .* box_sizes
# return r_ij
# end
58 changes: 17 additions & 41 deletions src/finite_diff_ifc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,7 @@ function second_order_finite_diff(sys_eq::SuperCellSystem{3}, pot::PairPotential

energy_unit = zero(potential(pot,1u"Å"))

#Make mutable #& change SimpleCrystals to not use SVector
posns = [Vector(a) for a in positions(sys_eq)]
posns = positions(sys_eq)

if atom_idxs[1] != atom_idxs[2]

Expand Down Expand Up @@ -61,25 +60,23 @@ function third_order_finite_diff(sys_eq::SuperCellSystem{3}, pot::PairPotential,

@assert length(atom_idxs) == 3
@assert length(cartesian_idxs) == 3
@assert !all(atom_idxs .== atom_idxs[1]) "Cannot check self terms with finite difference"

N_atoms = n_atoms(sys)
N_atoms = n_atoms(sys_eq)

z = zero(pot.σ)
combos = [[z,h,-h],[z,h,h],[z,-h,h],[z,-h,-h]]

force_unit = zero(force(pot,1u"Å"))
forces = zeros(4)*force_unit

#Make mutable #& change SimpleCrystals to not use SVector
posns = [Vector(a) for a in positions(sys)]
posns = positions(sys_eq)

for (c,combo) in enumerate(combos)

posns[atom_idxs[2]][cartesian_idxs[2]] += combo[2]
posns[atom_idxs[3]][cartesian_idxs[3]] += combo[3]

forces[c] = force_loop_j(pot, posns, force_unit, r_cut, sys.box_sizes_SC, N_atoms, atom_idxs[1], cartesian_idxs[1])
forces[c] = force_loop_j(pot, posns, force_unit, r_cut, sys_eq.box_sizes_SC, N_atoms, atom_idxs[1], cartesian_idxs[1])

posns[atom_idxs[2]][cartesian_idxs[2]] -= combo[2]
posns[atom_idxs[3]][cartesian_idxs[3]] -= combo[3]
Expand All @@ -98,16 +95,15 @@ function fourth_order_finite_diff(sys_eq::SuperCellSystem{3}, pot::PairPotential
@assert length(cartesian_idxs) == 4
@assert !all(atom_idxs .== atom_idxs[1]) "Cannot check self terms with finite difference"

N_atoms = n_atoms(sys)
N_atoms = n_atoms(sys_eq)

combos = [[h,h,h],[h,h,-h],[h,-h,h],[h,-h,-h],[-h,h,h],[-h,h,-h],[-h,-h,h],[-h,-h,-h]]

force_unit = zero(force(pot,1u"Å"))
forces = zeros(8)*force_unit

#Make mutable #& change SimpleCrystals to not use SVector
posns = [Vector(a) for a in positions(sys)]

posns = positions(sys_eq)

for (c,combo) in enumerate(combos)

posns[atom_idxs[1]][cartesian_idxs[1]] += combo[1]
Expand All @@ -120,7 +116,7 @@ function fourth_order_finite_diff(sys_eq::SuperCellSystem{3}, pot::PairPotential
for i in range(1,N_atoms)
if i != l
r = posns[i] .- posns[l]
nearest_mirror!(r, sys.box_sizes_SC)
nearest_mirror!(r, sys_eq.box_sizes_SC)
dist = norm(r)

#Make sure mirrored particle is in cuttoff
Expand Down Expand Up @@ -148,48 +144,28 @@ function fourth_order_finite_diff(sys_eq::SuperCellSystem{3}, pot::PairPotential
return (1/(8*(h^3)))*(forces[1] - forces[2] - forces[3] + forces[4] - forces[5] + forces[6] + forces[7] - forces[8])
end

#& CURRENTLY DIFFERS by -1 FROM MINE
function check_ifc(sys::SuperCellSystem{D}, ifc::DenseForceConstants{O}, pot::PairPotential, n_points;
r_cut = pot.r_cut, h = 0.04*0.5291772109u"Å") where {D,O}

N_atoms = n_atoms(sys)
ri = rand(1:N_atoms, (n_points,O)) #random indexes
ri = rand(1:N_atoms, (n_points,O)) #random atom indexes
rci = rand(1:D, (n_points, O)) #random cartesian indexes

finite_diff_funcs = [second_order_finite_diff, third_order_finite_diff, fourth_order_finite_diff]
finite_diff_func = finite_diff_funcs[O-1]

ifc_fd = zeros(n_points)*ifc.units
#& ONLY WORKS FOR SPECIFIC DIMENSION
ifc_actual = [ifc[D*(ri[i,1] - 1) + rci[i,1], D*(ri[i,2] - 1) + rci[i,2]] for i in 1:n_points]
# ifc_actual = [ifc[D*(ri[i,1] - 1) + rci[i,1], D*(ri[i,2] - 1) + rci[i,2], D*(ri[i,3] - 1) + rci[i,3]]
# for i in 1:n_points]
idxs = (ri_row, rci_row) -> [D*(ri_row[o] - 1) + rci_row[o] for o in 1:O]
ifc_actual = @views [ifc[idxs(ri[i,:], rci[i,:])...] for i in 1:n_points]

ifc_fd = zeros(n_points)*ifc.units
for i in 1:n_points
ifc_fd[i] = finite_diff_func(sys, pot, ri[i,:], rci[i,:]; r_cut = r_cut, h = h)
val = finite_diff_func(sys, pot, ri[i,:], rci[i,:]; r_cut = r_cut, h = h)
if ustrip.(abs(val)) > ifc.tol
ifc_fd[i] = val
end
end

return ifc_fd, (ifc_actual.*ifc.units)

end

function check_ifc(sys::SuperCellSystem{D}, ifc::SparseForceConstants{O}, pot::PairPotential, n_points;
r_cut = pot.r_cut, h = 0.04*0.5291772109u"Å") where {D,O}

random_idxs = rand(1:length(ifc), (n_points,))

finite_diff_funcs = [second_order_finite_diff, third_order_finite_diff, fourth_order_finite_diff]
finite_diff_func = finite_diff_funcs[O-1]

ifc_fd = zeros(n_points)*ifc.units
ifc_actual = value.(ifc[random_idxs])

for p in 1:n_points
idxs = idx(ifc[random_idxs[p]])
cart_idxs = ((idxs .- 1) .% D) .+ 1
atom_idxs = ((idxs .- cart_idxs) ./ D) .- 1
ifc_fd[p] = finite_diff_func(sys, pot, atom_idxs, cart_idxs; r_cut = r_cut, h = h)
end

return ifc_fd, (ifc_actual.*ifc.units)

end
16 changes: 8 additions & 8 deletions src/fourth_order.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ i = j and k = l (e.g. (i,i,j,j)) and terms where i = j = k (e.g. (i,i,i,j)). All
more than 2 unique indices will be 0.
"""
function fourth_order_sparse(sys::SuperCellSystem{D}, pot::PairPotential,
tol; float_type = Float32, r_cut = pot.r_cut, nthreads::Integer = Threads.nthreads()) where D
tol; r_cut = pot.r_cut, nthreads::Integer = Threads.nthreads()) where D

@assert all(r_cut .< sys.box_sizes_SC) "Cutoff larger than L/2"

Expand All @@ -21,7 +21,7 @@ function fourth_order_sparse(sys::SuperCellSystem{D}, pot::PairPotential,
#Create storage for each thread
χ_arrays = MultiVectorStorage(FC_val{float_type,4}, nthreads)

atoms_per_thd = cld(N_atoms,nthreads) #* DONT LIKE THIS, BAD BALANCING
atoms_per_thd = cld(N_atoms,nthreads) #* DONT LIKE THIS, BAD BALANCING, USE ITERATORS.PARTITION

Base.@sync for thd in 1:nthreads
Base.Threads.@spawn begin
Expand All @@ -43,7 +43,7 @@ function fourth_order_sparse(sys::SuperCellSystem{D}, pot::PairPotential,
for β in range(1,D)
for γ in range(1,D)
for δ in range(1,D)
val = float_type(ustrip(ϕ₄(pot, dist, r, α, β, γ, δ)))
val = ustrip(ϕ₄(pot, dist, r, α, β, γ, δ))

if abs(val) > tol
# i,i,j,j terms
Expand Down Expand Up @@ -73,7 +73,7 @@ function fourth_order_sparse(sys::SuperCellSystem{D}, pot::PairPotential,
for β in range(1,D)
for γ in range(1,D)
for δ in range(1,D)
val = -float_type(ustrip(ϕ₄(pot, dist, r, α, β, γ, δ)))
val = -ustrip(ϕ₄(pot, dist, r, α, β, γ, δ))

if abs(val) > tol
# j,j,j,i terms
Expand All @@ -97,16 +97,16 @@ function fourth_order_sparse(sys::SuperCellSystem{D}, pot::PairPotential,
end

#Concat arrays calculated by separate threads
χ = convert(Vector{FC_val{float_type, 4}}, χ_arrays)
χ = convert(Vector{FC_val{Float64, 4}}, χ_arrays)

#Acoustic Sum Rule (i,i,i,i terms) #& how to parallelize this part?
for i in range(1,N_atoms)
#Acoustic Sum Rule (i,i,i,i terms)
Threads.@threads for i in range(1,N_atoms)
for α in range(1,D)
for β in range(1,D)
for γ in range(1,D)
for δ in range(1,D)
ii = D*(i-1) + α; jj = D*(i-1) + β; kk = D*(i-1) + γ; ll = D*(i-1) + δ
val = float_type(ϕ₄_self(sys, pot, i, α, β, γ, δ, r_cut))
val = ϕ₄_self(sys, pot, i, α, β, γ, δ, r_cut)
if abs(val) > tol
push!(χ,FC_val(val, ii, jj, kk, ll))
end
Expand Down
Loading

0 comments on commit 4c05d88

Please sign in to comment.