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 p2/(ε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 p2/(ε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 with dipole

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: