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.

Sign up for GitHub

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

add changeable short range cutoff #6

Merged
merged 2 commits into from
Apr 19, 2024
Merged

add changeable short range cutoff #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d229995
adding changable rc
ArrogantGao Apr 15, 2024
13b91b8
add changable r_c
ArrogantGao Apr 19, 2024
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
6 changes: 3 additions & 3 deletions src/FSSoGInteraction.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ function FSSoG_naive(L::NTuple{3, T}, n_atoms::Int, r_c::T, k_c::T; preset::Int

@assert preset ≤ length(preset_parameters)

b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]), T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]) * r_c, T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
uspara = USeriesPara(b, σ, ω, M)
k_set = k_set_2D(k_c, L)
return FSSoG_naive{T}(b, σ, ω, M, r_c, k_c, ϵ, L, k_set, uspara, n_atoms)
Expand Down Expand Up @@ -78,7 +78,7 @@ function FSSoGInteraction(

@assert preset ≤ length(preset_parameters)

b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]), T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]) * r_c, T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])

return FSSoGInteraction(b, σ, ω, M, L, n_atoms, r_c, Q_short, r_min, N_real_mid, w, β, extra_pad_ratio, cheb_order, M_mid, N_grid_long, Q_long, Rz_0, Q_0, ϵ = ϵ)
end
Expand Down Expand Up @@ -117,7 +117,7 @@ function FSSoGThinInteraction(

@assert preset ≤ length(preset_parameters)

b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]), T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]) * r_c, T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])

return FSSoGThinInteraction(b, σ, ω, M, L, n_atoms, r_c, Q_short, r_min, N_real, R_z, w, β, cheb_order, Taylor_Q, Rz_0, Q_0, ϵ = ϵ)
end
16 changes: 8 additions & 8 deletions src/U_series.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,12 @@ function BSA(x::T, b::T, σ::T, M::Int) where{T}
end

const preset_parameters = [
tuple(2.0, 5.027010924194599, 0.994446492762232252, 18),
tuple(1.62976708826776469, 3.633717409009413, 1.00780697934380681, 39),
tuple(1.48783512395703226, 2.662784519725113, 0.991911705759818, 45),
tuple(1.32070036405934420, 2.277149356440992, 1.00188914114811980, 78),
tuple(1.21812525709410644, 1.774456369233284, 1.00090146156033341, 144),
tuple(1.14878150173321925, 1.370684570284264 , 1.000036834835822, 233)]
tuple(2.0, 0.5027010924194599, 0.994446492762232252, 18),
tuple(1.62976708826776469, 0.3633717409009413, 1.00780697934380681, 39),
tuple(1.48783512395703226, 0.2662784519725113, 0.991911705759818, 45),
tuple(1.32070036405934420, 0.2277149356440992, 1.00188914114811980, 78),
tuple(1.21812525709410644, 0.1774456369233284, 1.00090146156033341, 144),
tuple(1.14878150173321925, 0.1370684570284264 , 1.000036834835822, 233)]

function U_series(x::T, b::T, σ::T, M::Int) where{T}
U_series_value = zero(T)
Expand All @@ -41,10 +41,10 @@ function USeriesPara(b::T, σ::T, ω::T, M::Int) where{T}
return USeriesPara(sw)
end

function USeriesPara(preset::Int, T::DataType = Float64)
function USeriesPara(preset::Int; r_c::T = 10.0) where{T}
@assert preset ≤ length(preset_parameters)

b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]), T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]) * r_c, T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])

return USeriesPara(b, σ, ω, M)
end
Expand Down
8 changes: 4 additions & 4 deletions src/energy/energy_short.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,10 @@ function F0_cal(b::T, σ::T, ω::T, M::Int) where{T}
return - log(b) / sqrt(2π) / σ * (ω + (one(T) - b^(-M)) / (b - one(T)))
end

function F0_cal(;preset::Int = 1, T::DataType = Float64)
function F0_cal(;preset::Int = 1, r_c::T = 10.0) where{T}
@assert preset ≤ length(preset_parameters)

b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]), T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])
b, σ, ω, M = T(preset_parameters[preset][1]), T(preset_parameters[preset][2]) * r_c, T(preset_parameters[preset][3]), Int(preset_parameters[preset][4])

return F0_cal(b, σ, ω, M)
end
Expand All @@ -18,9 +18,9 @@ function Es_USeries_Cheb(uspara::USeriesPara{T}, r_min::T, r_max::T, Q::Int) whe
end

