Delocalization Indices in Solids
Critic2 can calculate Bader’s localization and delocalization indices (DI) in solids using the pseudopotentials/plane waves approach. The DIs are a measure of electron delocalization (sharing) between atoms, related to the covalent bond order. There are two ways of calculating the DIs in solids using critic2. The simplest way makes use of a transformation of the oneelectron KS states to Wannier functions, in such a way that the resulting states are maximally localized. These maximally localized wannier functions (MLWF) are useful because they allow discarding most of the atomic overlap integrals required for the DI calculation. The details of the algorithm and some examples are described in JCTC 14 (2018) 4699. Alternatively, it is possible to calculate the DIs using Bloch states. This is much slower than using Wannier functions but requires fewer steps (because you do not need to use wannier90) and can be used for metals. An example is given below for metallic sodium.
In the following examples, we use Quantum
ESPRESSO (QE) to run the SCF
calculations, wannier90 to compute the
transformation to MLWF and the development version of critic2 to
obtain the DIs. The tool we need to extract the KS states from the QE
run (pw2critic.x
) was introduced in version ~6.4, so either this or
a more recent version is required. Any version of wannier90 from 2.0
onwards works.
Each DI calculation is carried out using a sequence of steps. In general, you need to:

Run a PAW SCF calculation with
pw.x
. This is done in order to calculate the allelectron density required for the correct determination of the atomic basins. The allelectron density is written to a cube file (rhoae.cube
) withpp.x
. 
Run a normconserving SCF calculation using
pw.x
and the sameecutrho
as in step 1. This will generate the KS states that we will transform into our MLWFs. 
Calculate the pseudovalence electron density from the converged normconserving SCF calculation and write it to a cube file (
rho.cube
) withpp.x
. This is done only for consistency checks as, in practice, the pseudodensity will always be available as a sum over the squares of the MLWFs. 
Use
open_grid.x
to unpack the kpoint symmetry from the normconserving calculation and prepare the Wannier run. Usingopen_grid.x
saves time but a nonselfconsistent calculation with a list of all kpoints in the uniform grid can also be used. 
Extract the KS coefficient, structure, and kpoint mapping data to a
.pwc
file usingpw2critic.x
. 
Run wannier90 in the usual way to get its checkpoint file (
.chk
). In spinpolarized cases, this needs to be done twice, once for each spin component. 
Read the structure and the
.pwc
and the.chk
files as a field in critic2 and calculate the DIs. The calculation of the DIs is activated using INTEGRABLE and the DELOC keyword.
In the output, the localization indices as well as all interatomic DIs
are given as a function of distance from the reference atom. For an
\(m \times n \times l\) grid, all DIs between all pairs of atoms
inside a \(m\times n\times l\) supercell will be obtained. Hence, the
higher m, n, and l, the longer the calculation will take. The cost of
the DI calculation also depends on the ecutrho
of both SCF
calculations. In molecular crystals, molecular localization and
delocalization indices are also calculated.
Magnesium Oxide (MgO)
The PAW calculation (step 1) with pw.x
is straightforward:
&control
title='crystal',
prefix='crystal',
pseudo_dir='../zz_psps',
/
&system
ibrav=0,
celldm(1)=1.0,
nat=2,
ntyp=2,
ecutwfc=80.0,
ecutrho=320.0,
/
&electrons
conv_thr = 1d8,
/
ATOMIC_SPECIES
o 15.999400 o_paw.UPF
mg 24.305000 mg_paw.UPF
ATOMIC_POSITIONS crystal
mg 0.500000000 0.000000000 0.000000000
o 0.000000000 0.500000000 0.500000000
K_POINTS automatic
4 4 4 1 1 1
CELL_PARAMETERS cubic
0.000000000 3.901537427 3.901537427
3.901540775 3.901536617 0.000000811
3.901540775 0.000000811 3.901536617
This is the (rhombohedral) primitive cell of magnesium oxide, with
only two atoms in it. After this calculation is converged, we extract
the allelectron density with pp.x
and the input:
&inputpp
prefix='crystal',
plot_num=21,
/
&plot
iflag = 3,
output_format=6,
fileout='rhoae.cube',
/
This generates rhoae.cube
. We will use this as the reference field
for our atomic integrations, since it is the only density that can
give the correct Bader basins. Next is the normconserving calculation
(step 2):
&control
title='crystal',
prefix='crystal',
pseudo_dir='../zz_psps,
outdir='.',
/
&system
ibrav=0,
celldm(1)=1.0,
nat=2,
ntyp=2,
ecutwfc=80.0,
ecutrho=320.0,
/
&electrons
conv_thr = 1d8,
/
ATOMIC_SPECIES
o 15.999400 o.UPF
mg 24.305000 mg.UPF
ATOMIC_POSITIONS crystal
mg 0.500000000 0.000000000 0.000000000
o 0.000000000 0.500000000 0.500000000
K_POINTS automatic
4 4 4 0 0 0
CELL_PARAMETERS cubic
0.000000000 3.901537427 3.901537427
3.901540775 3.901536617 0.000000811
3.901540775 0.000000811 3.901536617
There are several things to note in this input:

