# Periodic boundary conditions and dipoles

Plane-wave DFT codes such as Castep generally work with 3D periodicity. When working with crystals, this is precisely what one wants. However, some systems lack 3D periodicity. Obvious examples are isolated molecules (0D), isolated nanowires or nanotubes (1D), and surfaces (2D). Even if one's system is fully periodic, there are still interesting questions, for the dipole moment of a periodic system is not well defined.

If one's simulation cell is large enough that the unwanted periodic images of one's system are sufficiently far away that they do not interact, then all is fine. This can require quite a large cell, which is expensive in both memory and CPU time.

For uncharged systems, the worst interaction, in that it decays
slowest with increasing cell size, is the dipole-dipole
interaction. This leads to an energy term of approximately
p^{2}/(εV), where p is the dipole moment
and V the cell volume. It is not only the energy which is
changed. Because the electrostatic potential is forced to be periodic,
and a dipole would usually introduce a step change in the potential,
the periodic boundary conditions effectively introduce a fictitious
electric field. This field will then cause the electrons to react,
and will produce forces on ions, and thus perturb the system from the
state it would occupy in the absence of such a field.

For isolated systems which can be padded with vacuum (i.e. not defects in a bulk crystal), there are two approaches to correcting for the interaction of the dipole moment with its fictitious images.

The first, and best, is to remove the fictitious field described above by placing a step in the potential in the vacuum region. This is referred to as a self-consistent correction, and should produce the best convergence. Castep can do this for surfaces (i.e. slabs) from version 19. One simply needs to give the direction perpendicular to the slab.

The second approach is to add on a post hoc correction to the
energy - the p^{2}/(εV) term, but taking care to get
the correct prefactor. This is not self-consistent, in that the
energy is being corrected, but the field is not. There are
arguments that the "correct" value for 1/ε is the value
averaged over the cell, which may contain a significant fraction of
material, not simply the vacuum value
1/ε_{0}. The average value is much harder
to calculate, so the vacuum value is generally used.

For any correction, it is wise to plot quantities such as energy and dipole moment for increasing vacuum, and to observe their convergence. The post hoc correction should improve the convergence of the energy, but will not improve the convergence of properties such as the system's dipole moment. The dipole moment will in general be enhanced by the presence of the images.

*Energy convergence for an
artificial system consisting of a slab of K atoms and a slab of Cl
atoms, separated by 2.7A, with the cell extent in the non-periodic
direction as shown, and a fixed 6A in the perpendicular
directions.*

The above graph shows that even a simple post hoc correction can significantly improve the convergence rate. It was calculated with:

$ c2x -Dc= -c --null KCl.check Ed: 0.000008 -0.000015 -1.000498 Total charge: -0.000000 Total dipole (eA): (0.000045,-0.000088,-1.909958) magnitude 1.909958 Calculated dipole energy correction (c axis): 0.458406 eV Corrected energy -451.005193 + 0.458406 = -450.546787 eV

The above being for the 20A cell length. Note that it is necessary
to tell c2x to read the charge density (`-c`) in order for it to be
able to calculate a dipole moment or correction. The
syntax `-Dc=` is shorthand for `-Dc=0.5,0.5,0.5`
and means calculate dipole (`-D`), include energy
correction for slab geometry with the c axis being the aperiodic one
(`c`), and with the point about which the dipole is
calculated being 0.5,0.5,0.5 (`=` with 0.5,0.5,0.5
assumed).

For an isolated molecule in a cubic box, the post hoc correction can be
calculated with `-Dm=` (see Makov & Payne and
Kantorovich below). From version 2.27 of c2x, this is
extended slightly to isolated molecules in tetragonal boxes with
the dipole moment parallel to the c axis (see Rutter below).

It should be stressed that the use of such corrections does not remove the need to check convergence with cell size. There will be other effects which these corrections do not address. For the post hoc corrections in c2x the most significant is usually the system's response to the extra electric field imposed, but even with a fully self-consistent correction there may still be terms due to quadrapole and higher moments, and due to wavefunction overlap.

Note too that the Rutter paper below describes a trick of using a tetragonal box with c/a=1.5 for calculations on molecules, arguing that this significantly reduces the extra electric field should a fully self-consistent correction be unavailable.

The theory of the corrections is not trivial, for it involves conditionally convergent sums and careful consideration of the boundary conditions at infinity. Exact results for the prefactor are known for particular geometries only. As starting points for further reading, one might consider:

*Periodic boundary conditions in ab initio calculations,*G Makov and MC Payne, Phys Rev B**51**4014 (1995)*The electrostatic potential in multipole lattices,*FW De Wette and BRA Nijboer, Physica**24**1105 (1958)*Adsorbate-substrate and adsorbate-adsorbate interactions of Na and K adlayers on Al(111),*J Neugebauer and M Scheffler, Phys Rev B**46**16067 (1992)*Elimination of the long-range dipole interaction in calculations with periodic boundary conditions,*LN Kantorovich, Phys Rev B**60**15476 (1999)*Molecules with dipoles in periodic boundary conditions in a tetragonal cell,*MJ Rutter J. Phys.: Condens. Matter**31**335901 (2019)