EwaldSummations

Build Status Build Status Coverage

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.

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.

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.

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

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)

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)

# consider 100 * 100 periodic images, and 20 layers of image charges
sys_q2d = SysQ2D((0.9, 0.9), (L, L, L), 50, 10, ϵ = 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)

How to Contribute

If you find any bug or have any suggestion, please open an issue.