function Es_Cheb_precompute(preset::Int, r_min::T, r_max::T, Q::Int) where{T}
uspara = USeriesPara(preset, T)
uspara = USeriesPara(preset, r_c = r_max)
uspara_cheb = Es_USeries_Cheb(uspara, r_min, r_max, Q)
F0 = F0_cal(preset = preset, T = T)
F0 = F0_cal(preset = preset, r_c = r_max)
return uspara_cheb, F0
end

Expand Down
2 changes: 2 additions & 0 deletions src/types.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@ struct USeriesPara{T}
sw::Vector{NTuple{2, T}}
end

Base.show(io::IO, uspara::USeriesPara{T}) where{T} = print(io, "USeriesPara($T), M=$(length(uspara.sw))")

struct FSSoG_naive{T} <: ExTinyMD.AbstractInteraction
b::T
σ::T
Expand Down
4 changes: 2 additions & 2 deletions test/U_series.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
@testset "compare the BSA formula" begin
@info "testing the BSA formula"
for (preset, accuracy) in zip([1, 2, 3, 4, 5], [1e-1, 1e-3, 1e-4, 1e-6, 1e-8])
for (preset, accuracy) in zip([1, 2, 3, 4, 5, 6], [1e-1, 1e-3, 1e-4, 1e-6, 1e-8, 1e-12])
@testset "preset = $preset" begin
b, σ, ω, M = FastSpecSoG.preset_parameters[preset]
for x in range(0.2, 10.0, length=100)
for x in range(1.0, 10.0, length=100)
@test abs(BSA(x, b, σ, M) * x - 1.0) < accuracy
end
end
Expand Down
96 changes: 49 additions & 47 deletions test/energy.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
@info "testing the energy cube"

n_atoms = 100
L = 100.0
L = 50.0
boundary = ExTinyMD.Q2dBoundary(L, L, L)

atoms = Vector{Atom{Float64}}()
Expand All @@ -17,36 +17,37 @@

info = SimulationInfo(n_atoms, atoms, (0.0, L, 0.0, L, 0.0, L), boundary; min_r = 1.0, temp = 1.0)

Ewald2D_interaction = Ewald2DInteraction(n_atoms, 5.0, 0.2, (L, L, L), ϵ = 1.0)
Ewald2D_interaction = Ewald2DInteraction(n_atoms, 5.0, 0.25, (L, L, L), ϵ = 1.0)
Ewald2D_neighbor = CellList3D(info, Ewald2D_interaction.r_c, boundary, 1)
energy_ewald = energy(Ewald2D_interaction, Ewald2D_neighbor, info, atoms)

fssog_naive = FSSoG_naive((L, L, L), n_atoms, 10.0, 2.0, preset = 2)
fssog_neighbor = CellList3D(info, fssog_naive.r_c, boundary, 1)
energy_sog_naive = energy_naive(fssog_naive, fssog_neighbor, info, atoms)
for r_c in [10.0, 15.0]
fssog_naive = FSSoG_naive((L, L, L), n_atoms, r_c, 4.0, preset = 3)
fssog_neighbor = CellList3D(info, fssog_naive.r_c, boundary, 1)
energy_sog_naive = energy_naive(fssog_naive, fssog_neighbor, info, atoms)

N_real = (128, 128, 128)
w = (16, 16, 16)
β = 5.0 .* w
extra_pad_ratio = 2
cheb_order = 10
preset = 2
M_mid = 3
N_real = (128, 128, 128)
w = (16, 16, 16)
β = 5.0 .* w
extra_pad_ratio = 2
cheb_order = 10
preset = 3
M_mid = 3

N_grid = (32, 32, 32)
Q = 24
R_z0 = 32
Q_0 = 32
r_c = 10.0
N_grid = (32, 32, 32)
Q = 48
R_z0 = 32
Q_0 = 32

fssog_interaction = FSSoGInteraction((L, L, L), n_atoms, r_c, Q, 0.5, N_real, w, β, extra_pad_ratio, cheb_order, M_mid, N_grid, Q, R_z0, Q_0; preset = preset, ϵ = 1.0)
fssog_interaction = FSSoGInteraction((L, L, L), n_atoms, r_c, Q, 0.5, N_real, w, β, extra_pad_ratio, cheb_order, M_mid, N_grid, Q, R_z0, Q_0; preset = preset, ϵ = 1.0)

fssog_neighbor = CellList3D(info, fssog_interaction.r_c, fssog_interaction.boundary, 1)
energy_sog = ExTinyMD.energy(fssog_interaction, fssog_neighbor, info, atoms)
fssog_neighbor = CellList3D(info, fssog_interaction.r_c, fssog_interaction.boundary, 1)
energy_sog = ExTinyMD.energy(fssog_interaction, fssog_neighbor, info, atoms)

