Loading Scalar Fields

Loading a Field (LOAD)

Scalar fields are loaded with the LOAD keyword. A number of different field formats are supported:

LOAD file.cube
LOAD file.bincube
LOAD file_{DEN|PAWDEN|ELF|ELF|POT|VHA|VHXC|VXC|GDEN1|
           GDEN2|GDEN3|LDEN|KDEN}
LOAD [file.]{CHGCAR|CHG|AECCAR0|AECCAR1|AECCAR2} [block.i|RHO|SPIN|MAGX|MAGY|MAGZ]
LOAD {[file.]ELFCAR} [block.i|RHO|SPIN|MAGX|MAGY|MAGZ]
LOAD file.qub
LOAD file.xsf
LOAD file_fmt
LOAD file.grid
LOAD file.{clmsum|clmup|clmdn} file.struct
LOAD file.{RHO,BADER,DRHO,LDOS,VT,VH}
LOAD file.001
LOAD file.OUT
LOAD STATE.OUT GEOMETRY.OUT
LOAD STATE.OUT GEOMETRY.OUT OTHER.OUT
LOAD file1.ion {nat1.i/at1.s} file2.ion ...
LOAD file.xml file.bin file.hsd
LOAD file.wfn
LOAD file.wfx
LOAD file.fchk [READVIRTUAL]
LOAD file.molden [READVIRTUAL] [ORCA|PSI4]
LOAD file.molden.input [READVIRTUAL]
LOAD file.pwc [file.chk [filedn.chk]] [SPIN spin.i] [KPT k1.i k2.i...]
              [BAND b1.i b2.i ...] [ERANGE emin.r emax.r]
LOAD COPY id.s [TO id2.s]
LOAD PROMOLECULAR
LOAD PROMOLECULAR [FRAGMENT file.xyz]
LOAD [WIEN|ELK|PI|CUBE|BINCUBE|ABINIT|VASP|VASPNOV|QUB|XSF|FMT|ELKGRID|SIESTA|FPLO|DFTB|
      WFN|WFX|MOLDEN|MOLDEN_ORCA|MOLDEN_PSI4|FCHK|PWC] file
LOAD ... [NEAREST|TRILINEAR|TRISPLINE|TRICUBIC|SMOOTHRHO [NENV nenv.i] [FDMAX fdmax.r]]
         [EXACT|APPROXIMATE] [RHONORM|VNORM] [NOCORE] [NUMERICAL|ANALYTICAL]
         [TYPNUC {-3,-1,1,3}] [NORMALIZE n.r] [{NAME|ID} id.s]
         [NOTESTMT] [ZPSP at1.s q1.r...]
LOAD AS "expression.s" [n1.i n2.i n3.i|SIZEOF id.s|GHOST]
LOAD AS PROMOLECULAR {n1.i n2.i n3.i|SIZEOF id.s}
        [FRAGMENT file.xyz]
LOAD AS CORE {n1.i n2.i n3.i|SIZEOF id.s} {ZPSP at1.s q1.r ...}
LOAD AS FFT [GX|GY|GZ|HXX|HXY|HXZ|HYY|HYZ|HZZ|GMOD|LAP|POT] id.s
LOAD AS RESAMPLE id.s n1.i n2.i n3.i
LOAD AS CLM {ADD id1.s id2.s|SUB id1.s id2.s}

Critic2 loads scalar fields in “field slots”: integer identifiers that represent the field throughout the run. The NAME (or ID) keyword can be used to assign a string identifier to a field as well, which is useful in cases when many fields are loaded at the same time. At the beginning of the run, field number 0 is the promolecular density, which is automatically loaded right after the structure is succesfully read. Successive LOAD commands assign field numbers in increasing order: the first LOAD will create field number 1, the second LOAD will be field number 2, and so on. By default, the first loaded field (other than the promolecular density in slot 0) becomes the reference field (see REFERENCE).

The simplest usage of load is:

LOAD file.ext

The scalar field information is read from file.ext. Critic2 uses the extension to decide which format should be used for the reading:

.cube for Gaussian cube files (grid field). Cube files can also be generated by programs other than Gaussian, for instance, Quantum ESPRESSO. They may correspond to a molecule or a periodic crysatl.
Binary cube files (.bincube) contain the same information as cube files but use a binary format. They can be generated by critic2.
_DEN, _PAWDEN, _ELF,… for abinit grids.
CHGCAR, AECCAR*, CHG, or ELFCAR for VASP grids.
.qub for aimpac grids.
.xsf for xcrysden grids.
_fmt for CASTEP _fmt grids.
.grid for elk’s grids (this requires a patch to elk, e-mail to request it).
.clmsum, .clmup, and .clmdn for WIEN2k.
.RHO, .BADER, .DRHO, .LDOS, .VT, and .VH for SIESTA’s grids.
.001 for FPLO grids.
.OUT for elk’s STATE.OUT, ELF3D.OUT, RHO3D.OUT, etc.
.ion for aiPI ion files.
.xml for DFTB+’s detailed.xml.
.wfn, .wfx, and .fchk for Gaussian wavefunction files. fchk files can also be written by psi4.
.molden for molecular wavefunction files (molden-style format). At present, psi4, orca, and molden are supported. Both STO and GTO molecular wavefunctions are supported. The molden.input extension works the same as a .molden file, but forces the use of the orca dialect (orca generates .molden.input files by default).
.pwc for the electron density plus Kohn-Sham state information from Quantum ESPRESSO, generated using the pw2critic.x program in PP/. Optionally, one or two checkpoint files (.chk) from wannier90 are read.

Specific details about each of these field types are given below. In some cases, the extension of the file to be loaded may be different from the one that critic2 expects. For instance, WIEN2k creates files with extensions other than clmsum for different fields (e.g. .clmup for the spin-up density) but the same format. In those cases, the keywords WIEN/ELK/… can be used to force critic2 to read a file with a specific format.

Cube Files (cube, bincube)

Cube files can be generated by Gaussian’s cubegen utility, but also by other programs such as Quantum ESPRESSO. If the cube file corresponds to a crystal structure, it is possible that it contains a non-orthogonal grid (despite its name). Binary cube files (.bincube) are created by critic2.

Abinit Density Grids (DEN)

Density files generated by abinit come in a binary DEN format (file_DEN or file_PAWDEN) and can be loaded by simply passing them to LOAD. These files contain a grid of field values. Alternatively, other files can be written using the “prt*” options in abinit, containing different scalar fields in the same format, including the electron localization function (prtelf, _ELF suffix), the Laplacian of the electron density (prtlden, _LDEN), the exchange-correlation potentiatl (prtvxc, _VXC), the kinetic energy density (usekden, prtkden, _KDEN), the Hartree potential (prtvha, _VHA), the Hartree plus xc potential (prtvhxc, _VHXC), and the total potential (prtpot, _POT).

VASP Files (CHGCAR, AECCAR*, CHG, ELFCAR)

VASP fields can come in two varieties. The CHGCAR, CHG, and AECCAR* files give grid values in a higher precision and multiplied by the cell volume. These files can be read directly with LOAD. The other format is the ELFCAR (containing the values of the ELF function). In this case, the grid values are not be multiplied by the cell volume. Both can be loaded with:

LOAD CHGCAR
LOAD ELFCAR

If you have a file in this format that does not conform to those names, you can use either the VASP (CHGCAR, AECCAR*) or the VASPNOV (ELFCAR) keyword to force using one of the two formats:

LOAD VASP STRANGE_CHG_FILE_NAME.elfcar
LOAD VASPNOV STRANGE_CHG_FILE_NAME.chgcar

The CHGCAR files are relatively safe to use, but CHGs data points may or may not be multiplied by the unit cell volume, and if you use the latter you may have to tell critic2 which one it is. If you have a CHG file that is multiplied by the volume, then you can force critic2 to load it in the same way as a CHGCAR with:

LOAD VASP CHG

I have never seen an ELFCAR that is multiplied by the volume, which would not make much sense anyway. So the default behavior for ELFCAR files is to assume that it is not.

In spin-polarized and non-collinear systems, the CHGCAR and CHG contain more than one grid. In spin-polarized calculations, the first grid is the total electron density, which is follwed by the spin density. In non-collinear calculations, the first grid is the total electron density followed by the spin magnetizations in the x-, y-, and z- directions. You can select which of these grids to read by appending an integer after the file name (iblock.i). Optionally, use the keyword RHO (iblock.i = 1), SPIN (iblock.i = 2), MAGX (iblock.i = 2), MAGY (iblock.i = 3), or MAGZ (iblock.i = 4). Similarly, the ELFCAR in a spin-polarized calculation contains two grids: the spin-up ELF (first block) and the spin-down ELF (second block). If no extra integer or keyword is passed, then the first grid in the file is read.

Note that the pseudopotential charges (ZPSP) are NOT read from the POTCAR.

CASTEP Files (fmt)

CASTEP density, potential, and ELF files ending in _fmt can be read by critic2. They contain a grid of values. Important: the cube and xsf files generated by CASTEP are not compatible with critic2 because they repeat the first point as last. Read the _fmt files directly instead.

QE Wavefunction (pwc) and Wannier90 Checkpoint Files (chk)

The .pwc file format is generated by the pw2critic utility in the post-processing (PP) package of Quantum ESPRESSO. The pwc file contains the crystal structure, k-point and plane-wave information, the Kohn-Sham state coefficients, as well as the various mappings necessary to use them. A field loaded using a .pwc file corresponds to the pseudo-valence density of the calculation, but the individual Kohn-Sham states are read as well, which allows additional calculations in critic2, such as the extraction of individual Bloch functions.

The LOAD command for the pwc file can be supplemented with a second file (with extension .chk) generated by Wannier90 during the calculation of the maximally localized Wannier functions of the system. A field loaded in this way can be used in the calculation of delocalization indices in periodic systems ¹. An additional file also generated by wannier90 is necessary to provide the MLWF orbital rotations in spin-polarized systems. In that case, the first checkpoint file comes from the up-spin wannier90 calculation and the second checkpoint file from the down-spin wannier90 calculation. See the FeO case in the DI example.

Instead of the total electron density, the contribution to the density from specific states can be selected using the SPIN, KPT, BAND, and ERANGE keywords. The SPIN keywords selects the total electron density (0), the spin-up density (1), or the spin-down density (2). This only works for a spin-polarized calculation. The KPT and BAND keywords allow selecting specific band and k-point contributions from the list. If these keywords are not present, all k-points and bands are used. The ERANGE keyword can be used to select only the states whose band energy is between emin.r and emax.r. Both energies must be given in electronvolts.

WIEN2k (clmsum, etc.)

In the case of WIEN2k, not all the necessary information is encapsulated in the .clmsum (for instance, the muffin tin radii are missing), so it is necessary to provide a second file: the .struct file.

If the file extension is not .clmsum, .clmup, or .clmdn, you can force critic2 to use the WIEN2k interpreter by preceding the file name with the keyword WIEN:

LOAD WIEN file.vtot

Some of these files (e.g. potentials) may have different normalizations for the l=0 m=0 components in the muffins if they do not represent electron densities (or spin densities). See RHONORM and VNORM in the additional options.

Elk (STATE.OUT and Other OUT Files)

There are three ways of loading a scalar field generated by elk. If a single .OUT file is passed to critic2, then it is assumed to be a grid generated by plot3d (tasks 33, 43, 53, and 63). The loaded scalar field will be a three-dimensional grid of values (i.e. the same as a cube file).

If two .OUT files are passed to critic2, then it is assumed that the first is the STATE.OUT and the second is the GEOMETRY.OUT for the current elk calculation, and the scalar field will be calculated analytically by evaluating the augmented plane-waves, same as in WIEN2k. Similar to WIEN2k, elk’s STATE.OUT is missing the structural information required to calculate the electron density, which is why the GEOMETRY.OUT is necessary.

Elk’s STATE.OUT is version-dependent. The following versions of elk have been tested and work with critic2: 5.2.14, 5.2.10, 4.3.6, 4.0.15, 3.3.17, 3.3.15, 3.1.12, 3.0.18, 3.0.4, 2.3.22, 2.3.16, 2.2.10, 2.2.9, 2.2.8, 2.2.5, 2.2.1, 2.1.25, 2.1.22, 1.4.22, 1.4.18, 1.4.5, 1.3.31, 1.3.30, 1.3.24, 1.3.22, 1.3.20, 1.3.15, 1.3.2, 1.2.20, 1.2.15, 1.1.4, 1.0.17, 1.0.16, and 1.0.0. If your version is not supported (most likely because it is newer; I tend to check this part of the code every aeon or so), please contact me.

Using a patched version of elk it is possible to generate other files (OTHER.OUT) containing the spherical harmonics/plane waves representation of other fields, such as the elf and the Coulomb potential. See the tools/elk_mod directory for the patch source files. When the OTHER.OUT (e.g. ELF.OUT) file is passed as the third argument to the LOAD command, then that field is loaded instead of the density.

aiPI ion Files (ion)

For aiPI inputs, several .ion files are necessary. The .ion files can be associated to atoms either by using their non-equivalent atom number (nat1.i) or the same atomic name as used in the crystal structure specification (at1.s).

