# Integration of Atomic Charges and Volumes with Grids

In this example, we examine how to compute Bader atomic charges and atomic volumes in periodic solids using different quantum chemistry software that uses plane-waves/pseudopotentials approach. For the integration, we will use the Yu-Trinkle method (YT keyword) but the same can be done with Henkelman et al.’s method (BADER keyword).

As example, we will use the urea crystal. When it crystallizes, urea forms an orthorhombic molecular crystal with 16 atoms in the unit cell (2 molecules), shown on the right. There are chains of double hydrogen bonds all along the c axis, and all molecules in the crystal are equivalent by symmetry.

The input files necessary for this example can be generated from critic2’s internal library of structures. For instance, to create a Quantum ESPRESSO input, we do:

```
crystal library urea
write urea.scf.in
```

This generates an input file *template*. The particular keywords this
template has are most likely not what you want, but the structure,
which is after all the difficult part to write, is correct. Likewise,
we can generate inputs for VASP or abinit with:

```
crystal library urea
write POSCAR
write urea.abin
```

Output formats for many other programs are supported as well.

## Quantum ESPRESSO

### How it Works

Quantum ESPRESSO does DFT
calculations in periodic solids using plane waves and
pseudopotentials. A typical calculation of the atomic charges has two
steps. In the first step, you carry out the SCF calculation with the
`pw.x`

program on a `.scf.in`

input file, such as the one generated in
the example above. This generates the converged wavefunction, which is
stored in QE’s internal format.

To extract the densities you use the `pp.x`

(post-process) program to
ask QE to write the density on a three-dimensional grid. This is not
an arbitrary choice of representation for the density: In a
plane-waves code the density and wavefunctions are expressed naturally
on uniform three-dimensional grids. `pp.x`

generates two types of
file formats for the 3d density grid: Gaussian cube files (`.cube`

)
and xcrysden’s `.xsf`

files. We will use the former, but critic2 can
also understand and use `.xsf`

files. Note that both cube and xsf
files contain also the structural information, so you can use them to
read the crystal structure into critic2 as well as the density at any
time.

There are several ways of calculating the atomic charges, depending on the type of pseudopotential you use. Plane-wave/pseudopotentials codes use pseudopotentials for two reasons: 1) remove the core electrons and 2) smooth out the system’s wavefunctions close to the atoms, in order to prevent having to use an extremely large number of plane waves to represent the oscillations close to the atomic nuclei. Because of this the density peaks close to the nuclei are not as high as they should be and, in systems with significant charge transfer, the maxima at the nuclei may be missing altogether. Therefore, the electron smoothed density used in QE (the pseudo-density) cannot be used to calculate the QTAIM atomic basins, and we must have a way to recover the all-electron density from the calculation. Depending on the method and pseudopotential employed, this may be easy or hard.

Ideally, you want to use the PAW method because the information about
the conversion to the all-electron wavefunctions is not lost, and
therefore you keep the option of reconstructing the all-electron
density. In recent versions of QE, this option is given by the
`plot_num=21`

option in `pp.x`

. Here is an example input for `pp.x`

that writes the reconstructed all-electron density to a cube file
(`rhoae.cube`

):

```
&inputpp
prefix='urea',
outdir='.',
plot_num=21,
/
&plot
iflag=3,
output_format=6,
fileout='rhoae.cube',
/
```

If your version of QE is not very recent, then `plot_num=21`

may not
be available to you. In that case, you can still get the all-electron
*valence* density with `plot_num=17`

, and then augment this density
with the core contribution calculated using critic2’s internal density
tables (see below).

The all-electron density is too steep close to the nuclei to be
efficiently represented by a grid and, consequently, it does not sum
to the correct number of electrons. Therefore, even though it gives
the correct atomic basins, integrating the all-electron density to
obtain the atomic charges is a poor idea. Instead, we integrate the
pseudo-density, which is the density actually used in the SCF
calculation and it is normalized to the number of valence
electrons. To obtain the pseudo-density from `pp.x`

, you use the same
input as above but with `plot_num=0`

. In the case of a non-PAW
calculation, using either norm-conserving or ultrasoft
pseudopotentials, `plot_num=0`

will be the only option available to
you and, in that case, the pseudo-density is used to calculate the
shape of the atomic basins as well. Although not ideal, experience has
shown that the calculated atomic charges are about the same as with
the correct all-electron density, provided core augmentation is
utilized.

To summarize, to calculate atomic volumes and charges you do:

### PAW calculation

If you ran a PAW calculation, get the all-electron density
(`plot_num=21`

, `rhoae.cube`

) and the pseudo-density (`plot_num=0`

,
`rho.cube`

) and integrate the latter in the basins of the former
with:

```
crystal rhoae.cube
load rhoae.cube
load rho.cube
integrable 2
yt
```