The pseudopotentials are naturally different from those in the PAW run but the crystal geometry is the same.

We use the
outdir
variable to have it write the.save
directory and the.wfc
files to the current directory, so there are no mishaps when we attempt to extract the KS states from them. 
The
ecutrho
is the same as in the PAW calculation. This will result in therhoae.cube
and the KS states being represented by grids with the same number of points (which is a prerequisite for the DI calculation in critic2). 
If symmetry is used in
pw.x
, then the density obtained from therho.cube
generated bypp.x
and that calculated as the sum of the squares of the MLWFs are slightly different. This is because the former undergoes an additional symmetrization step inside QE. 
We use an automatic uniform grid but we eliminate the kpoint shifts, which would result in kpoint positions that critic2 would not know how to handle.
The pseudodensity cube file is generated next, using:
&inputpp
prefix='crystal',
plot_num=0,
/
&plot
iflag = 3,
output_format=6,
fileout='rho.cube',
/
The rho.cube
file will only be used for checks, since the same
information is contained in the .pwc
file.
Now we run open_grid.x
(step 4) on this very simple input:
&inputpp
outdir='.',
prefix='crystal',
/
This will unpack the kpoints from the last SCF calculation and prepare
the wannier90 run. The execution of open_grid.x
generates a list of
kpoints that we will use later in the wannier90 input:
Writing output data file crystal_open.save/
Grid of qpoints
Dimensions: 4 4 4
Shift: 0 0 0
List to be put in the .win file of wannier90: (already in crystal/fractionary coordinates):
0.000000000000000 0.000000000000000 0.000000000000000 0.0156250000
0.000000000000000 0.000000000000000 0.250000000000000 0.0156250000
[...]
0.250000000000000 0.250000000000000 0.500000000000000 0.0156250000
0.250000000000000 0.250000000000000 0.250000000000000 0.0156250000
After unpacking the kpoint grid, we need to extract the KS state
coefficients from the QE files by using pw2critic.x
(step 5):
&inputpp
outdir = '.',
prefix='crystal_open',
seedname = 'mgo',
smoothgrid = .true.,
/
The crystal_open
prefix is automatically created by open_grid.x
from the original prefix by appending _open
. Note that pw2critic.x
needs to be run without MPI. Using mpirun
on it will not work. The
execution of pw2critic.x
creates a .pwc
file can be read by
critic2 and contains the reciprocalspace coefficients of the
KohnSham states. The smoothgrid
option ensures that the grid in the
.pwc
file has the same dimensions as the electron density cubes.
The second piece of information we need for the DI calculation is the
rotation of the KS states to yield the MLWF. We calculate this with
wannier90
(step 6). The input is:
num_wann = 4
num_iter = 20000
conv_tol = 1e4
conv_window = 3
begin unit_cell_cart
bohr
0.000000000 3.901537427 3.901537427
3.901540775 3.901536617 0.000000811
3.901540775 0.000000811 3.901536617
end unit_cell_cart
begin atoms_frac
mg 0.500000000 0.000000000 0.000000000
o 0.000000000 0.500000000 0.500000000
end atoms_frac
begin projections
random
end projections
search_shells = 24
kmesh_tol = 1d3
mp_grid : 4 4 4
begin kpoints
0.000000000000000 0.000000000000000 0.000000000000000 0.0156250000
[...]
0.250000000000000 0.250000000000000 0.250000000000000 0.0156250000
end kpoints
The num_wann
is the number of bands in the system, equal to the
number of KohnSham states
in the normconserving SCF output. The
geometry is the same as in the QE calculations. The kmesh is the same
as in the normconserving SCF calculation (4x4x4) and the list of
kpoints is the one generated in the output of open_grid.x
(shortened here for convenience). To run wannier90, first we generate
the list of items that it needs from QE by doing:
$ wannier90.x pp mgo.win
and then we run the pw2wannier90.x
utility on the input:
&inputpp
prefix='crystal_open',
seedname = 'mgo',
write_mmn = .true.,
write_amn = .true.,
outdir='.',
/
This creates the files containing the integrals required for the
wannier90 run. Note that we use the prefix for the files created by
open_grid.x
. Finally, we do:
$ wannier90.x mgo.win
and this runs the MLWF calculation and writes the rotation matrix to
the checkpoint (.chk
) file.
Now that we have the .pwc
file (with the KS coefficients) and the
.chk
file (with the orbital rotation), we can finally put everything
together and use critic2 to calculate the DIs (step 7). The input is:
crystal mgo.pwc
load rhoae.cube
load mgo.pwc mgo.chk
integrable 2
integrable 2 deloc
yt
The crystal structure is read from the .pwc
file (equivalently, it
can be read from any of the cube files or pw.x
input/output files.
They should all be the same.) The allelectron density is loaded first
and set as the reference fields so that it is used to calculate the
atomic basins. Then, the .pwc
and .chk
files are loaded together
into field number two, which is set as integrable and we activate the
calculation of the DIs with the DELOC keyword. Finally, we run the
YuTrinkle integration with YT.
After a few seconds, the DIs are done (it will take longer for larger crystals or crystals with more bands or denser density and kpoint grids). For each atom in the unit cell a table like this is written:
# Attractor 1 (cp=1, ncp=1, name=mg, Z=12) at: 0.5000000 0.0000000 0.0000000
# Id cp ncp Name Z Latt. vec.  Cryst. coordinates  Distance LI/DI
Localization index....................................................................... 0.01220619
26 2 2 o 8 0 1 0 0.0000000 0.5000000 0.5000000 3.9015366 0.10953819
58 2 2 o 8 1 1 0 1.0000000 0.5000000 0.5000000 3.9015366 0.10956477
40 2 2 o 8 1 0 1 1.0000000 0.5000000 0.5000000 3.9015366 0.10953816
8 2 2 o 8 0 0 1 0.0000000 0.5000000 0.5000000 3.9015366 0.10951576
2 2 2 o 8 0 0 0 0.0000000 0.5000000 0.5000000 3.9015408 0.10937694
64 2 2 o 8 1 1 1 1.0000000 0.5000000 0.5000000 3.9015408 0.10937112
27 1 1 mg 12 0 1 1 0.5000000 1.0000000 1.0000000 5.5176049 0.00114987
[...]
43 1 1 mg 12 1 3 1 1.5000000 3.0000000 1.0000000 13.5153222 0.00000393
85 1 1 mg 12 2 2 2 1.5000000 2.0000000 2.0000000 15.6061465 0.00000055
Total (atomic population)................................................................ 0.36083759
The numbers on the last column are the localization and delocalization
indices. The first is the LI, and then the DIs with the other atoms in
the system, ordered by distance to the atom that generates the
table. For a given row, the other fields are, in order: the
attractor ID, completelist ID, and nonequivalent list ID of the
atom, the atomic name, atomic number, and lattice vector translation
to its position relative to the main cell, its
crystallographic coordinates, and its distance to the atom that
generates the table. At the very end of the table, the sum of the
localization index and 0.5 times the sum of all the DIs is taken. This
number is the atomic population (average number of electrons in the
basin) and should coincide with the value given in the integration
table above (column $2
):
* Integrated atomic properties
# Id cp ncp Name Z mult Volume Pop Lap $2
1 1 1 mg 12  2.98869287E+01 1.55039497E+01 3.29194756E+01 3.60837590E01
2 2 2 o 8  8.88914848E+01 1.06091782E+01 3.29194756E+01 7.63916241E+00