DFTB+ Files (xml File and Others)

To load a field from the calcualted wavefunction using DFTB+, three files are required: the detailed.xml file (written by the WriteDetailedXML option), the eigenvec.bin file (using the WriteEigenVectors options), and the file containing the wavefunction coefficients, which has been provided by the authors of the DFTB+ package and can be obtained online The DFTB type is detected by the .xml extension of the first file and, therefore, the name of these files need not be detailed.xml or eigenvec.bin.

Two different types of systems can be run in DFTB+: crystals under periodic boundary conditions, and molecules in the gas-phase. In the first case, if using k-points other than Gamma, the wavefunction is complex, and the evaluation of the density is relatively slow. Molecules use real wavefunctions and they tend to be less crowded than crystals, so using molecular DFTB fields should be considerably faster.

See the DFTB+ example for worked out cases.

Gaussian-style Wavefunction Files (wfn/wfx)

Gaussian-style wavefunction files can be loaded by passing a .wfn or .wfx file to critic2. This option is normally used in combination with MOLECULE on the same file in order to obtain the same structure from the isolated molcule. The evaluation of wfn/wfx densities and their derivatives is analytical (no grids involved). If your file does not have the wfn/wfx extension, you can force this format preceding the file name with the WFN/WFX keyword.

A number of field modifiers can be used with molecular wavefunction fields to select derivatives of the electron density, a derivative field, or a particular molecular orbital. For instance, if a wfn/wfx file is loaded in field 1, then $1:4 gives the value of molecular orbital number 4. Other labels such as HOMO and LUMO are also available.

Supported wavefunction types: restricted (RHF), unrestricted (UHF), and fractional occupation (MP2, etc.). Note that Gaussian-style wfn and wfx files do not contain information about the virtual orbitals.

For molecular calculations that use effective core potentials (ECPs), Gaussian provides the core density as a sum over atom-centered Gaussian functions (“electron density functions” or EDFs). These functions are only present in a usable form in a wfx file, but not in the wfn format. Hence, densities for molecules containing ECPs will be correct and contain the core contribution only if the wavefunction is provided to critic2 using the wfx format. If a wfx file is not available, core augmentation using critic2’s internal density tables is an option.

The development notes contain a list of programs and versions whose wfn and wfx files are known to work with critic2, and some particularities regarding these formats. Orca wfn and wfx files can be read by critic2.

Gaussian Formatted Checkpoint Files (fchk)

Similar to wfn/wfx, formatted checkpoint files (.fchk) contain the structure and molecular wavefunction information. They are generated from a Gaussian checkpoint file (.chk) using the formchk utility but can also be written by psi4. They should be used preferably in combination with a preceding MOLECULE command on the same file. The calculation of the density and related properties is analytical. Field modifiers can be used to select particular molecular orbitals, same as in the wfn/wfx case. If your file does not have a .fchk extension, you can force this format by preceding the file name with the FCHK keyword.

Supported wavefunction types: restricted (RHF) and unrestricted (UHF).

Some parts of critic2 (the MOLCALC keyword and orbital selection via field modifiers in arithmetic expressions) may use virtual orbitals if available. By default critic2 will not read orbitals with zero occupation. The virtual orbitals can be read by passing the READVIRTUAL keyword to LOAD. This keyword only makes sense in formats that contain virtual orbital information (fchk and molden).

The development notes contain a list of programs and versions whose fchk files are known to work with critic2, and some particularities regarding these formats.

Molecular Molden-Style Files (molden and molden.input)

Similar to wfn/wfx, these files contain the molecular wavefunction information. They should be used preferably in combination with a preceding MOLECULE command on the same file. The calculation of the density and related properties is analytical (i.e. no numerical derivatives). Field modifiers can be used to select particular molecular orbitals, same as in the wfn/wfx case. If your file does not have the molden extension, you can force this format by preceding the file name with the MOLDEN keyword for the generic molden format or MOLDEN_ORCA or MOLDEN_PSI4 for the particular dialects those programs use.

Molden files generated by different programs use slightly different (and incompatible) versions of the same format, particularly regarding normalization of the primitives and other details. The following programs work with critic2:

Versions of psi4 from around 2016 onwards (current tests are based on the 1.4 version).
ADF wavefunctions in terms of STO basis functions (implemented and tested in 2017).
ORCA, version 4.0.2 and probably others, implemented and tested in 2020. Orca actually generates files with .molden.input extension. These are also understood by critic2, and passing them directly informs critic2 that the molden file uses the ORCA dialect.

