Miscellaneous tools
Molecular Calculations (MOLCALC)
The MOLCALC keyword allows the (limited) calculation of molecular properties:
MOLCALC
MOLCALC expr.s
MOLCALC PEACH
mo1a [>] mo1r k1
mo2a [>] mo2r k2
[...]
ENDMOLCALC/END
MOLCALC HF
MOLCALC ... [ASSIGN var.s]
MOLCALC integrates a field or an expression in a Beckestyle mesh. The MOLCALC keyword works mostly with molecular wavefunction (wfn/wfx/fchk/molden) fields but crystals can also be used with some options.
Without any other keywords, MOLCALC integrates the reference field. If
an expression (expr.s
) is used, the scalar field generated by this
expression is integrated on the mesh.
The PEACH keyword calculates the index for the characterization of
electronic excitations based on orbital overlap defined as uppercase
lambda in Peach et al., J. Chem. Phys. 128 (2008) 044118.
Small values (close to zero) correspond to little orbital overlap and
Rydberg excitations. Large values (close to one) are caused by high
orbital overlap and correspond to local excitations. Charge transfer
excitations can have both high or low PEACH value. Each line gives the
initial (mo1a
) and final (mo1r
) molecular orbital and the
oscillator strength (k1
) for each component of the excitation. If
Gaussian is used, copying and pasting the output of the TD
keyword
is recommended. Using the PEACH keyword requires the virtual orbitals,
which in turn requires using a fchk or molden file format with the
READVIRTUAL keyword (see the LOAD
keyword).
MOLCALC HF calculates the HartreeFock energy from the basis set
information in the reference field using molecular integrals, which
are in turn calculated by the libcint library.
Using this keyword requires compiling with critic2 with libcint. At
this moment, this keyword serves mostly for testing, debug, and
development purposes, as it is very inefficient (and not very
useful). The routine called by this keyword shows how to use libcint
inside critic2. Only molecular wavefunctions can be used with
MOLCALC HF
and only file formats that provide basis set information
(right now, only Gaussian fchk, but moldenstyle files could be
implemented). Gaussian wfn/wfx do not provide basis set shells, only
primitives, so they cannot be used with MOLCALC HF.
The ASSIGN keyword can be used to save the result of the computation
to a variable (var.s
).
Characterization of SigmaHoles in Molecules (SIGMAHOLE)
The SIGMAHOLE keyword calculates the properties of a \(\sigma\)hole in a molecule. The \(\sigma\)hole is the region of positive electrostatic potential on an electronegative atom, opposite a \(\sigma\)bond. Typically, but not always, the electronegative atom is a halogen, in which case the molecule can engage in a halogen bond. The syntax of the SIGMAHOLE keyword is:
SIGMAHOLE ib.i ix.i [NPTS nu.i nv.i] [ISOVAL rho.r] [MAXANG ang.r]
The calculated properties are those described in Kolar and Hobza, Chem. Rev. 116 (2016) 5155 (see Figure 7). Namely:

dist
: The distance between the electronegative atom (X) and the \(\sigma\)hole. 
\(l_{\sigma}\) (linearity): the angle between the \(\sigma\)hole, the electronegative atom (X), and the atom to which X is bonded (the basal atom, B).

\(m_{\sigma}\) (magnitude): the value of the electrostatic potential at the \(\sigma\)hole.

\(r_{\sigma}\) (range): the distance between the electronegative atom and the point at which the electrostatic potential changes sign, in the direction of the \(\sigma\)hole. Only defined in case of a \(\sigma\)hole with positive electrostatic potential.