The basins are given by the reference field (`rhoae.cube`

, the first
field loaded) and we define the second field (`rho.cube`

) as the
integrand. The YT keyword launches the calculation of atomic
properties.

### PAW calculation (old QE version)

If you ran a PAW calculation but your QE version is old and you do
not have `plot_num=21`

, then write the reconstructed valence density
(`plot_num=17`

, `rhoval.cube`

) and the pseudo-density (`plot_num=0`

,
`rho.cube`

) and integrate the latter in the basins of the former,
but using core augmentation to account for the missing core
electrons:

```
crystal rhoval.cube
load rhoval.cube core zpsp c 4 o 6 n 5 h 1
load rho.cube
integrable 2
yt
```

Core augmentation is activated with the CORE option to LOAD, and we
need to give to critic2 the pseudopotential charges for each
atom. This value is equal to the atomic number minus the number of
core electrons removed in the calculation. In QE, this value can be
obtained by summing the numbers in the `occ`

column under ```
Valence
configuration
```

in the corresponding `.UPF`

file. The core
contribution for the given number of core electrons is added by
critic2 every time the field is evaluated using critic2’s internal
density tables.

### NC or US pseudopotential calculation

If you ran a US or NC pseudopotential calculation, then you only
have the pseudodensity (`plot_num=0`

, `rho.cube`

). In this case,
load the pseudo-density twice and activate core augmentation the
first time. Use the core-augmented field to generate the atomic
basins and the non-augmented field as the integrand:

```
crystal rho.cube
load rho.cube core zpsp c 4 o 6 n 5 h 1
load rho.cube
integrable 2
yt
```

## VASP

To calculate Bader charges in VASP, you need to run the SCF
calculation in the usual way, but including the following tag in the
`INCAR`

file:

```
LAECHG = .TRUE.
```

This will generate two additional files: the `AECCAR0`

, containing the
core density, and the `AECCAR2`

, containing the reconstructed valence
density. The sum of the two needs to be used to generate the atomic
basins, while the pseudo-density, contained in the `CHGCAR`

file, is
integrated in them.

These operations can be all done within critic2 cheaply and without generating any additional files, using the following input:

```
crystal CHGCAR
load AECCAR0
load AECCAR2
load as "$1+$2"
reference 3
load CHGCAR
integrable 4
yt
```

As with the cube files, any of the `CHGCAR`

, `AECCAR0`

, etc. files
specifies the crystal structure inside, so you can have it critic2
read it from any of them. (You can also use the `POSCAR`

and the
`POTCAR`

.) We load the `AECCAR0`

as field 1, and the `AECCAR2`

as
field 2, and create field 3 as the sum of those two files. Field 3 now
contains the all-electron density, and we set it as reference in order
to have it generate the atomic basins. Finally, we load the `CHGCAR`

file as field 4 and mark it as an integrand for our YT calculation.

## Abinit

If you are using PAW datasets, then the reconstructed valence density can be obtained in abinit using:

```
prtden 2
```

in the input file. This creates two density files: the pseudo-density
(with suffix `_DEN`

) and the reconstructed valence density
(`_PAWDEN`

). Currently, there is no way to obtain the all-electron
density from abinit directly. In the case of a non-PAW calculation,
using either norm-conserving or ultrasoft pseudopotentials, then
`prtden 2`

is not available and you have to use:

```
prtden 1
```

which writes the `_DEN`

file only.

If you have the `_PAWDEN`

file, the calculation of the atomic charges
is straightforward:

```
crystal urea_o_DEN
load urea_o_PAWDEN core zpsp h 1 c 4 n 5 o 6
load urea_o_DEN
integrable 2
yt
```

The `_DEN`

and `_PAWDEN`

files both contain the structural information
about the crystal, so they can be used to read the crystal structure
in critic2 through the CRYSTAL keyword. The `_PAWDEN`

density is
loaded as field 1 and is augmented with core electrons corresponding
to the pseudopotential charges (atomic number minus number of replaced
core electrons). The all-electron density obtained in this way is used
to generate the atomic basins. The pseudo-density (`_DEN`

) is then
loaded as the second field and integrated inside the basins.

If you ran a non-PAW calculation, then replace the `_PAWDEN`

file with
the `_DEN`

file in the input above.

## Output for urea

The outputs from all the calculations above are very similar. The main table containing the results of the integration appears right at the end of the output:

```
* Integrated atomic properties
# (See key above for interpretation of column headings.)
# Integrable properties 1 to 4
# Id cp ncp Name Z mult Volume Pop Lap $2
1 1 1 C_ 6 -- 3.22216552E+01 4.21263334E+00 1.02937536E-01 2.19098664E+00
2 2 1 C_ 6 -- 3.22216552E+01 4.21263334E+00 1.02937536E-01 2.19098664E+00
3 3 2 H_ 1 -- 2.00349124E+01 5.13744359E-01 -3.64792425E-04 5.13705108E-01
4 4 2 H_ 1 -- 2.00349124E+01 5.13744359E-01 -3.64792425E-04 5.13705108E-01
5 5 2 H_ 1 -- 2.00349124E+01 5.13744359E-01 -3.64792425E-04 5.13705108E-01
6 6 2 H_ 1 -- 2.00349124E+01 5.13744359E-01 -3.64792425E-04 5.13705108E-01
7 7 3 H_ 1 -- 2.05322160E+01 5.11380307E-01 -1.96443849E-02 5.11335066E-01
8 8 3 H_ 1 -- 2.05322160E+01 5.11380307E-01 -1.96443849E-02 5.11335066E-01
9 9 3 H_ 1 -- 2.05322160E+01 5.11380307E-01 -1.96443849E-02 5.11335066E-01
10 10 3 H_ 1 -- 2.05322160E+01 5.11380307E-01 -1.96443849E-02 5.11335066E-01
11 11 4 O_ 8 -- 1.18754676E+02 9.40395235E+00 -1.26560690E-01 7.23176246E+00
12 12 4 O_ 8 -- 1.18754676E+02 9.40395235E+00 -1.26560690E-01 7.23176246E+00
13 13 5 N_ 7 -- 1.26002600E+02 8.28697652E+00 3.18207545E-02 6.26367524E+00
14 14 5 N_ 7 -- 1.26002600E+02 8.28697652E+00 3.18207545E-02 6.26367524E+00
15 15 5 N_ 7 -- 1.26002600E+02 8.28697652E+00 3.18207545E-02 6.26367524E+00
16 16 5 N_ 7 -- 1.26002600E+02 8.28697652E+00 3.18207545E-02 6.26367524E+00
------------------------------------------------------------------------------------------------
Sum 9.68231575E+02 6.44815761E+01 -1.12894416E-12 4.80003598E+01
```

In this table, critic lists the 16 attractors it found, which in this
case are all associated to atoms in the system (there are no
non-nuclear maxima). In addition, it gives the following columns: the
calculated atomic volumes (`Volume`

), the all-electron density
integrated in its basins (`Pop`

), the Laplacian of the all-electron
density (`Lap`

), and the integral of field number 2 (the
pseudo-density) in the all-electron basins (`$2`

). The last column
gives the valence electron population for each atom. The atomic charge
equals the pseudopotential charge minus the valence electron
population. For instance, the charge for carbon in this system would
be \(Z_{\rm psps} - N = 4 - 2.19 = 1.81\).

Because critic2 detects that urea is a molecular crystal, it will also helpfully list for you the integrated properties of the molecules as a whole:

```
* Integrated molecular properties
# (See key above for interpretation of column headings.)
# Integrable properties 1 to 4
# Mol Volume Pop Lap $2
1 4.84115787E+02 3.22407881E+01 -1.12265058E-12 2.40001799E+01
2 4.84115787E+02 3.22407881E+01 -6.32133235E-15 2.40001799E+01
```

The interpretation of the columns is the same, but this time the numbers refer to each of the two urea molecules in the unit cell. Because they are equivalent by symmetry, their volumes and electron populations are the same and the charge is zero (the sum of the pseudopotential charges in each molecule is 24).

## Example files package

Files: example_11_01.tar.xz. Run the examples as follows:

- QE/PAW calculation (
`qe_paw`

):`## generate the pseudopotentials first ld1.x < h.in ld1.x < c.in ld1.x < n.in ld1.x < o.in ## run the SCF calculation pw.x < urea.scf.in > urea.scf.out ## get the densities pp.x < urea.rho.in > urea.rho.out pp.x < urea.rhoval.in > urea.rhoval.out pp.x < urea.rhoae.in > urea.rhoae.out ## run critic2 critic2 urea.cri urea.cro`

If the

`rhoae.cube`

was not generated, you are using an old version of QE. In that case, edit`urea.cri`

, comment out the first block and uncomment the second block. - QE/NC or US calculation (
`qe_nc`

and`qe_us`

):`## generate the pseudopotentials first ld1.x < h.in ld1.x < c.in ld1.x < n.in ld1.x < o.in ## run the SCF calculation pw.x < urea.scf.in > urea.scf.out ## get the density pp.x < urea.rho.in > urea.rho.out ## run critic2 critic2 urea.cri urea.cro`

- VASP calculation: first, create the POTCAR file by concatenating the
POTCARs for C, H, O, and N, in that order. Then:
`vasp critic2 urea.cri urea.cro`

- abinit:
`abinit ## or abinis/abinip if you are running a really old version critic2 urea.cri`