Cleanup

docs/math/ewald_review.lyx

@@ -1719,7 +1719,7 @@ Recall that the Gaussian charge cloud self interactions
 \begin_layout Standard
 \begin_inset Formula 
 \begin{align*}
-E_{\mathrm{self}} & =\frac{1}{4\pi\epsilon_{0}}\frac{1}{2}\sum_{i,j}\delta_{ij}q_{i}q_{j}\frac{1}{r_{ij\mathbf{0}}}\mathrm{erf}\left(\frac{r_{ij\mathbf{0}}}{\sqrt{2}\,\sigma}\right),
+E_{\mathrm{self}} & =\frac{1}{4\pi\epsilon_{0}}\frac{1}{2}\sum_{i,j}\delta_{ij}q_{i}q_{j}\frac{1}{r_{ij}}\mathrm{erf}\left(\frac{r_{ij}}{\sqrt{2}\,\sigma}\right),
 \end{align*}
 
 \end_inset
@@ -1835,7 +1835,7 @@ status open
 \begin_layout Plain Layout
 \begin_inset Formula 
 \begin{equation}
-E_{\mathrm{self}}=\frac{1}{4\pi\epsilon_{0}}\left(\frac{1}{\sqrt{2\pi}\,\sigma}\sum_{i}q_{i}^{2}+\frac{1}{3}\frac{1}{\sqrt{2\pi}\,\sigma^{3}}\sum_{i}p_{i}^{2}\right).
+E_{\mathrm{self}}=\frac{1}{4\pi\epsilon_{0}}\left(\frac{1}{\sqrt{2\pi}\,\sigma}\sum_{i}q_{i}^{2}+\frac{1}{3\sqrt{2\pi}\,\sigma^{3}}\sum_{i}p_{i}^{2}\right).
 \end{equation}
 
 \end_inset
@@ -1903,7 +1903,7 @@ Consider, for simplicity, an ordered pair of dipoles
  The energy for this periodic, pairwise interaction is
 \begin_inset Formula 
 \begin{equation}
