rewrite the zeroth order

Now the 0th order term is calculated via mixed talyor-chebyshev FGT

src/FFCT/precompute.jl

@@ -4,8 +4,17 @@
     return fill!(a, zero(T))
 end
 
+function zero_paras_gen(L_z::T, Q_0::Int) where{T}
+
+    r_z0 = Vector{T}([L_z / 2 * ( 1.0 - cos((2i - 1)*π / 2Q_0) ) for i in 1:Q_0])
+    grids0 = zeros(T, Q_0)
+    chebcoefs0 = zeros(T, Q_0)
+
+    return r_z0, grids0, chebcoefs0
+end
+
 # precompute the parameters to be used
-function long_paras_gen(L::NTuple{3, T}, N_grid::NTuple{3, Int}, n_atoms::Int) where{T}
+function long_paras_gen(L::NTuple{3, T}, N_grid::NTuple{3, Int}) where{T}
     L_x, L_y, L_z = L
     N_x, N_y, N_z = N_grid
     k_x = Vector{T}([2π * i / L_x for i in -N_x:N_x])
@@ -19,11 +28,8 @@ function long_paras_gen(L::NTuple{3, T}, N_grid::NTuple{3, Int}, n_atoms::Int) w
 
     phase_x = zeros(Complex{T}, 2N_x + 1)
     phase_y = zeros(Complex{T}, 2N_y + 1)
-
-    sort_z = zeros(Int, n_atoms)
-    z = zeros(T, n_atoms)
 
-    return k_x, k_y, r_z, H_r, H_c, phase_x, phase_y,  sort_z, z
+    return k_x, k_y, r_z, H_r, H_c, phase_x, phase_y
 end
 
 @inbounds function revise_phase_neg!(phase_x::AbstractArray{Complex{T}, 1}, phase_y::AbstractArray{Complex{T}, 1}, k_x::Vector{T}, k_y::Vector{T}, x::T, y::T) where{T}

src/FFCT/zeroth_order.jl

@@ -1,17 +1,45 @@
-function FGT1d(qs::Vector{T}, poses::Vector{NTuple{3, T}}, s::T, L_z::T, Q_0::Int) where{T}
+function proper_N(accuracy::T, η::T) where{T}
+    N = 0
+    t = one(T)
+    if η ≥ one(T)
+        return N
+    end
+
+    while t > η * accuracy
+        N += 1
+        t *= η / N
+    end
+    return N
+end
+
+function FGT1d(qs::Vector{T}, poses::Vector{NTuple{3, T}}, s::T, L_z::T, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}) where{T}
     n_atoms = length(qs)
-    grids = zeros(T, Q_0)
-    r_z = [L_z / 2 * ( 1.0 - cos((2i - 1)*π / (2 * Q_0)) ) for i in 1:Q_0]
+    Q_0 = length(r_z)
+    # grids = zeros(T, Q_0)
+    # r_z = [L_z / 2 * ( 1.0 - cos((2i - 1)*π / (2 * Q_0)) ) for i in 1:Q_0]
+
+    set_zeros!(grids)
+    set_zeros!(chebcoefs)
+
+    # (L_z/s)^2N / factorials(N) < 1e-16
+    N = proper_N(1e-16, (L_z / s)^2)
 
     for i in 1:n_atoms
         q = qs[i]
         for j in 1:Q_0
             z = poses[i][3]
-            grids[j] += q * exp(- (z - r_z[j])^2 / s^2)
+            t = - (z - r_z[j])^2 / s^2
+            if N == 0
+                grids[j] += q * exp(t)
+            else
+                for k in 1:N
+                    grids[j] += q * special_powern(t, k)
+                end
+            end
         end
     end
 
-    chebcoefs = zeros(T, Q_0)
+    # chebcoefs = zeros(T, Q_0)
 
     for k in 1:Q_0, l in 1:Q_0
         cheb_temp = k == 1 ? 0.5 : chebpoly(k - 1, r_z[l] - L_z / T(2), L_z / T(2))