The development notes contain a list of programs and versions whose molden files are known to work with critic2, and some particularities regarding these formats.

If no dialect is provided by the user, critic2 will try to “guess” the format of the molden file from its contents. The psi4 or orca dialects can be forced by using the ORCA or PSI4 keywords after the name of the molden file. If you try molden files from any other program, please let me know how it goes so I can update this manual.

The supported wavefunction types are restricted (RHF) and unrestricted (UHF), both with STO and GTO primitives. Virtual orbitals can be read using the READVIRTUAL keyword. Cartesian and spherical primitives are supported, but only up to g.

Making Copies of a Field

The LOAD COPY keyword can be used to make a copy of the field id.s to the next available slot, or to id2.s if it is given with the keyword TO.

Loading the Promolecular Density

The keyword PROMOLECULAR is used to load a promolecular density field (same as field number 0). If a finite fragment of the crystal is passed to LOAD PROMOLECULAR as an xyz file (file.xyz) using the optional FRAGMENT keyword, then the sum of atomic densities is built using only the atoms in that fragment.

Additional LOAD Options

The definition of a field can be supplemented by additional optional keywords that control its behavior. The additional options are passed to LOAD after the file name. The list of options follows.

NEAREST|TRILINEAR|TRISPLINE|TRICUBIC|SMOOTHRHO [NENV nenv.i] [FDMAX fdmax.r]

Choose the grid interpolation mode in a grid field. NEAREST, use the field at the nearest grid point (zero first and second derivatives). TRILINEAR, trilinear interpolation (zero second derivatives). TRISPLINE, 3d-spline interpolation (adapted from the abinit code, code by A. Lherbier). If some derivatives are not available (first and second in NEAREST, second in TRISPLINE), they are taken as zero. TRICUBIC, tricubic interpolation using local information, see Lekien and Marsden, Int. J. Numer. Meth. Eng., 63 (2005) 455-471. SMOOTHRHO, polyharmonic spline interpolation with smoothing of electron density, see Otero-de-la-Roza, J. Chem. Phys. 156, 224116 (2022).

By default, TRICUBIC is used. TRISPLINE may require a lot of memory in very large cube files (say, 400^3 or larger). However, the memory allocation for the TRISPLINE interpolation only happens if the parts of critic2 that require derivatives or interpolation at points outside of the grid are called. If only grid-based algorithms (e.g. YT, BADER, NCIPLOT with default NSTEP,…) are used then no additional memory is used.

TRICUBIC and TRISPLINE are sensitive to rapid variations in the field and numerical noise, which may result in small oscillations in the field value.

If your field is an (all-electron) electron density, then the best interpolation method by far is SMOOTHRHO, and you should use it if this is your case. SMOOTHRHO eliminates the noise associated with the interpolation near the nuclei and allows finding the critical points of the electron density (see Otero-de-la-Roza, J. Chem. Phys. 156, 224116 (2022)). There are two options to SMOOTHRHO. NENV sets the number of grid nodes in the stencil for the polyharmonic spline interpolation (default: 512). FDMAX sets the distance for the interpolant combination function that ensures the continuity of the global interpolant. By default, fdmax.r is 1.75, which means that all grid nodes within 1.75 times the longest Voronoi-relevant vector are used in the combination. The higher NENV and FDMAX, the smoother the interpolation will be, but also more costly. For examples of how to find critical points of the electron density when it is given on a grid, see the corresponding example.

(Applies to: grids. Default: TRICUBIC.)

NOCORE

Deactivate the core contribution to this field (see ZPSP below for more details).

EXACT|APPROXIMATE

The calculation of the electron density in aiPI and DFTB+ fields is relatively expensive. Using the APPROXIMATE keyword, the atomic contributions to the density are calculated by interpolating from a radial grid, that is precomputed at the beginning of the run. EXACT calculates the fully analytical values of the field with no approximations.

Applies to: PI and DFTB+ fields. Default: APPROXIMATE.

RHONORM|VNORM

In WIEN2k, the clmsum and other files representing the electron density come in two different normalizations: The l=0 m=0 radial component may or may not divided by $\sqrt{4\pi}$ . RHONORM = use the same normalization as the density (divided by $\sqrt{4\pi}$ ). VNORM = use the same normalization as the potentials (no division).

Applies to: WIEN2k. Default: RHONORM.

NOTESTMT