\(s_{\sigma}\) (size): the are of the density isosurface associated with the \(\sigma\)hole that has positive electrostatic potential. Only defined in case of a \(\sigma\)hole with positive electrostatic potential.
The SIGMAHOLE keyword can only be used in molecular systems. The calculation of the density and the electrostatic potential is carried out using the reference field, which should be a molecular wavefunction with basis set information. In addition, the calculation of the electrostatic potential requires compiling critic2 with the LIBCINT library.
The SIGMAHOLE keyword must be followed by two integers, which identify
the basal (ib.i
, B) and the electronegative (ix.i
, X) atoms in the
molecule. A number of rays starting at the electronegative atom are
traced, distributed in such a way that each occupies the same amount
of solid angle. The rays subtend an angle with the BX line that
ranges between 180 degrees and 180 minus ang.r
degrees, where
ang.r
is the option to MAXANG. By default, its value is 90 degrees
(i.e. a hemisphere). The number of rays traced in the u (azimuthal)
and v (polar) parameters is controlled by the NPTS keyword. These
should be high enough that the calculation of the relevant properties
is essentially converged (default: 101 each).
The point is found along each ray where the electron density is equal
to the rho.r
isovalue. This value can be controlled with the ISOVAL
keyword; it is 0.001 a.u. by default. The electrostatic potential is
calculated at this point, and a twodimensional mapping of the
electrostatic potential on the density contour is assembled. From
this mapping, the properties above are calculated. Note that more than
one \(\sigma\)hole can be found. In the case of multiple
\(\sigma\)holes, if the positivepotential regions are in contact,
the calculation of \(s_{\sigma}\) is not done. Also, it is important
that the \(\sigma\)hole does not exceed the maximum angle (MAXANG)
and that NPTS is sufficiently high to have an accurate \(s_{\sigma}\)
value.
In addition to the list of \(\sigma\)hole properties, SIGMAHOLE also
writes two files for plotting the electrostatic potential mapped on
the density isosurface. The root_sigmahole.dat
file (root
is the
default prefix of the run) contains
the isosurface data. Each row is one traced ray. The columns are:
values of the u and v parameters, \(\theta\) and \(\phi\) spherical
coordinates, the Cartesian coordinates of the point (in the local
Cartesian frame), the distance of the isosurface to the
electronegative atom, and the value of the electrostatic potential. A
template gnuplot script, root_sigmahole.gnu
, is generated to aid in
the visualization of the data file. It is recommended that this plot
is always checked to make sure that the calculated properties in the
table make sense.
The ExchangeHole Dipole Moment Dispersion Model (XDM)
The XDM keyword calculates the dispersion energy using the exchangehole dipole moment (XDM) model. See J. Chem. Phys. 127, 154108 (2007), J. Chem. Phys. 136, 174109 (2012), and J. Chem. Phys. 138, 204109 (2013) for more details. The syntax of the XDM keyword is:
XDM GRID [RHO irho.s] [TAU itau.s] [ELF ielf.s]
[PDENS ipdens.s] [CORE icor.s] [LAP ilap.s]
[GRAD igrad.s] [RHOAE irhoae.s] [XB ib.s]
[XA1 a1.r] [XA2 a2.r] [ONLYC] [UPTO {6810}]
XDM [QEPOSTG] [FILE file.s] [BETWEEN at1.i at1.i ... AND at1.i at2.i ...]
[NOC6] [NOC8] [NOC10] [SCALC6 s6.r] [SCALC8 s8.r] [SCALC10 s10.r]
[C9] [SCALC9 s9.r]
[DAMP a1.r a2.r] [DAMP3 a3.r a4.r] [DAMP3BJN 36sqrt6]
XDM a1.r a2.r chf.s
There are three modes of operation for the XDM keyword.
In QE and POSTG, the volumes, moments, and coefficients are read from a Quantum ESPRESSO output file (QE) or a postg output file (POSTG). The output file is the one given via the FILE keyword or the loaded using CRYSTAL or MOLECULE, if no FILE is given. Using the dispersion coefficients read from the file, the XDM energy is recalculated. If BETWEEN and AND are given only the dispersion interaction between those pairs of atoms is calculated. This keyword is used mostly for testing purposes. A number of additional options allow controlling the expression of the XDM energy. If NOC6, NOC8, or NOC10 are given, omit the corresponding contribution to the pariwise dispersion energy. If C9 is given, calculate the threebody dispersion term using BJ damping. If SCALC6, SCALC8, SCALC9, or SCALC10 are given, scale the corresponding term by the real number directly following the keyword. If DAMP is given, recalculate the pairwise dispersion energy with the corresponding damping function parameters. If DAMP3 is given, use these parameters for the 3body damping. Use DAMP3BJN to choose the 3body damping function: a product of three BJ3 (3), three BJ6 (6), or the square of three BJ6 (sqrt6) damping functions.
In the GRID mode, the information necessary to calculate the XDM dispersion energy and related quantities is provided using grid fields.
If neither QE/POSTG nor GRID are given, a molecular XDM calculation is
assumed, with damping function coefficients a1.r
and a2.r
and
functional selector chf.s
. The latter can be either a keyword for a
functional (blyp
, b3lyp
, etc.) or a number between 0 and 1
indicating the fraction of exact exchange. The molecular XDM code can
be used only with wavefunction fields and DFTB fields.
The rest of the information in this section applies to the XDM GRID keyword. XDM GRID uses the electron density (RHO), the kinetic energy density (TAU), the Laplacian (LAP), and the gradient of the electron density (GRAD) to compute the exchangehole dipole moment in the BeckeRoussel model (B). The promolecular density (PDENS) and the core density (CORE) are used to calculate a Hirshfeld partitioning of the unit cell. All of these fields must be available or are calculated when running XDM. The corresponding keywords accept an integer, corresponding to a previously loaded field. During the XDM run, cubes for all of these fields are generated so they can be loaded in subsequent runs.
The list of requirements is:

