zeroth with cheb

Project.toml

@@ -1,6 +1,6 @@
 name = "FastSpecSoG"
 uuid = "ebb58456-9453-4d17-87e2-7100c1ea036c"
-authors = ["Xuanzhao Gao <[email protected]> and contributors"]
+authors = ["Xuanzhao Gao <[email protected]>, Jiuyang Liang and Qi Zhou"]
 version = "1.0.0-DEV"
 
 [deps]
@@ -11,10 +11,8 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastChebInterp = "cf66c380-9a80-432c-aff8-4f9c79c0bdde"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
-OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-SumOfExpVPMR = "2c69873a-e0bb-44e1-90b4-d15ac3b7e936"
 
 [compat]
 julia = "1"

src/FFCT/Chebyshev.jl

@@ -1,3 +1,22 @@
+function USeries_direct_0(z::T, L_z::T, uspara::USeriesPara{T}, M_mid::Int) where{T}
+    t = zero(T)
+
+    for l in M_mid + 1:length(uspara.sw)
+        sl, wl = uspara.sw[l]
+        N = proper_N(1e-16, (L_z / sl)^2)
+
+        if N == 0
+            t += π * wl * sl^2 * exp(-z^2 / sl^2)
+        else
+            for k in 1:N
+                t += π * wl * sl^2 * special_powern(- z^2 / sl^2, k)
+            end
+        end
+    end
+
+    return t
+end
+
 function USeries_direct(z::T, k_xi::T, k_yj::T, uspara::USeriesPara{T}, M_mid::Int) where{T}
     t = zero(T)
     k2 = k_xi^2 + k_yj^2
@@ -10,6 +29,14 @@ function USeries_direct(z::T, k_xi::T, k_yj::T, uspara::USeriesPara{T}, M_mid::I
     return t
 end
 
+function ChebUSeries_0(L_z::T, uspara::USeriesPara{T}, M::Int, Q::Int) where{T}
+
+    f = z -> USeries_direct_0(z, L_z, uspara, M)
+    x = chebpoints(Q, zero(T), L_z)
+
+    return chebinterp(f.(x), zero(T), L_z)
+end
+
 function ChebUSeries_k(k_xi::T, k_yj::T, L_z::T, uspara::USeriesPara{T}, M::Int, Q::Int) where{T}
 
     f = z -> USeries_direct(z, k_xi, k_yj, uspara, M)

src/FFCT/zeroth_order.jl

@@ -12,30 +12,19 @@ function proper_N(accuracy::T, η::T) where{T}
     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}
+function FGT1d(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L_z::T, chebuseries::ChebPoly{1, T, T}, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}) where{T}
     n_atoms = length(qs)
     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]
-            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
+            grids[j] += q * chebuseries(abs(z - r_z[j]))
         end
     end
 
@@ -58,32 +47,12 @@ function FGT1d(qs::Vector{T}, poses::Vector{NTuple{3, T}}, s::T, L_z::T, r_z::Ve
     return E_k0
 end
 
-function FGT1d_naive(qs::Vector{T}, poses::Vector{NTuple{3, T}}, s::T) where{T}
-    n_atoms = length(qs)
-    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( - big(zi - zj)^2 / s^2)
-    end
-
-    return Float64(E_k0)
-end
 
-
-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}
+function zeroth_order(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, chebuseries::ChebPoly{1, T, T}, 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, r_z, chebcoefs, grids)
-    end
+    E_k0 = FGT1d(qs, poses, L_z, chebuseries, r_z, chebcoefs, grids)
 
-    return π * E_k0 / (2 * L_x * L_y)
+    return E_k0 / (2 * L_x * L_y)
 end

src/FSSoGInteraction.jl

@@ -42,7 +42,8 @@ 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, Q_0::Int; 
+    N_grid_long::NTuple{3, Int}, Q_long::Int, 
+    Rz_0::Int, Q_0::Int; 
     ϵ::T = one(T)) where{T}
 
     uspara = USeriesPara(b, σ, ω, M)
@@ -61,30 +62,33 @@ function FSSoGInteraction(
     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)
 
-    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Q_0)
+    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Rz_0)
+    chebuseries = ChebUSeries_0(L[3], uspara, M_mid, 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)
+    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, chebuseries)
 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, Q_0::Int; 
+    N_grid_long::NTuple{3, Int}, Q_long::Int, 
+    Rz_0::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, Q_0, ϵ = ϵ)
+    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
 
 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, Q_0::Int; 
+    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int,
+    Rz_0::Int, Q_0::Int; 
     ϵ::T = one(T)) where{T}
 
     uspara = USeriesPara(b, σ, ω, M)
@@ -98,20 +102,22 @@ 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)
 
-    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Q_0)
+    r_z0, grids0, chebcoefs0 = zero_paras_gen(L[3], Rz_0)
+    chebuseries = ChebUSeries_0(L[3], uspara, 0, 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, Q_0, r_z0, chebcoefs0, grids0)
+    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, chebuseries)
 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, Q_0::Int; 
+    N_real::NTuple{2, Int}, R_z::Int, w::NTuple{2, Int}, β::NTuple{2, T}, cheb_order::Int, Taylor_Q::Int,
+    Rz_0::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, Q_0, ϵ = ϵ)
+    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

src/FastSpecSoG.jl

@@ -1,11 +1,11 @@
 module FastSpecSoG
 