Augmented plane-wave scalar fields (WIEN2k and elk) are checked by default for discontinuities on the muffin tin surface, which can affect the location of critical points and the tracing of gradient paths through the muffin tin surface. This test may take some time if the number of atoms in the system is large. NOTESTMT skips the discontinuity test.

Applies to: WIEN2k and elk. Default: test the muffin tin discontinuity.

ANALYTICAL/NUMERICAL

NUMERICAL = Calculate the derivatives of the field numerically. ANALYTICAL = calculate the derivatives of the field analytically. Mostly used for testing the implementation of new field types.

Applies to: all. Default: ANALYTICAL.

TYPNUC {-3,3}

Controls whether the nuclei are maxima (-3) or minima (3) of the field.

Applies to: all. Default: -3 (maxima).

NORMALIZE n.r

Normalize the grid integral over the unit cell to n.r. In the case of an electron density, n.r is the number of electrons per cell.

Applies to: grids.

{NAME|ID} id.s

Name the field using a string identifier. The identifier can be used in place of the field number in any expression (e.g. $name instead of $1). The label for this field will written in the part of the output describing the field. The string should be a single word containing letters and numbers. You can list all available field identifiers using the LIST keyword.

ZPSP at1.s q1.r...

The electron density from a pseudopotential/plane-waves calculation, given on a grid, only represents valence electrons. In order to get an approximation to the all-electron density, it can be augmented by summing the corresponding core contributions to the electron density at the atomic sites. The ZPSP keyword is used to determine how many electrons are added by each nucleus. For atom at1.s, the value passed to ZPSP (q1.r) must correspond to the number of valence electrons for that atom. The use of ZPSP in LOAD activates the calculation of the core contribution, which is added to the field’s value. Most of the time, core augmentation is not entirely satisfactory, and is best avoided if possible (e.g. by writing the all-electron density using the PAW transformation, if PAW is being used).

Applies to: grids, although it can be used with any other field. Default: no core augmentation.

Field Arithmetics

New fields can be defined as combinations of the existing ones using the LOAD AS keyword. This allows the creation of fields that are not directly computed by the electronic structure program, and is very useful in combination with the CUBE keyword to calculate and visualize new grids. For instance, if rhoup.cube and rhodn.cube are the spin-up and spin-down densities, it is possible to get the total density and the spin density using:

LOAD rhoup.cube
LOAD rhodn.cube
LOAD AS "$1+$2"
LOAD AS "$1-$2"

The spin-up density is loaded as field 1, rhodn.cube is field 2, field 3 is the sum of both and field 4 is the spin density, which can, for instance, be graphically represented or integrated in the atomic basins to get the atomic magnetic moments (see the INTEGRABLE keyword).

The list of arithmetic operations that can be performed on fields is quite large. If no additional keywords are present between LOAD AS and the arithmetic expression (like FFT, CLM, or RESAMPLE, see below), then the resulting new field can be of two types. If the arithmetic expression involves at least one grid, the resulting field will be a grid with the same number of points. If more than one grid appears in the expression, then the new grid size is the maximum of the number of points in every dimension.

If the expression contains no grid fields, then the resulting field is a “ghost field”. A ghost field is just an arithmetic expression that is parsed and processed every time the field is evaluated. The analyitcal derivatives for a ghost field are not available, so the NUMERICAL option is the default for this kind of field.

Regardless of the types of fields in the expression passed to LOAD AS, a grid field can be enforced by explicitly specifying the size of the grid, either by giving the number of points in each direction (n1.i, n2.i, n3.i) or by adopting the size of another grid (SIZEOF, the grid used as mold is field id.s). If grid fields appear in the expression but you still want a ghost field, then use the GHOST keyword.

The PROMOLECULAR option to LOAD AS allows the creation of a grid out of the promolecular density, with n1.i, n2.i, n3.i points in each direction. Alternatively, the size of grid field id.s can be used with the SIZEOF keyword. Similarly, CORE creates a grid using only the core densities, as specified by the number of valence electrons (ZPSP) of the atoms, with the same syntax as the ZPSP option to LOAD (q1.r is the number of valence electrons of atom at1.s). In this case, the use of ZPSP is mandatory. If a fragment of the crystal is passed as an xyz file to any of those keywords (with the optional FRAGMENT keyword), then only the atoms in the fragment contribute to the sum of atomic (or core) densities.

The FFT, CLM, and RESAMPLE keywords of LOAD AS apply only to specific types of fields:

