• Arithmetic Operations with Grids

    Introduction

    In critic2, scalar fields such as the electron density or the ELF can be loaded using the LOAD keyword, and then manipulated using arithmetic expressions. In the case of scalar fields given as grids (e.g. from a cube file), new grids can be generated by operating and combining existing grids. The grids generated are automatically integrated over the unit cell if the structure is a periodic crystal. This example presents several cases demonstrating the usefulness of this feature.

    Calculating Energy Components in a Quantum ESPRESSO Calculation

    The Quantum ESPRESSO code reports the different components of the total energy (kinetic energy, Ewald, exchange-correlation,…) when the pw2casino option is used in the course of a calculation with norm-conserving pseudopotentials. For instance, this is an input file for rutile (TiO2) at its experimental geometry:

    &control
 title='crystal',
 prefix='crystal',
 pseudo_dir='.',
/
&system
 ibrav=0,
 nat=6,
 ntyp=2,
 ecutwfc=60.0,
 ecutrho=240.0,
/
&electrons
 conv_thr = 1d-8,
/
ATOMIC_SPECIES
Ti   47.867000 ti.UPF
O    15.999400 o.UPF

ATOMIC_POSITIONS crystal
Ti    0.00000000     0.00000000     0.00000000
Ti    0.50000000     0.50000000     0.50000000
O    0.30530000     0.30530000     0.00000000
O    0.69470000     0.69470000     0.00000000
O    0.19470000     0.80530000     0.50000000
O    0.80530000     0.19470000     0.50000000

K_POINTS automatic
2 2 3 0 0 0

CELL_PARAMETERS bohr
        8.680891628421     0.000000000000     0.000000000000
        0.000000000000     8.680891628421     0.000000000000
        0.000000000000     0.000000000000     5.590036668212

    where the k-point grid has been made deliberately small for ease of calculation. This input is then run with pw.x:

    pw.x -pw2casino < rutile.scf.in > rutile.scf.out

    which produces a decomposition of the total energy:

     Energies determined by pw2casino tool
 -------------------------------------
 Kinetic energy      55.812056357004273       au  =     111.62411271400855       Ry
 Local energy       -40.509441848372816       au  =    -81.018883696745633       Ry
 Non-Local energy   -29.658425859702600       au  =    -59.316851719405200       Ry
 Ewald energy       -58.019117159310184       au  =    -116.03823431862037       Ry
 xc contribution    -16.583656552046751       au  =    -33.167313104093502       Ry
 hartree energy      18.029992555875470       au  =     36.059985111750940       Ry
 Total energy       -70.928592506552619       au  =    -141.85718501310524       Ry

    We now generate the cube files necessary for recalculating these energy components within critic2 (except for the non-local energy, which requires information from the pseudopotential files not available to critic2). First, we need the (pseudo) electron density, which can be obtained with the input:

    &inputpp
 prefix='crystal',
 plot_num=0,
 outdir='.',
/
&plot
 iflag = 3,
 output_format=6,
 fileout='rho.cube',
/

    and using pp.x:

    pp.x < rutile.rho.in > rutile.rho.out

    This operation generates a rho.cube file which can be loaded into critic2. Similarly, we generate cube files for the kinetic energy density (plot_num=22, rutile.kin.in, kin.cube), bare ionic component of the potential (plot_num=2, rutile.vlocal.in, vlocal.cube), and the bare ionic plus Hartree components of the potential (plot_num=11, rutile.vel.in, vel.cube).

    Now we write the critic2 input file to read all the cube files, combine them, and calculate the corresponding integrals over the unit cell. First, we read the crystal structure:

    CRYSTAL rutile.scf.in

    (Alternatively, we could have read the structure from any of the cube files; the result would be the same.) Then, we read the electron density:

    LOAD rho.cube ID rho

    This loads the grid as scalar field number 1. Critic2 reports in the output:

    + Field number 1
  Name: rho.cube
  Source: rho.cube
  Type: grid
  Grid dimensions : 48  48  30
[...]
  Cell integral (1) = 31.99999600
[...]
  Alias for this field (2): $1, $rho

    The output contains the field number (1, because it is the first field we loaded), the file it was read from, the field type (in this case, a grid), the grid dimensions (number of points) and the cell integral. The integral over the unit cell is simply the sum of the field over the whole grid times the cell volume divided by the number of points. Naturally, the density integrates to the total number of valence electrons (32 for this system). Lastly, critic2 reports the list of names (alias) which can be used to address this field: $1 (because it was the first field loaded) and $rho. This last alias was given by us using the ID option to LOAD.

    Similarly, loading the kinetic energy density:

    LOAD kin.cube ID kin

    gives:

      Cell integral (2) = 111.62416577

    which is the pw2casino kinetic energy in Rydberg (there is a slight disagreement because of how this field is calculated in pp.x). The bare ionic potential can be loaded as:

    LOAD vlocal.cube ID vlocal

    and the corresponding energy is calculated with:

    LOAD AS "$rho * $vlocal"

    This last command loads a grid taken as the product of the electron density times the local potential. Note the use of the corresponding aliases (rho for the density and vlocal for the potential) that match the corresponding ID options when those fields were loaded. The integral is the “local energy”:

      Cell integral (4) = -81.01892005

    which matches the corresponding result from the pw2casino output (in Rydberg).

    Next we load the Hartree plus bare ionic potential from pp.x:

    LOAD vel.cube ID vel

    To calculate the Hartree contribution only, we need to subtract the bare ionic potential first, and the corresponding energy would be 1/2 of the integral. Therefore:

    LOAD AS "0.5 * $rho * ($vel - $vlocal)"

    which gives:

      Cell integral (6) = 36.06001905

    again matching the pw2casino Hartree energy in Rydberg.

    The Ewald energy contribution can be calculated using the EWALD keyword, which requires specifying the atomic charges for each atom in the system. In this calculation, the Ti pseudopotential provides 4 valence electrons and the O pseudopotential 6 electrons. We first indicate these charges to critic2 using the Q keyword, and then calculate the Ewald energy:

    Q Ti 4 O 6
