Periodic boundary conditions and charges
It can readily be shown that the energy, per unit cell, of an infinite 3D tiling of charged cells is infinite. A bigger problem than that which arises with PBCs and dipole moments. So to get a meaningful, i.e. finite, answer, the cell must be neutral.
The simplest, and parameter-free, way of achieving this is to add a uniform neutralising charge density throughout the cell. It turns out that the usual way of calculating electrostatic energies in 3D periodic codes does this automatically. The g=0 Fourier components of the components (ionic and electronic) charge density are ignored and assumed to sum to zero, whereas the other Fourier components are treated correctly. If the g=0 components did not sum to zero, the effect is to add a uniform, neutralising, charge density.
This extra charge density varies with cell size, causes an energy term which varies with cell size, and causes extraneous electric fields which vary with cell size. It is best removed in a self-consistent fashion.
However, for some geometries there are known post hoc corrections to the energy, and those for systems whose true periodicity is 0D (isolated systems) and 2D (slabs and surfaces) are known to c2x version 2.40.
The 0D corrections for an isolated system modelled in a cubic (or other shaped) box are described by Makov and Payne, and include a term which depends on the charge, and one which depends on the quadrupole moment. The uncorrected energy converges not as 1/V, but merely as 1/L where L is the length of a lattice vector. Adding the charge-dependent term improves the convergence to 1/V, and the quadrupole term improves it to 1/L5. These corrections can be complicated by an incorrect treatment of the non-Coulombic g=0 term of the pseudopotential energy.
The figures below show the energy of a Li+ ion in a cubic box of the given size after correction by the charge, and both the charge and quadrupole, terms of the Makov-Payne correction. The first figure uses a common, but incorrect, calculation of the non-Coulombic psuedopotential energy, and the second a correct calculation.
$ c2x -q --null 10/Li.castep_bin (or 10/Li.check) Madelung energy correction: 2.042804 eV Original energy: -194.410917 eV Corrected energy: -192.368113 eV Quadrupole correction: -0.008019 eV Final corrected energy: -192.360094 eV
The c2x command line for the 10A calculation is given above, and it can be seen that the completely uncorrected energy is -194.4eV, and would be far off-scale on these graphs.
A corresponding correction for 2D systems, again including a charge and a quadrupole term, is given by Rutter. Applying the correction in c2x is very similar, save that the non-periodic axis (a, b or c) must be stated.
$ c2x -qc --null 10/graphene.castep_bin Charged slab correction: -57.649244 eV Original energy: -276.411603 eV Corrected energy: -334.060846 eV Total charge: 2.000000 Total dipole (eA): (-0.000000,-0.000000,0.000000) magnitude 0.000000 Quadrupole correction: 2.826319 eV Final corrected energy: -331.234527 eV
Here the c axis is specified to be the non-periodic one.
The correction is large because, in the non-periodic case of an isolated charged sheet, the field in the vacuum is a constant. The energy density of a constant field is ½ε0E2, so the energy is linear in the volume of the cell considered. This issue does not arise in the 0D case. In the 2D case, c2x corrects for it not by trying to extrapolate to infinite volume, and hence infinite energy, but to zero volume.
In most cases self-consistent corrections are much superior to any post hoc ones which c2x could apply. As one example, consider a sample which is conducting. Conductors cannot support a macroscopic field, for their macroscopic potential is constant. It follows that their macroscopic internal charge density must be zero, and that any net charge must reside on their surfaces.
If a constant charge density is added to a cell, in regions where this overlaps with a conductor charge will move from the surface of the conductor into its bulk. Thus this extra charge, lying where a net charge is prohibited, is itself neutralised. This can significantly change surface charges and fields, and do so in a way which is strongly volume dependent.
So one needs a self-consistent approach, such as those described by Rutter, or by Chagas da Silva at al., or one of the many other methods such as truncating the Coulomb potential, or reformulating the electrostatic energy calculation in a lower dimension, etc.
Negatively charged systems give rise to an extra complication. It is likely that some electrons will prefer to lie in a puddle in the vacuum rather than in the negatively-charged specimen. For the 2D geometry this is inevitable. For a 0D geometry, some ions are stable, such as Cl-, but most, including most multiply-charged ions, are not. If the basis states extend into the vacuum, as will be the case with plane wave codes, electrons are likely to escape.
(Much has been written on this topic. The two 2021 papers below should be considered not only for their own merits, but also for the useful lists of references to other approaches which they provide.)
- Periodic boundary conditions in ab initio calculations, G Makov and MC Payne, Phys Rev B 51 4014 (1995)
- Charged surfaces and slabs in periodic boundary conditions, MJ Rutter, Electron. Struct. (2021)
- Self-Consistent Potential Correction for Charged Periodic Systems, M Chagas da Silva et al., Phys Rev Lett 126 076401 (2021)