FFT applies only to grid fields (except FFT LAP, see next bullet point). The FFT creates a new grid by transforming the original grid using a fast Fourier transform plus some operation. The operation is determined by the next keyword:
- GX, GY, GZ: the new grid is the first derivative with respect to the corresponding Cartesian coordinate.
- GMOD: the new grid contains the gradient norm of the original grid.
- HXX, HXY, HXZ, HYY, HYZ, HZZ: the new grid is the second derivative with respect to the indicated coordinates.
- LAP: the new grid contains the Laplacian of the original grid.
- POT: if the original grid is a density in atomic units (electrons/bohr^3), POT makes a new grid containing the electrostatic potential ( $V$ ) associated with that density ( $\rho$ ), calculated by solving the Poisson equation ( $\nabla^2V = -4\pi\rho$ ), and also in atomic units (Hartree/electron).
The FFT LAP keyword can also be applied to WIEN2K and elk fields. In both cases, the keyword creates a field of the same type containing the Laplacian of the field id.s. The Laplacian is calculated analytically both inside the muffin tins and in the interstitial using the appropriate transformation.
CLM defines a new field as the addition (id1.s+id2.s) or substraction (id1.s-id2.s) of fields of the WIEN2K or elk type by using the ADD and SUB keywords respectively. Note that the muffin tin radii, number of plane-waves, etc. have to be the same for both source fields.
RESAMPLE is used to create larger grids from smaller grids via Fourier transform. The RESAMPLE keyword requires a source grid field (id.s) and the size of the new grid (the number of points in each direction, n1.i, n2.i, and n3.i). Note that for non-smooth source grids this operation may induce oscillations in the regions with low field value.

Changing the Field Options After LOAD (SETFIELD)

The options of a given field can be changed anywhere in the input after it has been loaded using the SETFIELD keyword:

SETFIELD {id.s} [NEAREST|TRILINEAR|TRISPLINE|TRICUBIC|SMOOTHRHO [NENV nenv.i] [FDMAX fdmax.r]]
   [EXACT|APPROXIMATE] [RHONORM|VNORM] [NOCORE]
   [NUMERICAL|ANALYTICAL] [TYPNUC {-3,-1,1,3}]
   [NORMALIZE n.r] [ZPSP at1.s q1.r...]

SETFIELD changes the properties of field with identifier id.s. The keywords have the same meaning as above.

Unloading a Field (UNLOAD)

UNLOAD {id.s|ALL}

The UNLOAD keyword unloads the field id.s or all fields (ALL).

The Reference Field (REFERENCE)

One of the loaded fields is chosen as the reference field. The reference field is used as the primary field in all computations. For instance, when BASINPLOT is used, the basins of the reference field are plotted. If one uses an integration method (INTEGRALS, QTREE, YT, BADER), the reference field provides the attraction basins in which other scalar fields are integrated. If AUTO is used, the critical points of the reference field are found. And so on.

The first field loaded becomes the reference field. If no fields have been loaded, then the reference is the promolecular density (field 0 or rho0). In order to change the reference field, the REFERENCE keyword can be used:

REFERENCE id.s

REFERENCE sets field id.s as reference.

Typical Output for the LOAD Keyword

The output when loading a field is composed of two parts: a first part containing information particular to the type of field loaded, and a second part containing the field flags.

For instance, loading a cube file (or any grid) gives an output whose first part looks like this:

+ Field number 1
  Source: urea.rho.cube
  Type: grid
  Grid dimensions : 81  81  72
  First elements... 1.084595807000E-03  1.081705980000E-03  1.085704517000E-03
  Last elements... 1.121839312000E-03  1.114200414000E-03  1.106977600000E-03
  Sum of elements... 2.316327381350E+04
  Sum of squares of elements... 6.142727758913E+03
  Cell integral (grid SUM) = 47.99979210
  Min: 3.81287906E-04
  Average: 4.90340095E-02
  Max: 1.01178629E+00
  Interpolation mode (1=nearest,2=linear,3=spline,4=tricubic): 4

The output for the other grid file formats (abinit DEN, VASP, qub, etc.) is essentially the same. This output gives the field integer identifier (1 in this case), the source file, the type of field (grid), and a number of information specific to grid fields. This includes the number of points on the grid (grid dimensions), the first and last three elements, and the sum and sum of the squares of the elements. The Cell integral line is particularly important because it gives the integral of the grid field over the unit cell (i.e. $\frac{(\sum_i f_i) V}{N}$ with $V$ the unit cell volume and $N$ the total number of grid points). If the grid field represents the electron density, it will approximately sum to the number of electrons in the cell (in the example above, 48). The output also provides the minimum, maximum, and average field value on the grid, and the chosen interpolation mode.

