diff --git a/Project.toml b/Project.toml index f2124a7..38a8e06 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastSpecSoG" uuid = "ebb58456-9453-4d17-87e2-7100c1ea036c" -authors = ["Xuanzhao Gao and contributors"] +authors = ["Xuanzhao Gao , 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" diff --git a/src/FFCT/Chebyshev.jl b/src/FFCT/Chebyshev.jl index 7f9b655..713b127 100644 --- a/src/FFCT/Chebyshev.jl +++ b/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) diff --git a/src/FFCT/zeroth_order.jl b/src/FFCT/zeroth_order.jl index 18fa23a..be85781 100644 --- a/src/FFCT/zeroth_order.jl +++ b/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 \ No newline at end of file diff --git a/src/FSSoGInteraction.jl b/src/FSSoGInteraction.jl index 6ed04c3..6fd8cff 100644 --- a/src/FSSoGInteraction.jl +++ b/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 \ No newline at end of file diff --git a/src/FastSpecSoG.jl b/src/FastSpecSoG.jl index 568bb2e..180f88c 100644 --- a/src/FastSpecSoG.jl +++ b/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 diff --git a/src/energy/energy_long.jl b/src/energy/energy_long.jl index 91ac70c..287cefc 100644 --- a/src/energy/energy_long.jl +++ b/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 diff --git a/src/types.jl b/src/types.jl index 0ee904a..d547087 100644 --- a/src/types.jl +++ b/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 \ No newline at end of file diff --git a/test/energy.jl b/test/energy.jl index 5da1d65..8a57e27 100644 --- a/test/energy.jl +++ b/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) diff --git a/test/energy_long.jl b/test/energy_long.jl index e0e9742..1a6d726 100644 --- a/test/energy_long.jl +++ b/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 diff --git a/test/zeroth_order.jl b/test/zeroth_order.jl index c9d7dee..c3bed5d 100644 --- a/test/zeroth_order.jl +++ b/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