@@ -32,29 +60,29 @@ end
 
 function FGT1d_naive(qs::Vector{T}, poses::Vector{NTuple{3, T}}, s::T) where{T}
     n_atoms = length(qs)
-    E_k0 = zero(T)
+    E_k0 = big(zero(T))
 
     for i in 1:n_atoms, j in 1:n_atoms
         qi = qs[i]
         xi, yi, zi = poses[i]
         qj = qs[j]
         xj, yj, zj = poses[j]
-        E_k0 += qi * qj * exp( - (zi - zj)^2 / s^2)
+        E_k0 += qi * qj * exp( - big(zi - zj)^2 / s^2)
     end
 
-    return E_k0
+    return Float64(E_k0)
 end
 
 
-function zeroth_order(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, uspara::USeriesPara{T}, Q_0::Int) where{T}
+function zeroth_order(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, uspara::USeriesPara{T}, M_mid::Int, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}) where{T}
 
     L_x, L_y, L_z = L
 
     E_k0 = zero(T)
 
     for l in M_mid + 1:length(uspara.sw)
         sl, wl = uspara.sw[l]
-        E_k0 += wl * sl^2 * FGT1d(qs, poses, sl, L_z, Q_0)
+        E_k0 += wl * sl^2 * FGT1d(qs, poses, sl, L_z, r_z, chebcoefs, grids)
     end
 
     return π * E_k0 / (2 * L_x * L_y)

src/FSSoGInteraction.jl

@@ -42,7 +42,7 @@ function FSSoGInteraction(
     L::NTuple{3, T}, n_atoms::Int,
     r_c::T, Q_short::Int, r_min::T,
     N_real_mid::NTuple{3, Int}, w::NTuple{3, Int}, β::NTuple{3, T}, extra_pad_ratio::Int, cheb_order::Int, M_mid::Int, 
-    N_grid_long::NTuple{3, Int}, Q_long::Int, soepara::SoePara{Complex{T}}; 
+    N_grid_long::NTuple{3, Int}, Q_long::Int, Q_0::Int; 
     ϵ::T = one(T)) where{T}
 
     uspara = USeriesPara(b, σ, ω, M)
@@ -58,31 +58,33 @@ function FSSoGInteraction(
     gridinfo, gridbox, cheb_coefs, scalefactor = mid_paras_gen(N_real_mid, w, β, L, extra_pad_ratio, cheb_order, uspara, M_mid)
 
     # some of these term are for other methods, not used
-    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y, sort_z, z = long_paras_gen(L, N_grid_long, n_atoms)
+    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y = long_paras_gen(L, N_grid_long)
     cheb_mat = ChebUSeries(k_x, k_y, L[3], uspara, M_mid, Q_long)
 
-    return FSSoGInteraction{T}(b, σ, ω, M, ϵ, L, boundary, n_atoms, uspara, position, charge, uspara_cheb, r_c, F0, gridinfo, gridbox, cheb_coefs, scalefactor, M_mid, k_x, k_y, r_z, phase_x, phase_y, H_r, H_c, cheb_mat, z, sort_z, soepara)
+    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Q_0)
+
+    return FSSoGInteraction{T}(b, σ, ω, M, ϵ, L, boundary, n_atoms, uspara, position, charge, uspara_cheb, r_c, F0, gridinfo, gridbox, cheb_coefs, scalefactor, M_mid, k_x, k_y, r_z, phase_x, phase_y, H_r, H_c, cheb_mat, Q_0, r_z0, chebcoefs0, grids0)
 end
 
 function FSSoGInteraction(
     L::NTuple{3, T}, n_atoms::Int,
     r_c::T, Q_short::Int, r_min::T,
     N_real_mid::NTuple{3, Int}, w::NTuple{3, Int}, β::NTuple{3, T}, extra_pad_ratio::Int, cheb_order::Int, M_mid::Int, 
-    N_grid_long::NTuple{3, Int}, Q_long::Int, soepara::SoePara{Complex{T}}; 
+    N_grid_long::NTuple{3, Int}, Q_long::Int, Q_0::Int; 
     preset::Int = 1, ϵ::T = one(T)) 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])
 
