Merge pull request #6 from HPMolSim/cutoff

add changeable short range cutoff
HPMolSim · Apr 19, 2024 · 4194286 · 4194286
2 parents 0f20d4e + 13b91b8
commit 4194286
Show file tree

Hide file tree

Showing 7 changed files with 79 additions and 74 deletions.
diff --git 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
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git 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
-Original file 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 @@