From 98cc170597b0c55a04a25575ce6f134dadb637de Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sun, 1 Jan 2023 13:57:08 -0700 Subject: [PATCH] Cleanup --- docs/math/ewald_review.lyx | 13 ++++-- docs/src/index.md | 8 ++-- src/Ewalder.jl | 86 ++++++++++++++++++++------------------ test/runtests.jl | 2 +- 4 files changed, 59 insertions(+), 50 deletions(-) diff --git a/docs/math/ewald_review.lyx b/docs/math/ewald_review.lyx index a62b7bd..5100ef5 100644 --- a/docs/math/ewald_review.lyx +++ b/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 diff --git a/docs/src/index.md b/docs/src/index.md index 56e4198..8d2151c 100644 --- a/docs/src/index.md +++ b/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 diff --git a/src/Ewalder.jl b/src/Ewalder.jl index e132a61..2427fad 100644 --- a/src/Ewalder.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index ac20ead..94c77d9 100644 --- a/test/runtests.jl +++ b/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