@test abs(energy_ewald - energy_sog_naive) < 1e-4
@test abs(energy_ewald - energy_sog) < 1e-4
@test abs(energy_sog_naive - energy_sog) < 1e-6
@test abs(energy_ewald - energy_sog_naive) < 1e-4
@test abs(energy_ewald - energy_sog) < 1e-4
@test abs(energy_sog_naive - energy_sog) < 1e-6
end
end

@testset "energy thin" begin
Expand Down Expand Up @@ -74,28 +75,29 @@ end
Ewald2D_neighbor = CellList3D(info, Ewald2D_interaction.r_c, boundary, 1)
energy_ewald = energy(Ewald2D_interaction, Ewald2D_neighbor, info, atoms)

N_real = (128, 128)
R_z = 32
w = (16, 16)
β = 5.0 .* w
cheb_order = 16
preset = 2
Q = 24
Q_0 = 32
R_z0 = 32
Taylor_Q = 24
r_c = 10.0

fssog_interaction = FSSoGThinInteraction((Lx, Ly, Lz), n_atoms, r_c, Q, 0.5, N_real, R_z, w, β, cheb_order, Taylor_Q, R_z0, Q_0; preset = preset, ϵ = 1.0)

fssog_neighbor = CellList3D(info, fssog_interaction.r_c, fssog_interaction.boundary, 1)
energy_sog = ExTinyMD.energy(fssog_interaction, fssog_neighbor, info, atoms)

fssog_naive = FSSoG_naive((Lx, Ly, Lz), n_atoms, 10.0, 2.0, preset = 2)
fssog_neighbor = CellList3D(info, fssog_naive.r_c, boundary, 1)
energy_sog_naive = energy_naive(fssog_naive, fssog_neighbor, info, atoms)

@test abs(energy_ewald - energy_sog) < 1e-3
@test abs(energy_ewald - energy_sog_naive) < 1e-3
@test abs(energy_sog - energy_sog_naive) < 1e-6
for r_c in [10.0, 15.0]
N_real = (128, 128)
R_z = 32
w = (16, 16)
β = 5.0 .* w
cheb_order = 16
preset = 3
Q = 48
Q_0 = 32
R_z0 = 32
Taylor_Q = 24

fssog_interaction = FSSoGThinInteraction((Lx, Ly, Lz), n_atoms, r_c, Q, 0.5, N_real, R_z, w, β, cheb_order, Taylor_Q, R_z0, Q_0; preset = preset, ϵ = 1.0)

fssog_neighbor = CellList3D(info, fssog_interaction.r_c, fssog_interaction.boundary, 1)
energy_sog = ExTinyMD.energy(fssog_interaction, fssog_neighbor, info, atoms)

fssog_naive = FSSoG_naive((Lx, Ly, Lz), n_atoms, r_c, 3.0, preset = 3)
fssog_neighbor = CellList3D(info, fssog_naive.r_c, boundary, 1)
energy_sog_naive = energy_naive(fssog_naive, fssog_neighbor, info, atoms)

@test abs(energy_ewald - energy_sog) < 1e-3
@test abs(energy_ewald - energy_sog_naive) < 1e-3
@test abs(energy_sog - energy_sog_naive) < 1e-6
end
end
21 changes: 11 additions & 10 deletions test/energy_short.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,17 +37,18 @@ end

for preset in 1:5
@testset "preset = $preset" begin
uspara = USeriesPara(preset)
r_min = 0.5
r_max = 10.0
Q = 64
uspara_cheb, F0 = Es_Cheb_precompute(preset, r_min, r_max, Q)
interaction = FSSoG_naive((L, L, L), n_atoms, r_max, 3.0, preset = preset)
neighbor = CellList3D(info, interaction.r_c, boundary, 1)
for r_c in range(1.0, 1.0, 10)
r_min = 0.5
r_max = r_c
Q = 64
uspara_cheb, F0 = Es_Cheb_precompute(preset, r_min, r_max, Q)
interaction = FSSoG_naive((L, L, L), n_atoms, r_max, 3.0, preset = preset)
neighbor = CellList3D(info, interaction.r_c, boundary, 1)

Es_naive = short_energy_naive(interaction, neighbor, position, charge)
Es_Cheb = short_energy_Cheb(uspara_cheb, interaction.r_c, F0, boundary, neighbor, position, charge)
@test isapprox(Es_naive, Es_Cheb; atol = 1e-12)
Es_naive = short_energy_naive(interaction, neighbor, position, charge)
Es_Cheb = short_energy_Cheb(uspara_cheb, interaction.r_c, F0, boundary, neighbor, position, charge)
@test isapprox(Es_naive, Es_Cheb; atol = 1e-12)
end
end
end
end
Loading