-    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, soepara, ϵ = ϵ)
+    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, Q_0, ϵ = ϵ)
 end
 
 function FSSoGThinInteraction(
     b::T, σ::T, ω::T, M::Int,
     L::NTuple{3, T}, n_atoms::Int,
     r_c::T, Q_short::Int, r_min::T,
-    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int, soepara::SoePara{Complex{T}}; 
+    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int, Q_0::Int; 
     ϵ::T = one(T)) where{T}
 
     uspara = USeriesPara(b, σ, ω, M)
@@ -96,21 +98,20 @@ function FSSoGThinInteraction(
 
     gridinfo, pad_grids, cheb_coefs, scalefactors, H_r, H_c, cheb_value, r_z = thin_paras_gen(N_real, R_z, w, β, L, cheb_order, uspara, Taylor_Q)
 
-    sort_z = zeros(Int, n_atoms)
-    z = zeros(T, n_atoms)
+    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Q_0)
 
-    return FSSoGThinInteraction{T}(b, σ, ω, M, ϵ, L, boundary, n_atoms, uspara, position, charge, uspara_cheb, r_c, F0, r_z, H_r, H_c, gridinfo, pad_grids, scalefactors, cheb_coefs, cheb_value, z, sort_z, soepara)
+    return FSSoGThinInteraction{T}(b, σ, ω, M, ϵ, L, boundary, n_atoms, uspara, position, charge, uspara_cheb, r_c, F0, r_z, H_r, H_c, gridinfo, pad_grids, scalefactors, cheb_coefs, cheb_value, Q_0, r_z0, chebcoefs0, grids0)
 end
 
 function FSSoGThinInteraction(
     L::NTuple{3, T}, n_atoms::Int,
     r_c::T, Q_short::Int, r_min::T,
-    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int, soepara::SoePara{Complex{T}}; 
+    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int, Q_0::Int; 
     preset::Int = 1, ϵ::T = one(T)) 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])
 
-    return FSSoGThinInteraction(b, σ, ω, M, L, n_atoms, r_c, Q_short, r_min, N_real, R_z, w, β, cheb_order, Taylor_Q, soepara, ϵ = ϵ)
+    return FSSoGThinInteraction(b, σ, ω, M, L, n_atoms, r_c, Q_short, r_min, N_real, R_z, w, β, cheb_order, Taylor_Q, Q_0, ϵ = ϵ)
 end

src/FastSpecSoG.jl

@@ -11,8 +11,9 @@ export short_energy_naive, long_energy_naive, energy_naive
 
 export Es_Cheb_precompute, short_energy_Cheb
 
-export mid_paras_gen, energy_mid, energy_long, energy_short
-export long_paras_gen, interpolate_nu_cheb!, interpolate_nu_loop!
+export mid_paras_gen, long_paras_gen, zero_paras_gen
+export energy_mid, energy_long, energy_short
+export interpolate_nu_cheb!, interpolate_nu_loop!
 export real2Cheb!, Cheb2real!, gather_nu
 export energy_long_loop_k, energy_long_cheb_k, energy_long_0
 

src/energy/energy_long.jl

@@ -1,13 +1,11 @@
 function energy_long_0(
     qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, M_mid::Int, 
-    z::Vector{T}, sort_z::Vector{Int}, 
-    uspara::USeriesPara{T}, soepara::SoePara{Complex{T}}
+    uspara::USeriesPara{T}, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}
     ) where{T}
 
     @assert M_mid ≤ length(uspara.sw)
 
-    revise_z!(z, sort_z, poses)
-    E_0 = zeroth_order(qs, z, soepara, uspara, sort_z, L, M_mid)
+    E_0 = zeroth_order(qs, poses, L, uspara, M_mid, r_z, chebcoefs, grids)
 
     return E_0
 end
@@ -16,7 +14,7 @@ function energy_long(interaction::FSSoGInteraction{T}) where{T}
 
     E_k = energy_long_cheb_k(interaction.charge, interaction.position, interaction.L, interaction.k_x, interaction.k_y, interaction.r_z, interaction.phase_x, interaction.phase_y, interaction.H_r, interaction.H_c, interaction.cheb_mat)
 
