• 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 Becke-style 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 Hartree-Fock 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 molden-style 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 Sigma-Holes in Molecules (SIGMAHOLE)

    The SIGMAHOLE keyword calculates the properties of a -hole in a molecule. The -hole is the region of positive electrostatic potential on an electronegative atom, opposite a -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 -hole.

    • (linearity): the angle between the -hole, the electronegative atom (X), and the atom to which X is bonded (the basal atom, B).

    • (magnitude): the value of the electrostatic potential at the -hole.

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

    • (size): the are of the density isosurface associated with the -hole that has positive electrostatic potential. Only defined in case of a -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 B-X 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 hemi-sphere). 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 two-dimensional mapping of the electrostatic potential on the density contour is assembled. From this mapping, the properties above are calculated. Note that more than one -hole can be found. In the case of multiple -holes, if the positive-potential regions are in contact, the calculation of is not done. Also, it is important that the -hole does not exceed the maximum angle (MAXANG) and that NPTS is sufficiently high to have an accurate value.

    In addition to the list of -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, and 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 Exchange-Hole Dipole Moment Dispersion Model (XDM)

    The XDM keyword calculates the dispersion energy using the exchange-hole 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 {6|8|10}]
XDM [QE|POSTG] [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 3|6|sqrt6]
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 three-body 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 3-body damping. Use DAMP3BJN to choose the 3-body damping function: a product of three BJ-3 (3), three BJ-6 (6), or the square of three BJ-6 (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 exchange-hole dipole moment in the Becke-Roussel 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.

    • 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 all-electron density on a cube. If given, replaces the pseudo-density plus core in the calculation of the atomic volumes.

    • XB: the exchange-hole dipole moment in the Becke-Roussel model. Calculated from the above unless given.

    Only closed-shell (non-spinpolarized) systems can be calculated for now.

    The usual way of running XDM for the first time is:

    CRYSTAL rho.cube
LOAD rho.cube ZPSP C 4 O 6 H 1
LOAD elf.cube
XDM GRID elf 2

    This generates several cube files: root-tau.cube, root-pdens.cube, root-core.cube, root-lap.cube, root-grad.cube, and root-b.cube. Subsequent runs can circumvent the calculation of B, PDENS, and CORE by doing:

    CRYSTAL rho.cube
LOAD rho.cube ZPSP C 4 O 6 H 1
LOAD root-b.cube
LOAD root-pdens.cube
LOAD root-core.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 {6|8|10}

    Only calculate the contributions to the energy coming from the term (6), from the and terms (8) and from , , and (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 {EULER|HEUN|BS|RKCK|DP}] [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, Bogacki-Shampine’s embedded 2(3)th-order method with error estimation and first step as last. RKCK, Runge-Kutta-Cash-Karp embedded 4-5th order method with local extrapolation and error estimate. DP, Dormand-Prince 4-5th 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 = 1e-7, and maxerr.r = 1e-5.

    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 {GAULEG|QAGS|QNG|QAG}] [NR nr.i] [ABSERR err.r]
           [RELERR err.r] [ERRPROP prop.i] [PREC prec.r]

    The TYPE keyword selects the quadrature method:

    • GAULEG: Gauss-Legendre.

    • QAGS: quadpack’s dqags (general-purpose, extrapolation, globally adaptive, end-point singularities). All Q methods are sometimes unstable for heavy atoms and big beta-spheres, but this does not happen very often.

    • QNG: quadpack’s dqng (smooth integrand, non-adaptive, Gauss-Kronrod(Patterson))

    • QAG: quadpack’s dqag (general-purpose, integrand examiner, globally adaptive, Gauss-Kronrod)

    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: 7-15, 10-21, 15-31, 20-41, 25-51, 30-61. 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 = 1d-12. The default errprop is the field value and the Laplacian. delta.r = 1e-5 bohr.

    Integration Meshes (MESHTYPE)

    The MESHTYPE keyword chooses the type and quality of the molecular integration mesh:

    MESHTYPE {BECKE|FRANCHINI} [SMALL|NORMAL|GOOD|VERYGOOD|AMAZING]

    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: Becke-style integration grid. Only available for molecules.

    • FRANCHINI: a Becke-style 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:

    PRECISECUBE|STANDARDCUBE

    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.3). 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 {COV|VDW|} [at1.s|z1.i] rad1.r [[at2.s|z2.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)

    {RUN|SYSTEM} 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).

    Examples

    Updated: