diff --git a/src/FSSoGInteraction.jl b/src/FSSoGInteraction.jl index 6fd8cff..139e25d 100644 --- a/src/FSSoGInteraction.jl +++ b/src/FSSoGInteraction.jl @@ -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) @@ -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 @@ -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 \ No newline at end of file diff --git a/src/U_series.jl b/src/U_series.jl index c399c85..94f00cc 100644 --- a/src/U_series.jl +++ b/src/U_series.jl @@ -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) @@ -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 diff --git a/src/energy/energy_short.jl b/src/energy/energy_short.jl index c88fd5b..dbd93d1 100644 --- a/src/energy/energy_short.jl +++ b/src/energy/energy_short.jl @@ -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 @@ -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 diff --git a/src/types.jl b/src/types.jl index d547087..905e28d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -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 diff --git a/test/U_series.jl b/test/U_series.jl index 4adac0b..21ab4e6 100644 --- a/test/U_series.jl +++ b/test/U_series.jl @@ -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 diff --git a/test/energy.jl b/test/energy.jl index 8a57e27..53b126a 100644 --- a/test/energy.jl +++ b/test/energy.jl @@ -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}}() @@ -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 @@ -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 \ No newline at end of file diff --git a/test/energy_short.jl b/test/energy_short.jl index 3e6bf0b..fa0b12b 100644 --- a/test/energy_short.jl +++ b/test/energy_short.jl @@ -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 \ No newline at end of file