RHO: the electron density. By default,
irho.s
is the reference field. This field is required in XDM. It is also required to give the pseudopotential charges using the ZPSP keyword for all types of atoms in the system. 
TAU: the kinetic energy density. It is used in the calculation of B, and can be extracted from the ELF. Hence, it is required except if the ELF or the B is given. If the ELF is used instead of TAU, a cube file (
tau.cube
) is written. 
ELF: the electron localization function. Can be used in place of TAU. This is useful because most programs (e.g. QE, VASP) generate cubes for the ELF but not for the kinetic energy density.

PDENS: the promolecular density. It is generated by critic2 if not present in the XDM call, and written to a cube file (
pdens.cube
) for future use. 
CORE: the core density. It is generated by critic2 if not present in the XDM call and written to a cube (
core.cube
), unless the B and RHOAE are given, in which case it is ignored. The ZPSP of all atoms is required in order to calculate this quantity. 
LAP: the Laplacian of the electron density. It is generated by Fourier transform of RHO unless it is given or B is given, in which case it is not needed and its calculation is skipped.

GRAD: the gradient of the electron density. It is generated from RHO unless B is given. If B is available, the calculation of GRAD is skipped.

RHOAE: the allelectron density on a cube. If given, replaces the pseudodensity plus core in the calculation of the atomic volumes.