The second part of the output is common to all fields:

  Use core densities? F
  Numerical derivatives? F
  Nuclear CP signature: -3
  Number of non-equivalent critical points: 5
  Number of critical points in the unit cell: 16
  Alias for this field (1): $1

It contains the field flags: whether core densities are being used, if derivatives are calculated numerically or analytically, and the signature of the nuclear critical points. In addition, it also gives the number of critical points found for this particular field. Since we just loaded it, the number of critical points is the same as the number of atoms (which are always considered critical points). Finally, the output lists the aliases that can be used in arithmetic expressions to refer to this field.

For WIEN2k and elk fields, the first part of the output contains:

+ Field number 2
  Source: srtio3.clmsum
  Type: wien2k
  Complex?: F
  Spherical harmonics expansion LMmax: 6
  Max. points in radial grid: 781
  Total number of plane waves (new/orig): 350/3119
  Density-style normalization? T

The APW-specific information includes, whether the field contains complex coefficients in the expansions (it does if the crystal structure is not centrosymmetric), the maximum L for the spherical harmonics expansion inside the muffin tins, the number of radial points inside the muffin tins, the cutoff for discarding negligible plane-wave contributions in the interstitial, and the number of plane-waves remaining after pruning.

In the case of WIEN2k and elk fields, the muffin tins are tested for discontinuities:

* Muffin-tin discontinuity test
  Atom: 1 rmt= 1.8000000 RMS/max/min(fout-fin) = 0.000420  0.001878  -0.000006
  Atom: 2 rmt= 1.8000000 RMS/max/min(fout-fin) = 0.000116  0.000473  0.000006
  Atom: 3 rmt= 1.8000000 RMS/max/min(fout-fin) = 0.000314  0.002103  -0.000059
  Atom: 4 rmt= 1.8000000 RMS/max/min(fout-fin) = 0.000087  0.000410  -0.000006
  Atom: 5 rmt= 1.8000000 RMS/max/min(fout-fin) = 0.000738  0.005342  -0.000313
+ Assert - no spurious CPs on the muffin tin surface: T

This system has five muffin tins, all of them with 1.8 bohr radius. The test determines whether the discontinuity at the muffin flips the sign of the radial gradient of the field, which would create spurious critical points. In this case, the field passes the test, which means it can be used without problem in AUTO or GRDVEC. The NOTESTMT field option can be used to skip this test.

A typical LOAD output for molecular wavefunctions is:

+ Field number 3
  Source: benzene.wfx
  Type: molecular wavefunction
  Wavefunction type: restricted
  Number of MOs (total): 21
  Number of MOs (occupied): 21
  Number of primitives (GTO): 216
  Number of EDFs: 0
  Basis set data NOT available

It contains the type of wavefunction (restricted closed-shell in this case), the number of molecular orbitals, number of primitives, and number of EDFs (in case the calculation used ECPs). In fields loaded from a fchk and if critic2 is compiled with the libcint library, the basis set data will be avilable for molecular integrals with the MOLCALC keyword.

Other field types contain specific information for their types, which is usually self-explanatory. If the LOAD keyword is the first in the run, then the field becomes reference, which will be shown in the output:

* Field number 1 is now REFERENCE.

See the REFERENCE keyword for more information.

Typical Output for the REFERENCE Keyword

The output of the REFERENCE keyword is relatively simple. First, the field set as reference is given:

* Field number 2 is now REFERENCE.

Then, the output shows the list of the integrable properties.

* List of integrable properties (3)
#  Id  Type  Field  Name
    1   v        0  Volume
    2  fval      2  Pop
    3  lval      2  Lap

This is the list of quantities that will be integrated by the basin integration methods like INTEGRALS or YT. Chaning the reference field resets the integrable properties list to its default value, which is to integrate the basin volumes, the reference field, and the Laplacian of the reference field. If the reference field is an electron density, then the last two correspond to the atomic electron population and the electron density Laplacian. The list of integrable properties can be modified with the INTEGRABLE keyword.

Finally, critic2 lists the known critical points for this field:

* Summary of critical points for the reference field
  Non-equivalent critical points = 16
  Cell critical points = 52
  Topological class (n|b|r|c): 5(16) 6(22) 3(10) 2(4)
  Morse sum: 0

If a critical point search has not been conducted, then the summary will include only the atoms as critical points.

http://dx.doi.org/10.1021/acs.jctc.8b00549 ↩