-E_{ij}=\frac{1}{4\pi\epsilon_{0}}\frac{1}{2}\mathbf{p}_{i}\cdot\left[\frac{4\pi}{V}\sum_{\mathbf{k}\neq\mathbf{0}}\frac{e^{-\sigma^{2}k^{2}/2}}{k^{2}}(\mathbf{k}\otimes\mathbf{k})e^{-i\mathbf{k}\cdot\mathbf{r}_{ij}}+\sideset{}{^{'}}\sum_{\mathbf{n}}\overleftrightarrow{\mathcal{E}}_{\mathrm{dd}}(\mathbf{r}_{ij\mathbf{n}})+\frac{2\delta_{ij}}{3\sigma^{3}}\sqrt{\frac{1}{2\pi}}\,I\right]\cdot\mathbf{p}_{j},
+E_{ij}=\frac{1}{4\pi\epsilon_{0}}\mathbf{p}_{i}\cdot\left[\frac{1}{2}\sideset{}{^{'}}\sum_{\mathbf{n}}\overleftrightarrow{\mathcal{E}}_{\mathrm{dd}}(\mathbf{r}_{ij\mathbf{n}})+\frac{4\pi}{2V}\sum_{\mathbf{k}\neq\mathbf{0}}\frac{e^{-\sigma^{2}k^{2}/2}}{k^{2}}(\mathbf{k}\otimes\mathbf{k})e^{-i\mathbf{k}\cdot\mathbf{r}_{ij}}+\frac{\delta_{ij}}{3\sqrt{2\pi}\,\sigma^{3}}\sqrt{\frac{1}{2\pi}}\,I\right]\cdot\mathbf{p}_{j},
 \end{equation}
 
 \end_inset
@@ -1927,7 +1927,12 @@ noprefix "false"
 \end_inset
 
 .
- A factor of two is required to account for contributions from both 
+ When 
+\begin_inset Formula $i\neq j$
+\end_inset
+
+, a factor of two should be applied to account for contributions from both
+
 \begin_inset Formula $E_{ij}$
 \end_inset
 

docs/src/index.md

@@ -58,7 +58,7 @@ Ewald summation works by decomposing the Coulomb energy into two parts: (1) A
 real-space part, evaluated as a sum over point charges at short range, and (2) a
 Fourier space part, which captures the long-range interactions through a sum
 over $\mathbf k$-vectors, which label frequencies in the decomposition of total
-charge density. The `Ewalder.jl` implementation uses a straightforward algorithm
+charge density. The _Ewalder_ implementation uses a straightforward algorithm
 that scales as $N^{3/2}$ for increasing system sizes $N$ (assuming efficient
 calculation of the neighbor list). For very large scale simulations, it would be
 preferable to use instead a method that scales near-linearly with $N$, such as
@@ -68,7 +68,7 @@ preferable to use instead a method that scales near-linearly with $N$, such as
 
 The full Ewald method is summarized below. See our [mathematical
 writeup](https://raw.githubusercontent.com/SunnySuite/Ewalder.jl/main/docs/math/ewald_review.pdf)
-(PDF format) for the full derivation. 
+(PDF format) for the derivation. 
 
 The Ewald energy for a system of charges decomposes into three parts,
 
@@ -102,7 +102,7 @@ E_{\mathrm{self}}=\frac{1}{4\pi\epsilon_{0}}\frac{1}{\sqrt{2\pi}\,\sigma}\sum_{i
 ```
 which must be removed from the final energy output $E$.
 
-Both real- and Fourier- space sums are rapidly convergent. `Ewalder.jl` truncates these sums using the inequalities,
+Both real- and Fourier- space sums are rapidly convergent. _Ewalder_ truncates these sums using the inequalities,
 ```math
 r_{ij\mathbf{n}} \leq c_0 \sqrt{2} \sigma \\
 |\mathbf{k}_{\mathbf{m}}| \leq c_0 \sqrt{2} / \sigma
@@ -117,7 +117,7 @@ Apart from the small truncation error, the energy $E$ is mathematically independ
 
 where $c_1$ is a dimensionless parameter (defaulting to $c_1 = 2$) and $L = V^{1/3}$. This choice achieves, in principle, a total computational cost that scales as $O(N^{3/2})$, provided that the ion neighbor lists are calculated efficiently. For maximum efficiency, the user should provide a neighbor list. If one is not provided, this package will perform inefficient, brute-force calculation of the neighbor list.
 
-In addition to point charges, `Ewalder.jl` also allows specification of point dipoles, each of which generates an electrostatic potential that decays like $1/r^2$, and can similarly be handled using Ewald summation. The final formulas can be found in our [PDF note](https://raw.githubusercontent.com/SunnySuite/Ewalder.jl/main/docs/math/ewald_review.pdf).
+In addition to point charges, _Ewalder_ also allows specification of point dipoles, each of which generates an electrostatic potential that decays like $1/r^2$, and can similarly be handled using Ewald summation. The final formulas can be found in our [PDF note](https://raw.githubusercontent.com/SunnySuite/Ewalder.jl/main/docs/math/ewald_review.pdf).
 
 ## API
 

src/Ewalder.jl

@@ -1,7 +1,7 @@
 module Ewalder
 
 import StaticArrays: SVector, SMatrix, SA
-import LinearAlgebra: normalize, ⋅, ×, inv
+import LinearAlgebra: normalize, inv, ⋅, ×
 import SpecialFunctions: erfc
 
 const Vec3 = SVector{3,Float64}
@@ -66,7 +66,17 @@ end
 real_space_cutoff(sys::System)    = √2 * sys.c0 * sigma(sys)
 fourier_space_cutoff(sys::System) = √2 * sys.c0 / sigma(sys)
 
+# Number of cells in each direction to find all cells within a center-to-center
+# distance of rmax.
+function required_cell_displacement(latvecs, rmax)
+    recipvecs = eachrow(inv(reduce(hcat, latvecs)))
+    return map(latvecs, recipvecs) do a, b
+        round(Int, rmax / (a⋅normalize(b)) + 1e-6)
+    end
+end
 
+# Ensure that every position is within the parallelpiped unit cell defined by
+# sys.latvecs
 function wrap_positions!(sys::System; warn=false)
     # Lattice vectors as a 3x3 matrix
     A = reduce(hcat, sys.latvecs)
@@ -88,32 +98,26 @@ function wrap_positions!(sys::System; warn=false)
 end
 
 # Perform O(N^2) calculation of neighbor list.
-function get_neighbors(sys::System; rmax=0.0)
+function get_neighbors(sys::System)
     # Make sure positions are valid
     wrap_positions!(sys; warn=true)
 
-    # Real space cutoff
-    if iszero(rmax)
-        rmax = real_space_cutoff(sys)
-    end
-
-    # Maximum latvec displacements, in each direction
-    as = sys.latvecs
-    bs = recipvecs(sys)
-    nmax = map(as, bs) do a, b
-        round(Int, rmax / (a⋅normalize(b)) + 1e-6) + 1
-    end
-
-    N = length(sys.pos)
-
-    ret = Neighbor[]
-    for i = 1:N, j = 1:N
-        for n1 = -nmax[1]:nmax[1], n2 = -nmax[2]:nmax[2], n3 = -nmax[3]:nmax[3]
-            n = (n1, n2, n3)
-            neigh = Neighbor(; i, j, n)
-            if distance2(sys, neigh) <= rmax*rmax
-                if i != j || !iszero(Vec3(n))
-                    push!(ret, neigh)
+    # Target is to find all neighbors within distance rmax
+    rmax = real_space_cutoff(sys)
+
+    # Maximum latvec displacements in each direction. Shift by 1 to account for
+    # atoms that are at the edges of the cell.
+    nmax = required_cell_displacement(sys.latvecs, rmax) .+ 1
+
+    ret = Vector{Neighbor}[]
+    for i = eachindex(sys.pos)
+        push!(ret, Neighbor[])
+        for j = eachindex(sys.pos)
+            for n1 = -nmax[1]:nmax[1], n2 = -nmax[2]:nmax[2], n3 = -nmax[3]:nmax[3]
+                n = (n1, n2, n3)
+                neigh = Neighbor(; i, j, n)
+                if 0 < distance2(sys, neigh) <= rmax*rmax
+                    push!(ret[i], neigh)
                 end
             end
         end
@@ -127,9 +131,10 @@ end
 
 Calculate the Ewald energy in units of $1/4\pi\epsilon_0$ for a periodic system
 of `charges` and `dipoles`. If either charge or dipole parameters are omitted,
-they will be assumed zero.
+they will be assumed zero. A neighbors list can also be optionally provided; see
+the source code for details.
 """
-function energy(sys::System; neighbors=Neighbor[], charges=Float64[], dipoles=Vec3[])
+function energy(sys::System; neighbors=Vector{Neighbor}[], charges=Float64[], dipoles=Vec3[])
     N = length(sys.pos)
 
     # Find neighbors if not provided
@@ -146,16 +151,20 @@ function energy(sys::System; neighbors=Neighbor[], charges=Float64[], dipoles=Ve
     @assert abs(Q) < 1e-12
     charges .-= Q / N
 
+    # Precalculate constants
+    σ = sigma(sys)
+    σ² = σ*σ
+    σ³ = σ^3
+
     # Energy accumulator
     E = 0.0
 
     #####################################################
     ## Real space sum
-    σ = sigma(sys)
-    σ² = σ*σ
 
-    for neigh = neighbors
-        (; i, j) = neigh
+    for i = 1:N, neigh = neighbors[i]
+        @assert i == neigh.i
+        j = neigh.j
 
         qᵢ, pᵢ = (charges[i], dipoles[i])
         qⱼ, pⱼ = (charges[j], dipoles[j])
@@ -182,19 +191,15 @@ function energy(sys::System; neighbors=Neighbor[], charges=Float64[], dipoles=Ve
 
     #####################################################
     ## Fourier space sum
-    V = volume(sys)
-    as = sys.latvecs
-    bs = recipvecs(sys)
 
-    # Maximum recipvec displacements, in each direction
+    V = volume(sys)
+    recip = recipvecs(sys)
     kmax = fourier_space_cutoff(sys)
-    mmax = map(as, bs) do a, b
-        round(Int, kmax / (normalize(a)⋅b) + 1e-6) + 1
-    end
+    mmax = required_cell_displacement(recip, kmax)
 
     # Loop over grid of k points
     for m1 = -mmax[1]:mmax[1], m2 = -mmax[2]:mmax[2], m3 = -mmax[3]:mmax[3]
-        k = Vec3(m1, m2, m3)' * bs
+        k = Vec3(m1, m2, m3)' * recip
         k² = k⋅k
         if 0 < k² <= kmax*kmax
             ρk = 0.0im
@@ -205,12 +210,11 @@ function energy(sys::System; neighbors=Neighbor[], charges=Float64[], dipoles=Ve
         end
     end
 
-
     #####################################################
     ## Remove self energies
-    σ³ = σ^3
+
     for (q, p) = zip(charges, dipoles)
-        E += - q*q/(√(2π)*σ) - (p⋅p)/(√(2π)*(3σ³))
+        E += - (q*q)/(√(2π)*σ) - (p⋅p)/(3√(2π)*σ³)
     end
 
     return E

test/runtests.jl

@@ -7,7 +7,7 @@ using TestItemRunner
     using LinearAlgebra
 
     # Madelung constants of ionic crystals:
-    # - Y. Sakamoto, The J. of Chem. Phys., vol. 28, (1958), p. 164
+    # * Y. Sakamoto, J. Chem. Phys., vol. 28, (1958), p. 164
     E_NaCl = -1.7475645946331822
     E_CsCl = -1.76267477307099