XB: the exchangehole dipole moment in the BeckeRoussel model. Calculated from the above unless given.
Only closedshell (nonspinpolarized) systems can be calculated for now.
The usual way of running XDM for the first time is:
CRYSTAL rho.cube
ZPSP C 4 O 6 H 1
LOAD rho.cube
LOAD elf.cube
XDM GRID elf 2
This generates several cube files: roottau.cube
, rootpdens.cube
,
rootcore.cube
, rootlap.cube
, rootgrad.cube
, and
rootb.cube
. Subsequent runs can circumvent the calculation of B,
PDENS, and CORE by doing:
CRYSTAL rho.cube
ZPSP C 4 O 6 H 1
LOAD rho.cube
LOAD rootb.cube
LOAD rootpdens.cube
LOAD rootcore.cube
XDM GRID XB 2 PDENS 3 CORE 4
Note that passing RHO 1
is not necessary because rho.cube
is the
first field loaded (hence the reference) and assumed by default to be
the density by XDM. The ZPSP of all atoms are needed for an XDM
calculation. During the calculation, cubes for almost all
the properties above are generated so they can be reused in future
calculations. The other options are:
XA1 [a1.r]
The value of the a1 damping parameter (adimensional). Default: 0.6836 (PW86PBE parametrization for QE).
XA2 [a2.r]
The value of the a1 damping parameter (in angstrom). Default: 1.5045 (PW86PBE parametrization for QE).
ONLYC
Calculate the dispersion coefficients but not the dispersion energy, forces, and stress. By default, they are calculated.
UPTO {6810}
Only calculate the contributions to the energy coming from the \(C_6\) term (6), from the \(C_6\) and \(C_8\) terms (8) and from \(C_6\), \(C_8\), and \(C_{10}\) (10). The latter is the default.
Control Commands and Options
This section lists a number of keywords that are used to control the
operation of critic2. In general, they should be given before the
keywords that employ those options are used. For instance, ODE_MODE
should be given before integrating by bisection.
Gradient Path Tracing (ODE_MODE)
ODE_MODE [METHOD {EULERHEUNBSRKCKDP}] [MAXSTEP maxstep.r]
[MAXERR maxerr.r] [GRADEPS gradeps.r]
The keyword ODE_MODE
is used to control the gradient path
integration algorithm. EULER
selects plain explicit Euler (1
evaluation per step) with a poor man’s technique for the step size
adaptation. HEUN
, Heun’s method (2 evaluations), with the same size
adaptation as Euler’s. BS
, BogackiShampine’s embedded 2(3)thorder
method with error estimation and first step as last. RKCK
,
RungeKuttaCashKarp embedded 45th order method with local
extrapolation and error estimate. DP
, DormandPrince 45th order with
local extrapolation and error estimate (7 evaluations per
step). maxstep.r
is the initial (and maximum) step size (bohr in
crystals, angstrom in molecules), gradeps.r
is the
gradient norm termination criterion for the gradient path. maxerr.r
is
the maximum error in the trajectory (in bohr). MAXERR
only affects
methods that provide an error estimate for the predicted steps: BS
,
RKCK
, and DP
.
This keyword applies to all gradient paths except those in qtree. The
defaults are BS
method with maxstep.r
= 0.3 bohr, gradeps.r
= 1e7,
and maxerr.r
= 1e5.
Gradient Path Plot Pruning (PRUNE_DISTANCE)
Plots of gradient paths are pruned so that only one point is plotted every certain length. The PRUNE_DISTANCE keyword is used to control this length:
PRUNE_DISTANCE prune.r
Only one point is plotted every prune.r
distance (units: bohr in
crystals and angstrom in molecules, default: 0.1 bohr). Gradient paths
are forced to have steps of length less than or equal to
prune.r
. Therefore, a very small PRUNE_DISTANCE
will result in
slow gradient path tracing.
Radial Integration Method (INT_RADIAL)
The INT_RADIAL
keyword chooses the type of radial integration method
used inside spheres or atomic basins:
INT_RADIAL [TYPE {GAULEGQAGSQNGQAG}] [NR nr.i] [ABSERR err.r]
[RELERR err.r] [ERRPROP prop.i] [PREC prec.r]
The TYPE
keyword selects the quadrature method:

GAULEG
: GaussLegendre. 
QAGS
: quadpack’s dqags (generalpurpose, extrapolation, globally adaptive, endpoint singularities). All Q methods are sometimes unstable for heavy atoms and big betaspheres, but this does not happen very often. 
QNG
: quadpack’s dqng (smooth integrand, nonadaptive, GaussKronrod(Patterson)) 
QAG
: quadpack’s dqag (generalpurpose, integrand examiner, globally adaptive, GaussKronrod)
The number of radial integration points, if appropriate (GAULEG
,
QAG
) is selected with NR
. If the selected method is QAG
, the
number of points may vary from nr.r
. The allowed intervals are:
715, 1021, 1531, 2041, 2551, 3061. Critic2 selects the
appropriate interval by comparing the given nr.r
to the lower
limits of these intervals.
ABSERR
is the requested absolute error for QUADPACK
integrators. RELERR
is the requested relative error for QUADPACK
integrators. ERRPROP
controls the property for which the error is
estimated. If prop.i
corresponds to one of the integrable properties
then RELERR
and ABSERR
apply only to it. The selected property
guides the adaptive integration procedure. If prop.i
does not
represent one of the integrable properties, the maximum of the
absolute value of the properties vector is used. Note: the option
where max(abs(prop)) is used is unstable. Some spheres (usually
associated to heavy atoms) may integrate to nonsense, depending on the
optimization levels of the compiler. Therefore, this option is
disabled by default.
In the case of basin integrations, PREC
controls the precision in
the determination of the interatomic surface (bohr in crystal,
angstrom in molecules).
Default: GAULEG
, nr.r
= 50, aerr.r
= 1d0, rerr.r
= 1d12. The
default errprop is the field value and the Laplacian. delta.r
= 1e5
bohr.
Integration Meshes (MESHTYPE)
The MESHTYPE keyword chooses the type and quality of the molecular integration mesh:
MESHTYPE {BECKEFRANCHINI} [SMALLNORMALGOODVERYGOODAMAZING]
These meshes are used in molecules and crystals to calculate integrals over the whole space (cell in crystals, R3 in molecules). The MESHTYPE affects XDM integration in molecules, MOLCALC, and the MESH seeding option to AUTO.

