c2x and Potentials
From version 2.20, c2x can calculate electrostatic potentials. So herewith a few notes on how it does so.
To do so it needs to be able to read an electronic charge density and atomic charges and positions. These are found in a Castep .check file, and might also be found in a .cube file.
It should be remembered that the zero of the electrostatic potential for an infinitely periodic system is not well defined.
For the electronic contribution to the potential, c2x transforms the electron density to reciprocal space, divides by g2, and transforms back to real space. That is entirely conventional.
For the ionic contribution, c2x has two techniques. Considering the charge density of the ions to be delta functions is very unhelpful: it leads to poles in the potential. So conventionally the ions are considered to be Gaussian packets of charge. This c2x does, and one may vary the width of the Gaussian smearing. By default, the smearing in real space is exp(-μ2r2) with μ=4.2A-1. The smearing width can be varied.
Such Gaussian smearing of the density gives rise to a potential which varies like erf(μr)/r. This is a monotonic function, with a finite quadratic maximum at the origin. However, the maximum is sharper than the local part of most pseudopotentials, yet c2x does not read and cannot make use of pseudopotentials. So for atoms for which the atomic number and the charge differ, indicating that core electrons exist in the pseudopotential, c2x modifies the above ionic charge density.
The charge-neutral, localised, modification is to subtract a density which gives rise to a Gaussian potential, and to choose the Gaussian potential so as to flatten the peak of the total ionic potential. This generally gives a closer approximation to a genuine pseudopotential.
The above image, of potential in volts as a function of distance in Angstroms, shows the potentials arising from bare Gaussian densities with μ=4 and 6, and from a modified Gaussian density with μ=6. The true 1/r potential has a pole, and the unmodified Gaussian with μ=6 models this most closely. However, the modified Gaussian density with μ=6 has a similar peak value to the bare μ=4 density, usefully lower than the bare μ=6 density, but it follows the 1/r potential more closely at around ±0.3A than the bare μ=4 Gaussian density.
Smearing widths
The potential will be "wrong" inside the smearing region, just as it will be "wrong" inside a pseudopotential core region. A small smearing width will lead to large peak values, and introduce aliasing problems as the charge will not be well localised in reciprocal space. A large width may lead to "wrong" values in the area of interest. The potential should reach about 85% of its true value at r=1/μ, over 96% of its true value at r=1.5/μ, and 99.5% of its true value at r=2/μ.
One can vary μ on the command line by specifying, e.g., "-E=5", to set μ to 5. With -vv, c2x prints some data warning about aliasing:
Using smearing of 1/4.250 A and both bare and pseudo Coulomb potentials
Aliasing test in ion_recrho(): 0.00198131 0.00198131 0.00162353
-0.0062398 -0.0062398 -0.00532859
The first line gives an estimate of the fractional error introduced by aliasing along each of the three axes for the bare Coulomb potential, and the second line the same for the modified potential. The values above are probably at the high end of acceptable. Interpolating onto a finer grid can sharply reduce these figures.
Adjusting μ is likely to result in a shift of the potential everywhere, due to the convention that the average potential is zero.
One can turn off the use of the modified Coulomb potential, and simply use the bare Coulomb potential for all atoms. To do this, specify a negative μ, or, for the default value of μ, simply "-E=-".
Interpolation and Averaging
C2x can perform Fourier interpolation and averaging. Usually this is performed as a single step, but in the case of calculating potentials it is performed in two steps. The first step is to interpolate in any directions where a finer grid is requested, then the potential is calculated, and then any coarsening grid change is performed. As an example, if a calculation were performed on a 40x40x80 grid, and c2x asked to calculated the electrostatic potential and interpolate to a 1x1x100 grid (perhaps if wanting to see planar averages perpendicular to a slab), it would first interpolate to 40x40x100, then calculate the potential, then coarsen to 1x1x100.
Castep's .cst_esp file
Castep produces a potential in a .cst_esp file by default. This c2x can read, and, from version 2.20, it will by default rescale it from the original units of -Hartrees to +volts. One can specify "-R" to leave it unscaled (the behaviour before 2.20).
Castep's potential is quite different to that calculated by c2x. The biggest difference is that it includes the exchange-correlation potential. To test c2x's potential calculation against Castep's use of genuine pseudopotentials, it is necessary to run Castep with no XC potential, i.e. to place in the param file
XC_FUNCTIONAL = ZERO
Then reasonable agreement is obtained, although the wavefunction etc. will be nonsense.
Vasp
Vasp can produce a local potential in a LOCPOT if the
parameter lvpot is set to true in the INCAR
file. C2x gained the ability to read LOCPOT files in
version 2.34d. From Vasp 5.2.12 the default is to include the
exchange-correlation potential, but this can be excluded
if lvhar
is set to true.
An Example
A calculation was performed on a slab consisting of six layers of SiC, with hydrogen termination. The slab is polar. Running
c2x -c -D= --null -v SiC.check
reports
Total dipole (eA): (-0.026953,-0.000000,0.184254) magnitude 0.186215 Extra dipole field in V/A: (-0.024864,-0.000000,0.169974)
and running
c2x -E=3.5 -P='(0,0,0):(1,0,0):151' -i=0,1,1 SiC.check SiC.gnu
or
c2x -E=3.5 -P=a -i=0,1,1 SiC.check SiC.gnu
gives the following plot of potential in volts, averaged across the plane parallel to the slab, against position in Angstroms. A guide line of gradient -0.16997V/A has been added, so that one can see that the field in the vacuum region at the two ends of the cell is as predicted.
There are also examples of calculating potentials in nanotubes in several codes.