EWALD

    which gives:

    + Electrostatic potential at atomic positions
#id name mult    charge         Vel(Ha/e)
1    Ti   2     4.000000E+00  -2.7190702464E+00
2    O    4     6.000000E+00  -3.9285696807E+00

* Ewald electrostatic energy (Hartree) = -5.801911715453E+01

    which agrees with pw2casino.

    Lastly, the exchange-correlation component can be calculated using the LIBXC library integrated into critic2. This requires making the library available to critic2 in the compilation step. This calculation used PBE, which is a GGA functional, so we will need the norm of the gradient of the electron density. To calculate it, we use:

    LOAD AS FFT GMOD rho ID grad

    This command loads a new field (with alias grad) as the gradient norm of rho calculated using a fast Fourier transform, the exact same operation as used in Quantum ESPRESSO and similar programs. To obtain the list of functionals provided by the LIBXC library we use the LIBXC command:

    echo LIBXC | critic2 -q

    which gives a very long list of functionals and their associated codes. For PBE:

    * LIST of libxc functionals
1     lda_x                       x    LDA
2     lda_c_wigner                c    LDA
[...]
101   gga_x_pbe                   x    GGA
[...]
130   gga_c_pbe                   c    GGA

    and therefore we need codes 101 for PBE exchange and 130 for PBE correlation (the codes may change depending on the version of LIBXC used). The energy density from each of these components can be accessed in critic2 with the xc() function, which takes as the arguments the electron density, the norm of the density gradient, and the functional code. To build the PBE energy density, use:

    LOAD AS "2*(xc($rho,$grad,101)+xc($rho,$grad,130))" ID exc

    where the factor of 2 is used to convert into Rydberg. The result is:

      Cell integral (8) = -33.16729549

    which also matches the pw2casino result. (Important note: the energy calculated in this way only matches pw2casino if the non-linear core correction (NLCC) is not used in the pseudopotential generation. As in the case of the non-local energy, calculating the NLCC contribution requires information not available to critic2.)

    Calculating Crystal Void Volumes and Electron Populations

    We consider now how to estimate the size of the crystal voids and the number of electrons in them. For this, we use the same structure as in the previous example (rutile). First, the structure is read:

    CRYSTAL rutile.scf.in

    and then we define crystal voids, somewhat arbitrarily, as the part of the unit cell where the sum of atomic densities (the promolecular density) is less than 0.01 atomic units. To calculate the volume of this region, we load the field:

    LOAD AS "$rho0 < 0.01" 48 48 30

    Inside the arithmetic expression, $rho0 is the sum of the atomic densities, which is always available to critic2 provided a structure has been read, either as $rho0 or simply as $0. The field loaded in the previous command has zeros at the positions where the promolecular density is higher than 0.01 and ones where it is lower. The last three numbers (48 48 30) are used to create a uniform grid with those dimensions, so that we can integrate easily it over the unit cell and obtain the volume of the crystal voids. These values are the same as in the cube files generated by pp.x in the previous example, but others could have been used and they can be increased if more accuracy is required. The integral of the grid defined above is read from the output:

      Cell integral (1) = 1.92586872

    and, since the field has 1s where the density is low and 0s where it is high, this value corresponds to the volume of the region where the promolecular density is lower than 0.01. Compare this result with the cell volume, also from the output:

      Cell volume (bohr^3): 421.25331

    Therefore, about 0.45% of the cell is occupied by voids, according to our definition. Similarly, we can calculate the number of electrons inside the voids by doing:

    LOAD AS "($rho0 < 0.01)*$rho0" SIZEOF 1

    In this case, we multiply the comparison expression $rho0 < 0.01 by the value of the promolecular density. At the points where the density is higher than 0.01, the expression in parenthesis is zero, and therefore this new grid is also zero. If the density is lower than 0.01, the parenthesis expression equals one and therefore the new field has the same value as the promolecular density. The integral of this field:

      Cell integral (2) = 0.01890330

    gives the number of electrons inside the crystal voids. The SIZEOF keyword is used to copy the number of grid points from the first field we loaded (48 48 30).

    If a proper calculated electron density is available, we can use it to repeat the same calculation with it instead of the promolecular density. First, we read the density cube:

    LOAD rho.cube ID rho

    and then we calculate the crystal void volume and number of electrons:

    LOAD AS "$rho < 0.01"
LOAD AS "($rho < 0.01)*$rho"

    resulting in:

      Cell integral (4) = 47.64696721
  Cell integral (5) = 0.36062947

    When a self-consistent calculation is carried out, the electron density migrates from the crystal voids to the atomic and bonding regions, and the crystal voids are therefore much larger.

    Example Files Packages

    Files: example_02_01.tar.xz.

    Updated: