perilap3d: triply periodic electrostatic kernels for general unit cell

version 0.7 (11/17/19)

Author: Alex H Barnett

This python/numba library computes the potential and fields at a set of targets inside a given general unit cell, due to a triply-periodized set of charges (which must sum to zero) and dipoles in the unit cell, to a requested accuracy tolerance. The applications include molecular dynamics and density functional theory in crystals. The parallelepiped unit cell is general (described by three lattice vectors), although currently it may not have high aspect ratio. Potentials and fields are also available at the sources themselves, where the self-interaction (j=i) is excluded from the sum. Instead of the usual idea of handling summation over the infinite lattice, the problem is solved as a PDE with periodic (hence nonlocal) boundary conditions. In particular, it writes the potential as a direct sum over only the given charges/dipoles and their near images, plus a smooth solution to the Laplace equation with inhomogeneous periodic BCs. The latter BVP solve is done to high accuracy via an auxiliary proxy point (method of fundamental solutions or MFS) representation. The whole scheme is compatible with the fast multipole method (FMM), although currently only plain direct summation is used.

For N sources and targets, where the nonperiodic direct evaluation cost is N², the cost of the periodic evaluation is about c³N² + O(p²N), where c is a small constant around 2, and p scales linearly with the number of digits required. Typically p=8 to 16, and the second term is around 10³N. For N=1000, the periodic evaluation is 10x slower than the non-periodic,at 3 digits of accuracy, and 20x slower at 9 digits. (If the FMM were used, both of the N² in the above would be replaced by O(N), and there would be ways to replace c by close to 1.) There is also a precomputation phase with cost O(p⁶), which need be done only once for a given unit cell shape, independent of the source or target locations.

Dependencies

This code has been tested in python3. It requires numpy, scipy, numba, and, if plotting is needed, matplotlib.

We recommend the Intel Distribution for Python, in which all of our tests were done.

Testing and usage

Run test.py for a complete test of the library. It will produce output similar to the file test.out, which has timings for an i7 laptop.

Here is a simple example (see demo.py):

from numpy import array,random,mean
import perilap3d as l3p

L = array([[1,0,0],[-.3,1.1,0],[.2,-.4,.9]])    # each row is a lattice vec 
p = l3p.lap3d3p(L)    # make a periodizing object for this lattice
p.precomp(tol=1e-6)   # do expensive precomputation (0.6 sec)

ns = 300                              # how many sources
y = (random.rand(ns,3)-1/2).dot(L)    # randomly in unit cell
d = random.randn(ns,3)                # dipole strength vectors
c = random.randn(ns); c -= mean(c)    # compatible charge strengths
pot,grad = p.eval(y,d,None,c)         # grad contains negatives of E field

The first call to eval after importing will require a few seconds to jit-compile the numba code. With this done, the above eval call takes 11 ms on an i7.

For another simple example, see madelung.py which computes the Madelung constant 1.74756459463318... for the cubic NaCl lattice to 15-digit accuracy in 3 ms (after a lattice-dependent setup time of 0.2 s).

See the test codes in perilap3d.py for more usage examples.

To Do

• add quadrupole sources and Hessian outputs

• spherical harmonics (or better?) for aux rep for rapid empty BVP solve

• fix fat Q case with QR solve? (only needed if m<P)

• pass-fail accuracy test, wider range of unit cells?

• way to save the Q factors for later use for that unit cell, for >8 digit acc

References

• For a reference on a similar scheme in 2D due to the author, combining with the FMM, see:

A unified integral equation scheme for doubly-periodic Laplace and Stokes boundary value problems in two dimensions, Alex H. Barnett, Gary Marple, Shravan Veerapaneni, Lin Zhao, Comm. Pure Appl. Math. 71(11), 2334-2380 (2018). http://arxiv.org/abs/1611.08038

• For a recent use of this idea in 3D cubic unit cells combining with the FMM (and showing that the extra near-neighbor cost can essentially be removed), see:

Flexibly imposing periodicity in kernel independent FMM: A Multipole-To-Local operator approach, Wen Yan and Michael Shelley, J. Comput. Phys. 355, 214-232 (2018). http://arxiv.org/abs/1705.02043

• For an accurate value of the Madelung constant and a very advanced way to compute it, see:

Ten problems in experimental mathematics, David H. Bailey, Jonathan M. Borwein, Vishaal Kapoor, and Eric W. Weisstein, Amer. Math. Monthly 113(6), 481-509 (2006).

Changelog

version 0.6, 9/14/18, dipoles only
version 0.7 11/17/19, charges added & Madelung example