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 different treatments 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 incompatible, calculation of the non-Coulombic pseudopotential 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.
Z alpha, or g=0 non Coulomb term
The average value of a pseudopotential within its core radius may, and almost certainly does, differ from that of the equivalent Coulomb potential. The difference is known as the g=0 non-Coulomb term of the potential.
This difference enters the total energy of the system, in a term which is simply the sum of the g=0 non-Coulomb terms for all the potentials present, multiplied by the total charge and divided by the volume. This term is known as the "Z alpha" term from equations 19 and 25 of Ihm et al. (1979), also found in equation 2.14 of Payne, Teter et al.
The total charge was historically taken to be the total electronic charge, but, more recently, Bruneval et al. (2014), have shown that using the total ionic charge is more justifiable. Of course, for neutral systems there is no difference, but for charged systems there is a difference, one which scales as q/V (where q is the net charge).
So if one wishes to correct for terms which scale as q/V, one needs to be clear about which convention for the Z alpha term one is using. One resolution is to assign to each pseudopotential a quadrupole moment of 6ε0 times its g=0 non-Coulomb term. If these moments are added to the total quadrupole moment of the system, then the Makov Payne correction will be correct in the case of using the total ionic charge in the Z alpha term. Without it, the Makov Payne correction is correct if using the total electronic charge in the Z alpha term.
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.)
- Momentum-space formalism for the total energy of solids, J Ihm, A Zunger and ML Cohen, J. Phys. C: Solid State Phys. 21 4409 (1979)
- Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients MC Payne, MP Teter, et al., Rev. Mod. Phys. 64 1045 (1992)
- Periodic boundary conditions in ab initio calculations, G Makov and MC Payne, Phys Rev B 51 4014 (1995)
- Consistent treatment of charged systems within periodic boundary conditions: The projector augmented-wave and pseudopotential methods revisited F Bruneval et al., Phys. Rev. B 89 045116 (2014)
- Charged surfaces and slabs in periodic boundary conditions, MJ Rutter, Electron. Struct. 3 015002 (2021)
- Self-Consistent Potential Correction for Charged Periodic Systems, M Chagas da Silva et al., Phys Rev Lett 126 076401 (2021)