BECKE: Beckestyle integration grid. Only available for molecules.

FRANCHINI: a Beckestyle molecular mesh with weights and quality parameters adapted from Franchini et al. J. Comput. Chem. 34 (2013) 1819. Works with molecules and crystals.
The fineness of the integration grid is be chosen by using one of SMALL, NORMAL, GOOD, VERYGOOD, and AMAZING. No pruning of the mesh is done yet.
Default: FRANCHINI, GOOD.
Cube File Precision (PRECISECUBE/STANDARDCUBE)
The PRECISECUBE and STANDARDCUBE keywords control the precision of the numbers written to cube files:
PRECISECUBESTANDARDCUBE
The field values in cube files written by Gaussian are written in exponential format with six significant digits. This precision may not be enough for some applications, particularly involving energies. The default behavior of critic2 (corresponding to the PRECISECUBE keyword) is to write cubes with 14 significant digits. Some other programs that use Gaussian cube files may not like this, however, so the keyword STANDARDCUBE is provided to make critic2 write cube files in the default Gaussian format.
Bond Distance Factor (BONDFACTOR)
The BONDFACTOR keyword controls the bond detection subroutines in critic2:
BONDFACTOR bondfactor.r
Critic2 considers two atoms are covalently bonded if their distance is
less than the sum of their covalent radii times bondfactor.r
(default:
1.4). The maximum bondfactor allowed is 2.0.
Covalent and van der Waals Radii (RADII)
The keyword RADII allows changing the internal table of covalent and van der Waals radii:
RADII {COVVDW} [at1.sz1.i] rad1.r [[at2.sz2.i] rad2.r ...]
Without any additional arguments, RADII prints the list of covalent
and van der Waals radii to the output. With COV, the following
arguments modify the covalent radii. With VDW, modify the van der
Waals radii. A list of argument pairs follows, allowing the change of
several atomic radii using the same command. The first argument in
each pair can be either the atomic symbol at1.s
or the atomic number
z1.i
of the atom whose radius we want to change. Its radius is
changed to rad1.r
. (Units: by default, bohr in crystals, angstroms
in molecules, unless the
UNITS keyword is used).
Root of Files Created by Critic2 (ROOT)
The keyword ROOT:
ROOT root.s
changes the <root>
of the critic2 run. The root is used as prefix
for most of the auxiliary files written by critic2.
Grid Calculations (SUM,MAX,MIN,MEAN,COUNT)
These five keywords:
SUM [id.s]
MAX [id.s]
MIN [id.s]
MEAN [id.s]
COUNT id.s eps.r
are used to preform calculations using the field on a grid id.s. SUM calculates the sum of the grid point values (in the output, SUM) and the integral over the unit cell (in the output, INTEGRAL). The latter is calculated as the sum of all grid points times the grid volume divided by the number of points.
MAX calculates the maximum value of all points in the grid. MIN
calculates the minimum value. MEAN calculates the average value of all
grid points. COUNT counts the number of elements that are greater than
eps.r
. If no id.s
is given in the first four commands, then the
reference field is used, provided it is defined on a grid.
Benchmark Calculations on Fields (BENCHMARK)
The BENCHMARK keyword:
BENCHMARK [nn.i]
is used to determine the average speed of the evaluation of the
reference field using nn.i
points (default: 10000). Mostly for debug
and development purposes.
Run System Commands (RUN/SYSTEM)
{RUNSYSTEM} command.s
Execute a shell command.
Echo a Message (ECHO)
ECHO string.s
Write the string string.s
to the output. Useful for partitioning
long outputs.
Terminate the Run (END,EXIT)
END
EXIT
The END or EXIT keywords terminate the execution of critic2.
Expression Evaluation
expression.s
If the input is not identified with any of the reserved keywords, then evaluate the command as an arithmetic expression and print its value. Useful for simple calculations in the command line (with critic2 q).