Sum 1.18778413E+02 2.61131279E+01 4.81747975E12 8.00000000E+00
Final remark: because the DI calculation can take a long time,
checkpoint files are automatically generated that contain the atomic
overlap matrices (.pwcsij
) and the FAB integrals (.pwcfa
). These
files are automatically read in second and subsequent runs.
The files for this example can be found in the
example_11_10.tar.xz
package, mgo
subdirectory. The runit.sh
script automatizes the
steps above.
MgO Using a nonSCF Calculation Instead of open_grid.x
Exactly the same DI calculation can be carried out in a different way
by using a nonselfconsistent (nscf) calculation instead of the
open_grid.x
program. If the latter is available, there is really not
much point in carrying out the DI calculation in this manner, save for
testing purposes. The sequence of steps is the same as above except
instead of running open_grid.x
, we use pw.x
to run:
&control
title='crystal',
prefix='crystal',
pseudo_dir='../../data',
calculation='nscf',
wf_collect=.true.,
verbosity='high',
outdir='.',
/
&system
ibrav=0,
celldm(1)=1.0,
nat=2,
ntyp=2,
ecutwfc=80.0,
ecutrho=320.0,
/
&electrons
conv_thr = 1d8,
/
ATOMIC_SPECIES
o 15.999400 o.UPF
mg 24.305000 mg.UPF
ATOMIC_POSITIONS crystal
mg 0.500000000 0.000000000 0.000000000
o 0.000000000 0.500000000 0.500000000
K_POINTS crystal
64
0.00000000 0.00000000 0.00000000 1.562500e02
[...]
0.75000000 0.75000000 0.75000000 1.562500e02
CELL_PARAMETERS cubic
0.000000000 3.901537427 3.901537427
3.901540775 3.901536617 0.000000811
3.901540775 0.000000811 3.901536617
The calculation is nonselfconsistent (calculation='nscf'
) and uses
the converged wavefunction calculated in the selfconsistent step. In
addition, the list of kpoints is passed by hand, and it contains the
64 points corresponding to an unshifted 4x4x4 grid. This list can be
generated with the kmesh.pl
utility from the wannier90 package:
$ kmesh.pl 4 4 4
The rest of the DI calculation is exactly the same except all mentions
to the crystal_open
prefix (which correspond to the output of
open_grid.x
) are replaced by crystal
and the list of kpoints in
the wannier90 input must be the same as in the nonselfconsistent
calculation.
Graphite
The sequence of steps for calculating the DIs in more complex systems is exactly the same as in MgO. Because the Wannier transformation is illdefined in systems with partially occupied bands (metals), the calculation of DIs via Wannier functions does not work for those (see below for how to use Bloch states for this). However, semimetals such as graphite can be calculated without a problem.
In the graphite example included in the files package, we use a 8x8x2 kpoint grid. The resulting DIs clearly differentiate between DIs for atoms in the same graphene layer and in different layers:
+ Delocalization indices
Each block gives information about a single atom in the main cell.
First line: localization index. Next lines: delocalazation index
with all atoms in the environment. Last line: sum of LI + 0.5 * DIs,
equal to the atomic population. Distances are in bohr.
# Attractor 1 (cp=1, ncp=1, name=C, Z=6) at: 0.0000000 0.0000000 0.2500000
# Id cp ncp Name Z Latt. vec.  Cryst. coordinates  Distance LI/DI
Localization index....................................................................... 1.84476772
59 3 2 C 6 0 1 0 0.3333333 0.3333333 0.2500000 2.6795792 1.20502882
507 3 2 C 6 1 1 0 0.6666667 0.3333333 0.2500000 2.6795793 1.20297829
3 3 2 C 6 0 0 0 0.3333333 0.6666667 0.2500000 2.6795793 1.20496148
449 1 1 C 6 1 0 0 1.0000000 0.0000000 0.2500000 4.6411674 0.05704442
65 1 1 C 6 1 0 0 1.0000000 0.0000000 0.2500000 4.6411674 0.05704442
505 1 1 C 6 1 1 0 1.0000000 1.0000000 0.2500000 4.6411674 0.05696313
73 1 1 C 6 1 1 0 1.0000000 1.0000000 0.2500000 4.6411674 0.05696313
57 1 1 C 6 0 1 0 0.0000000 1.0000000 0.2500000 4.6411674 0.05679926
9 1 1 C 6 0 1 0 0.0000000 1.0000000 0.2500000 4.6411674 0.05679926
67 3 2 C 6 1 0 0 1.3333333 0.6666667 0.2500000 5.3591585 0.03623788
499 3 2 C 6 1 2 0 0.6666667 1.3333333 0.2500000 5.3591585 0.03538782
451 3 2 C 6 1 0 0 0.6666667 0.6666667 0.2500000 5.3591585 0.03534794
6 2 1 C 6 0 0 1 0.0000000 0.0000000 0.2500000 6.3268031 0.01795914
2 2 1 C 6 0 0 0 0.0000000 0.0000000 0.7500000 6.3268031 0.01796043
456 4 2 C 6 1 0 1 0.3333333 0.3333333 0.2500000 6.8708502 0.00629409
452 4 2 C 6 1 0 0 0.3333333 0.3333333 0.7500000 6.8708502 0.00629535
8 4 2 C 6 0 0 1 0.6666667 0.3333333 0.2500000 6.8708502 0.00629989
4 4 2 C 6 0 0 0 0.6666667 0.3333333 0.7500000 6.8708502 0.00629977
512 4 2 C 6 1 1 1 0.3333333 0.6666667 0.2500000 6.8708502 0.00627913
508 4 2 C 6 1 1 0 0.3333333 0.6666667 0.7500000 6.8708502 0.00627871
123 3 2 C 6 1 1 0 1.3333333 0.3333333 0.2500000 7.0895003 0.01039461
51 3 2 C 6 0 2 0 0.3333333 1.3333333 0.2500000 7.0895003 0.01034672
435 3 2 C 6 2 2 0 1.6666667 1.3333333 0.2500000 7.0895003 0.01030576
75 3 2 C 6 1 1 0 1.3333333 1.6666667 0.2500000 7.0895003 0.01038713
443 3 2 C 6 2 1 0 1.6666667 0.3333333 0.2500000 7.0895003 0.01031799
11 3 2 C 6 0 1 0 0.3333333 1.6666667 0.2500000 7.0895003 0.01035279
[...]
The first three atoms are the three inplane covalent bonds. Atoms in the same layer (z = 0.25) show distinctly higher DIs than atoms in different layers, even if they are farther away.
A Molecular Crystal: Urea
When a crystal composed of discrete units (a molecular crystal) is read into critic2, the program will automatically detect it and calculate the discrete fragment (molecules) that compose the crystal:
+ List of fragments in the system (2)
# Id = fragment ID. nat = number of atoms in fragment. Com = center of mass (bohr).
# Discrete = is this fragment finite?
# Id nat Center of mass Discrete
1 8 1.000000 0.500000 0.314384 Yes
2 8 0.500000 1.000000 0.685616 Yes
In the example used above, the urea crystal, there are two molecules in the unit cell each comprising 8 atoms.
The DI calculation in a molecular crystal follows the same sequence of steps as in MgO. However, critic2 detects the existence of discrete units and offers more information. Specifically, the program calculates the molecular localization indices (the sum of the LIs of all atoms in the molecule plus the intramolecular DIs):
* Integrated molecular properties
+ Localization indices
# Mol LI(A) N(A)
1 23.26310183 23.99999842
2 23.26310705 24.00000158
Then, critic2 also writes to the output the list of intermolecular delocalization indices. These are calculated as the sum of the DIs between all atoms in the two interacting molecules. There is a table of DIs for every molecule in the unit cell (two in this case) and the DI information is sorted by distance to this molecule, in much the same way as atomic DIs:
+ Delocalization indices
# Molecule 1 with 8 atoms at 1.000000 0.500000 0.314384
# Mol Latt. vec.  Center of mass (cryst)  Distance LI/DI
Localization index................................................... 23.26310183
2 1 0 0 1.5000000 1.0000000 0.6856165 8.1298272 0.16615726
2 0 0 0 0.5000000 1.0000000 0.6856165 8.1298272 0.16615840
2 1 1 0 1.5000000 0.0000000 0.6856165 8.1298272 0.16615814
2 0 1 0 0.5000000 0.0000000 0.6856165 8.1298272 0.16615731
1 0 0 1 1.0000000 0.5000000 0.6856165 8.8514772 0.46900724
2 1 1 1 1.5000000 0.0000000 0.3143835 9.2882471 0.07550946
2 1 0 1 1.5000000 1.0000000 0.3143835 9.2882471 0.07550956
2 0 0 1 0.5000000 1.0000000 0.3143835 9.2882471 0.07550951
2 0 1 1 0.5000000 0.0000000 0.3143835 9.2882471 0.07550956
1 1 0 0 0.0000000 0.5000000 0.3143835 10.5163259 0.01704341
1 0 1 0 1.0000000 0.5000000 0.3143835 10.5163259 0.01704340
1 1 0 1 0.0000000 0.5000000 0.6856165 13.7456087 0.00045912
1 0 1 1 1.0000000 0.5000000 0.6856165 13.7456087 0.00045912
1 1 1 0 0.0000000 0.5000000 0.3143835 14.8723308 0.00261252
1 1 1 1 0.0000000 0.5000000 0.6856165 17.3070757 0.00049917
Total (atomic population)............................................ 23.99999842
Note that the sum of LI plus half of all DIs is the average molecular electron population. In this case, both molecules are neutral, since they are equivalent by symmetry.
A SpinPolarized Case: FeO
Calculating the DIs in a spinpolarized calculation is a bit more convoluted, but follows essentially the same procedure. In this example, we calculate the DIs in iron(II) oxide (FeO), which has the same structure type as MgO (rocksalt) but is ferromagnetic. The PAW and NC SCF calculations and the extraction of the density cube files is done in the same way as in MgO, except that the calculation is spinpolarized. In PAW, this is done:
nspin=2,
starting_magnetization(1)=1.0,
occupations='smearing',
smearing='cold',
degauss=0.1,
and in the NC calculation it is essential that we have completely filled bands, so we fix the total magnetization based on the result of PAW:
nspin=2,
tot_magnetization=4.0,
This results in a total of 9 bands, as shown in the NC SCF output:
number of KohnSham states= 9
With verbosity=high
, we can verify that the spinup channel has 9
occupied bands and spindown has 5 occupied and 4 unoccupied bands:
[...]
 SPIN UP 