-    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, interaction.M_mid, interaction.z, interaction.sort_z, interaction.uspara, interaction.soepara)
+    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, interaction.M_mid, interaction.uspara, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
 
     return (E_k + E_0) / (4π * interaction.ϵ)
 end
@@ -25,7 +23,7 @@ function energy_long(interaction::FSSoGThinInteraction{T}) where{T}
 
     E_k = energy_long_thin_k(interaction.charge, interaction.position, interaction.L, interaction.r_z, interaction.H_r, interaction.H_c, interaction.gridinfo, interaction.pad_grids, interaction.scalefactors, interaction.cheb_coefs, interaction.cheb_value)
 
-    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, 0, interaction.z, interaction.sort_z, interaction.uspara, interaction.soepara)
+    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, 0, interaction.uspara, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
     return (E_k + E_0) / (4π * interaction.ϵ)
 end
 

src/energy/energy_long_naive.jl

@@ -79,19 +79,21 @@ function long_energy_sw_0(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{
 
     E = Atomic{T}(zero(T))
 
-    @threads for n in CartesianIndices((n_atoms, n_atoms))
-        i, j = Tuple(n)
+    @threads for i in 1:n_atoms
         qi = qs[i]
         xi, yi, zi = poses[i]
-        qj = qs[j]
-        xj, yj, zj = poses[j]
+        t = big(zero(T))
+        for j in 1:n_atoms
+            qj = qs[j]
+            xj, yj, zj = poses[j]
 
-        ϕ_ij = w * s^2 * exp( - (zi - zj)^2 / s^2)
+            t += qj * exp( - big(zi - zj)^2 / s^2)
+        end
 
-        atomic_add!(E, qi * qj * ϕ_ij)
+        atomic_add!(E, qi * T(t))
     end
 
-    return E[] * π / (2 * L[1] * L[2])
+    return w * s^2 * E[] * π / (2 * L[1] * L[2])
 end
 
 function long_energy_us_k(qs::Vector{T}, poses::Vector{NTuple{3, T}}, cutoff::Int, L::NTuple{3, T}, uspara::USeriesPara{T}, M_min::Int, M_max::Int) where{T}

src/types.jl

@@ -53,9 +53,10 @@ mutable struct FSSoGInteraction{T} <: ExTinyMD.AbstractInteraction
     cheb_mat::Array{ChebPoly{1, T, T}, 2}
 
     # for FFCT zero order term
-    z::Vector{T}
-    sort_z::Vector{Int}
-    soepara::SoePara{Complex{T}}
+    Q_0::Int
+    r_z0::Vector{T}
+    chebcoefs0::Vector{T}
+    grids0::Vector{T}
 end
 
 mutable struct FSSoGThinInteraction{T} <: ExTinyMD.AbstractInteraction
@@ -87,8 +88,9 @@ mutable struct FSSoGThinInteraction{T} <: ExTinyMD.AbstractInteraction
     cheb_coefs::NTuple{2, ChebCoef{T}}
     cheb_value::Vector{Array{T, 1}}
 
-    # for FFCT zero order term
-    z::Vector{T}
-    sort_z::Vector{Int}
-    soepara::SoePara{Complex{T}}
+    # for zero order term
+    Q_0::Int
+    r_z0::Vector{T}
+    chebcoefs0::Vector{T}
+    grids0::Vector{T}
 end

test/Chebyshev.jl

@@ -21,7 +21,7 @@ end
     uspara = USeriesPara(2)
     M_mid = 0
 
-    k_x, k_y, r_z, H_r1, H_c, phase_x, phase_y,  sort_z, z = long_paras_gen(L, N_grid, n_atoms)
+    k_x, k_y, r_z, H_r1, H_c, phase_x, phase_y = long_paras_gen(L, N_grid)
     cheb_mat = ChebUSeries(k_x, k_y, L[3], uspara, M_mid, 64)
 
     H_r2 = copy(H_r1)

test/energy.jl

@@ -35,9 +35,10 @@
 
     N_grid = (32, 32, 32)
     Q = 24
+    Q_0 = 32
     r_c = 10.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, SoePara16(); 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, 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)
@@ -79,10 +80,11 @@ end
     cheb_order = 16
     preset = 2
     Q = 24
+    Q_0 = 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, SoePara16(); preset = preset, ϵ = 1.0)
+    fssog_interaction = FSSoGThinInteraction((Lx, Ly, Lz), n_atoms, r_c, Q, 0.5, N_real, R_z, w, β, cheb_order, Taylor_Q, 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)

test/energy_long.jl

@@ -9,16 +9,16 @@
 
     N_grid = (32, 32, 32)
     uspara = USeriesPara(2)
-    soepara = SoePara16()
     M_mid = 5
 
-    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y, sort_z, z = long_paras_gen(L, N_grid, n_atoms)
+    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y = long_paras_gen(L, N_grid)
     cheb_mat = ChebUSeries(k_x, k_y, L[3], uspara, M_mid, 64)
+    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], 64)
 
     E_long_loop_k = energy_long_loop_k(qs, poses, L, M_mid, k_x, k_y, r_z, phase_x, phase_y, H_r, H_c, uspara)
     E_long_cheb_k = energy_long_cheb_k(qs, poses, L, k_x, k_y, r_z, phase_x, phase_y, H_r, H_c, cheb_mat)
 
-    E_long_0 = energy_long_0(qs, poses, L, M_mid, z, sort_z, uspara, soepara)
+    E_long_0 = energy_long_0(qs, poses, L, M_mid, uspara, r_z0, chebcoefs0, grids0)
 
     @info "running the direct summation for the long range part of the energy"
     # using the direct summation
@@ -76,7 +76,7 @@ end
     uspara = USeriesPara(1)
     M_mid = 0
 
-    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y,  sort_z, z = long_paras_gen(L, N_grid, n_atoms)
+    k_x, k_y, r_z, H_r, H_c, phase_x, phase_y = long_paras_gen(L, N_grid)
     cheb_mat = ChebUSeries(k_x, k_y, L[3], uspara, M_mid, 64)
 
     E_long_loop_k = energy_long_loop_k(qs, poses, L, M_mid, k_x, k_y, r_z, phase_x, phase_y, H_r, H_c, uspara)

test/interpolate.jl

@@ -11,11 +11,10 @@
     N_grid = (16, 16, 16)
     uspara = USeriesPara(2)
     Q = 64
-    soepara = SoePara16()
 
     for M_mid in 2:length(uspara.sw) - 1
         @testset "M_mid = $M_mid" begin
-            k_x, k_y, r_z, H_r, H_c, phase_x, phase_y,  sort_z, z = long_paras_gen(L, N_grid, n_atoms)
+            k_x, k_y, r_z, H_r, H_c, phase_x, phase_y = long_paras_gen(L, N_grid)
             cheb_mat = ChebUSeries(k_x, k_y, L[3], uspara, M_mid, Q)
 
             H_r_loop = copy(interpolate_nu_loop!(copy(H_r), qs, poses, L, k_x, k_y, phase_x, phase_y, r_z, uspara, M_mid))

test/runtests.jl

@@ -5,14 +5,14 @@ using Random
 Random.seed!(1234)
 
 @testset "FastSpecSoG.jl" begin
-    # include("U_series.jl")
-    # include("energy.jl")
-    # include("energy_short.jl")
-    # include("energy_mid.jl")
-    # include("interpolate.jl")
-    # include("energy_long.jl")
-    # include("Chebyshev.jl")
-    # include("Taylor.jl")
-    # include("scaling.jl")
+    include("U_series.jl")
+    include("energy.jl")
+    include("energy_short.jl")
+    include("energy_mid.jl")
+    include("interpolate.jl")
+    include("energy_long.jl")
+    include("Chebyshev.jl")
+    include("Taylor.jl")
+    include("scaling.jl")
     include("zeroth_order.jl")
 end