-using ExTinyMD, LinearAlgebra, SpecialFunctions, ChebParticleMesh, SumOfExpVPMR, LoopVectorization, OMEinsum, FastChebInterp, Polynomials, FFTW, DoubleFloats
+using ExTinyMD, LinearAlgebra, SpecialFunctions, ChebParticleMesh, FastChebInterp, Polynomials, FFTW, DoubleFloats
 using Base.Threads
 
 import FastChebInterp: ChebPoly
 
-export USeriesPara, U_series, BSA, ChebUSeries, proper_M
+export USeriesPara, U_series, BSA, ChebUSeries, ChebUSeries_0, proper_M
 export FSSoG_naive, FSSoGInteraction, FSSoGThinInteraction
 export short_energy_naive, long_energy_naive, energy_naive
 

src/energy/energy_long.jl

@@ -1,11 +1,8 @@
 function energy_long_0(
-    qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, M_mid::Int, 
-    uspara::USeriesPara{T}, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}
+    qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTuple{3, T}, chebuseries::ChebPoly{1, T, T}, r_z::Vector{T}, chebcoefs::Vector{T}, grids::Vector{T}
     ) where{T}
 
-    @assert M_mid ≤ length(uspara.sw)
-
-    E_0 = zeroth_order(qs, poses, L, uspara, M_mid, r_z, chebcoefs, grids)
+    E_0 = zeroth_order(qs, poses, L, chebuseries, r_z, chebcoefs, grids)
 
     return E_0
 end
@@ -14,7 +11,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.uspara, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
+    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, interaction.chebuseries, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
 
     return (E_k + E_0) / (4π * interaction.ϵ)
 end
@@ -23,7 +20,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.uspara, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
+    E_0 = energy_long_0(interaction.charge, interaction.position, interaction.L, interaction.chebuseries, interaction.r_z0, interaction.chebcoefs0, interaction.grids0)
     return (E_k + E_0) / (4π * interaction.ϵ)
 end
 

src/types.jl

@@ -57,6 +57,7 @@ mutable struct FSSoGInteraction{T} <: ExTinyMD.AbstractInteraction
     r_z0::Vector{T}
     chebcoefs0::Vector{T}
     grids0::Vector{T}
+    chebuseries::ChebPoly{1, T, T}
 end
 
 mutable struct FSSoGThinInteraction{T} <: ExTinyMD.AbstractInteraction
@@ -93,4 +94,5 @@ mutable struct FSSoGThinInteraction{T} <: ExTinyMD.AbstractInteraction
     r_z0::Vector{T}
     chebcoefs0::Vector{T}
     grids0::Vector{T}
+    chebuseries::ChebPoly{1, T, T}
 end

test/energy.jl

@@ -35,10 +35,11 @@
 
     N_grid = (32, 32, 32)
     Q = 24
+    R_z0 = 32
     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, 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)
@@ -81,10 +82,11 @@ end
     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, Q_0; 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, 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)

test/energy_long.jl

@@ -14,11 +14,12 @@
     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)
+    chebuseies = ChebUSeries_0(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)
     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, uspara, r_z0, chebcoefs0, grids0)
+    E_long_0 = energy_long_0(qs, poses, L, chebuseies, r_z0, chebcoefs0, grids0)
 
     @info "running the direct summation for the long range part of the energy"
     # using the direct summation

test/zeroth_order.jl

@@ -25,31 +25,6 @@ function long_energy_us_0_big(qs::Vector{T}, poses::Vector{NTuple{3, T}}, L::NTu
     return E0
 end
 
-@testset "FGT1d by mixed taylor chebshev" begin
-    @info "testing the FGT1d by mixed taylor chebshev"
-    n_atoms = 32
-    L = 20.0
-    Q_0 = 32
-
-    r_z, grids, chebcoefs = zero_paras_gen(L, Q_0)
-
-    qs = [(-1.0)^i for i in 1:n_atoms]
-    poses = [(rand() * L, rand() * L, rand() * L) for i in 1:n_atoms]
-
-    for preset in 1:6
-        uspara = USeriesPara(preset)
-
-        for l in 5:length(uspara.sw)
-            s, w = uspara.sw[l]
-            E_FGT = FGT1d(qs, poses, s, L, r_z, chebcoefs, grids)
-            E_naive = FGT1d_naive(qs, poses, s)
-            N = FastSpecSoG.proper_N(1e-16, (L / s)^2)
-            # @show l, N, E_FGT, E_naive, (E_FGT - E_naive) / E_naive
-            @test E_FGT ≈ E_naive
-        end
-    end
-end
-
 @testset "zeroth_order" begin
     @info "testing the zeroth_order energy"
     n_atoms = 32
@@ -64,7 +39,8 @@ end
     for preset in 1:6
         uspara = USeriesPara(preset)
         M_mid = 10
-        E_FGT = zeroth_order(qs, poses, (L, L, L), uspara, M_mid, r_z, chebcoefs, grids)
+        chebuseries = ChebUSeries_0(L, uspara, M_mid, Q_0)
+        E_FGT = zeroth_order(qs, poses, (L, L, L), chebuseries, r_z, chebcoefs, grids)
         E_naive = long_energy_us_0(qs, poses, (L, L, L), uspara, M_mid + 1, length(uspara.sw))
         E_naive_big = long_energy_us_0_big(qs, poses, (L, L, L), uspara, M_mid + 1, length(uspara.sw))
         # @show preset, E_FGT, E_naive, (E_FGT - E_naive) / E_naive