k = 0.0000 0.0000 0.0000 ( 1639 PWs) bands (ev):
8.7938 7.7292 7.8322 7.9838 8.1232 8.1283 8.1355 8.1770
8.1880
occupation numbers
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000
 SPIN DOWN 
k = 0.0000 0.0000 0.0000 ( 1639 PWs) bands (ev):
8.2542 8.7307 8.7479 8.7545 11.2368 11.3607 11.4147 12.2011
12.4321
occupation numbers
1.0000 1.0000 1.0000 1.0000 1.0000 0.0000 0.0000 0.0000
0.0000
[...]
An important difference relative to MgO happens when we calculate the
MLWFs. The way wannier90 works is that it requires two different
executions, one for each spin channel, where we need to indicate how
many bands are occupied. The spinup wannier90 input
(feo_up.win
) is simple because we want to use all the occupied bands
in the transformation. It contains:
num_wann = 9
num_bands = 9
[...]
and the corresponding pw2wannier90.x
file to generate the integrals
for wannier90 is:
&inputpp
prefix='crystal_open',
seedname = 'feo_up',
spin_component="up",
write_mmn = .true.,
write_amn = .true.,
/
The wannier90 calculation is run in the same way as before:
$ wannier90.x pp feo_up.win
$ pw2wannier90.x < feo.pw2wan.up.in  tee feo.pw2wan.up.out
$ wannier90.x feo_up.win
which generates the feo_up.chk
checkpoint file. We will read this
file into critic2.
The spindown case is a little more complicated because only a subset
of the bands are filled. There are two different spindow wannier90
inputs: one for before and one for after pw2wannier90.x
. The first
input is:
num_wann = 5
num_bands = 9
exclude_bands : 69
where we indicate that, even though there are 9 bands in the calculation, we want to exclude bands 6 to 9, which are unoccupied. We write the files that request the calculation of the integrals in QE:
$ wannier90.x pp feo_dn.win
and then we run pw2wannier90.x
with the input:
&inputpp
prefix='crystal_open',
seedname = 'feo_dn',
spin_component="down",
write_mmn = .true.,
write_amn = .true.,
/
In the second wannier run, where the actual MLWFs are calculated, we only have 5 bands available, so it contains:
num_wann = 5
num_bands = 5
To generate the MLWFs
$ pw2wannier90.x < feo.pw2wan.dn.in  tee feo.pw2wan.dn.out
$ wannier90.x feo_dn.win
This creates the feo_dn.chk
checkpoint file with the rotation for
the spindown MLWFs.
Finally, we indicate that this is a spinpolarized case by passing both checkpoint files to critic2:
crystal feo.pwc
load rhoae.cube id rho
load feo.pwc feo_up.chk feo_dn.chk
integrable 2
integrable 2 deloc
yt
The interpretation of the critic2 output is essentially the same as in the MgO case.
A Metal: Elemental Sodium
The DIs in systems that have partially filled bands (metals) cannot be calculated with maximally localized Wannier functions at present, because the Wannier transformation is illdefined. In this case, the Bloch states resulting from the SCF calculation are used directly. This results in a substantially higher computational cost of the DI calculation, but the calculation process is simpler because wannier90 needs not be used. The DIs can always be calculated using Bloch states, including the examples above, although there is no reason to do so if MLWFs can be used.
The steps to calculate DIs using Bloch states are:

Run a PAW SCF calculation with
pw.x
to obtain the allelectron density (rhoae.cube
). 
Run a normconserving SCF calculation using
pw.x
and the sameecutrho
as in step 1. This will generate the KS states. This calculation needs to be run without kpoint symmetry (i.e. withnosym=.true.
andnoinv=.true.
). 
Use
pw2critic.x
withsmoothgrid=.true.
to obtain the.pwc
file. 
The
.pwc
file is read in critic2 with an input like:
crystal na.pwc
load rhoae.cube
load na.pwc
integrable 2
integrable 2 deloc
bader
Note that the checkpoint files in this case are absent, but the input is otherwise the same.
DIs from a Checkpoint File
Carrying out the DI integration generates two checkpoint files: one
containing the atomic overlap matrices (with extension .pwcsij
) and
another with the exchangecorrelation atomic integrals (extension
.pwcfa
). Either of them can be used to regenerate the DIs simply
and quickly. The fa
files is actually calculated from the sij
file.
Because it tends to be big in some cases it is of interest to delete
the .pwc
file after a successful DI calculation. If you keep either
of the checkpoints, then it is possible to regenerate the DIs without
the .pwc
or the wannier90 checkpoint files. You do this by doing:
crystal rhoae.cube
load rhoae.cube
integrable deloc_sijchk na.pwcsij
bader
where the integrable property points to the checkpoint file
directly. The fa
checkpoint file can be used instead:
integrable deloc_fachk na.pwcfa
Example Files Package
Files: example_11_10.tar.xz.