EwaldSummations
EwaldSummations.jl
is an implementation of the basic Ewald summation algorithm, including Ewald3D and Ewald2D. This package also include the naive approach: sum directly over periodic directions.
Specially, for dielectric confined quasi2D systems, this package provide a set of tools based on the image charge method, combining with the Ewald2D or Ewald3DELC method.
This package is designed to provide standard/correct results of electrostatic interaction in charged particle systems, and can be used as a tool box when developing new algorithms.
For some of the implemented methods, the package support parallel computation by Base.Threads
, start that by julia t 8 your_script.jl
.
Getting Started
Add this package and its deps ExTinyMD.jl
in julia by typing ]
in Julia REPL and then
pkg> add ExTinyMD, EwaldSummations
to install the package.
Numerical Methods
In this package, we provide various methods to calculate the electrostatic interaction energy of charged particle systems. For system without dielectric mismatch, the following methods are available:
 Direction summation: sum directly over periodic directions (3D and quasi2D)
 Ewald3D: the standard Ewald summation for 3D systems
 Ewald2D: the standard Ewald summation for 2D systems
Then for quasi2D systems confined by dielectric mismatch, we combine the image charge method and the methods above, given by
 Image charge reflection: generating the image charge series according to the dielectric mismatch
 ICM + direct summation: combining the image charge method and the direct summation
 ICM + Ewald2D: combining the image charge method and the Ewald2D method
 ICM + Ewald3D + ELC: extending the quasi2D system as triply periodic, and calculate the energy via Ewald3D and electrostatic layer correction (ELC)
For detailed numerical formula, please refer to the reference^{1}.
Examples of Usage
Here are some examples of usage, for more details please see the test files.
Triply periodic systems
Here is an example for calculating energy of 3D systems, via Ewald3D or direction summation.
# Init the system
n_atoms = 100
L = 100.0
boundary = Boundary((L, L, L), (1, 1, 1))
atoms = Vector{Atom{Float64}}()
for i in 1:n_atoms÷2
push!(atoms, Atom(type = 1, mass = 1.0, charge = 1.0))
end
for i in n_atoms÷2 + 1 : n_atoms
push!(atoms, Atom(type = 2, mass = 1.0, charge =  1.0))
end
# randomly generate position of particles
info = SimulationInfo(n_atoms, atoms, (0.0, L, 0.0, L, 0.0, L), boundary; min_r = 1.0, temp = 1.0)
# take alpha = 0.2 and s = 3.0
Ewald3D_interaction = Ewald3DInteraction(n_atoms, 3.0, 0.2, (L, L, L), ϵ = 1.0, ϵ_inf = 1.0)
# setup cell list for short range interaction
neighbor = CellList3D(info, Ewald3D_interaction.r_c, boundary, 1)
# calculation interaction energy
energy_ewald = energy(Ewald3D_interaction, neighbor, info, atoms)
# direction summation, take 60 * 60 * 60 periodic images into account
sys_3d = Sys3D((L, L, L), (30, 30, 30), ϵ = 1.0)
coords = [p_info.position for p_info in info.particle_info]
charge = [atoms[p_info.id].charge for p_info in info.particle_info]
energy_direct = Energy_3D(sys_3d, coords, charge)
Doubly periodic systems without dielectric mismatch
Here is an example for calculating energy of Q2D systems, via Ewald3D or direction summation.
n_atoms = 100
L = 100.0
boundary = ExTinyMD.Q2dBoundary(L, L, L)
atoms = Vector{Atom{Float64}}()
for i in 1:n_atoms÷2
push!(atoms, Atom(type = 1, mass = 1.0, charge = 1.0))
end
for i in n_atoms÷2 + 1 : n_atoms
push!(atoms, Atom(type = 2, mass = 1.0, charge =  1.0))
end
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, 3.0, 0.2, (L, L, L), ϵ = 1.0)
neighbor = CellList3D(info, Ewald2D_interaction.r_c, boundary, 1)
energy_ewald = energy(Ewald2D_interaction, neighbor, info, atoms)
# consider 400 * 400 periodic images, with no boundary mismatch
sys_q2d = SysQ2D((0.0, 0.0), (L, L, L), 200, 0, ϵ = 1.0)
coords = [p_info.position for p_info in info.particle_info]
charge = [atoms[p_info.id].charge for p_info in info.particle_info]
ref_pos, ref_charge = SysQ2DInit(sys_q2d, coords, charge)
energy_icm = Energy_Q2D(sys_q2d, coords, charge, ref_pos, ref_charge)
Doubly periodic systems with dielectric mismatch
Here is an example for calculating energy of Q2D systems confined by dielectric mismatch, via direction summation together with image charge method.
n_atoms = 100
L = 100.0
boundary = ExTinyMD.Q2dBoundary(L, L, L)
s = 6.0
alpha = 2.0
γ = (0.9, 0.9)
N_img = 10
N_pad = 20
N_real = 50
atoms = Vector{Atom{Float64}}()
for i in 1:n_atoms÷2
push!(atoms, Atom(type = 1, mass = 1.0, charge = 1.0))
end
for i in n_atoms÷2 + 1 : n_atoms
push!(atoms, Atom(type = 2, mass = 1.0, charge =  1.0))
end
info = SimulationInfo(n_atoms, atoms, (0.0, L, 0.0, L, 0.0, L), boundary; min_r = 1.0, temp = 1.0)
# direction summations, consider 100 * 100 periodic images, and 20 layers of image charges
sys_q2d = SysQ2D(γ, (L, L, L), N_real, N_img, ϵ = 1.0)
coords = [p_info.position for p_info in info.particle_info]
charge = [atoms[p_info.id].charge for p_info in info.particle_info]
ref_pos, ref_charge = SysQ2DInit(sys_q2d, coords, charge)
energy_icm = Energy_Q2D(sys_q2d, coords, charge, ref_pos, ref_charge)
# ICMEwald2D, consider 20 layers of image charges
ICMEwald2D_interaction = IcmEwald2DInteraction(n_atoms, s, alpha, γ, (L, L, L), N_img)
energy_icmewald2d = ICM_Ewald2D_energy(ICMEwald2D_interaction, coords, charge)
# ICMEwald3D, consider 20 layers of image charges and 40 layers of padding
ICMEwald3D_interaction = IcmEwald3DInteraction(n_atoms, s, alpha, γ, (L, L, L), N_pad, N_img)
energy_icmewald3d = ICMEwald3D_energy(ICMEwald3D_interaction, coords, charge)
How to Contribute
If you find any bug or have any suggestion, please open an issue.

J. Yuan, H. S. Antila, and E. Luijten, Particle–particle particle–mesh algorithm for electrolytes between charged dielectric interfaces, J. Chem. Phys., 154 (2021), p. 094115. ↩