• Non-covalent interaction (NCI) plots

    Introduction

    The non-covalent interaction plots are graphical representations of the regions where the densities of two interacting molecules overlap. The accumulation of density in the intermolecular space, which is related to the amount of intermolecular repulsion, can be used to qualitatively infer the strength of the interaction. The implementation of NCI plots in critic2 works with molecules and with periodic solids. NCI plots can be created using the NCIPLOT keyword:

    NCIPLOT
  ONAME root.s
  CUTOFFS rhocut.r dimcut.r
  RHOPARAM rhoparam.r
  RHOPARAM2 rhoparam2.r
  CUTPLOT rhoplot.r dimplot.r
  SRHORANGE srhomin.r srhomax.r
  VOID void.r
  RTHRES rthres.r
  INCREMENTS x.r y.r z.r
  NSTEP nx.i ny.i nz.i
  ONLYNEG
  NOCHK
  CUBE x0.r y0.r z0.r x1.r y1.r z1.r
  CUBE file1.xyz file2.xyz ...
  MOLMOTIF
  FRAGMENT file.xyz
  FRAGMENT
   x.r y.r z.r # (in angstrom!)
   ...
  ENDFRAGMENT/END
ENDNCIPLOT/END

    The NCIPLOT keyword calculates the density () and the reduced density gradient: \begin{equation} s({\bf r}) = \frac{\lvert{\bf \nabla}\rho({\bf r})\rvert}{2(3\pi^2)^{1/3}\rho({\bf r})^{4/3}} \end{equation} on a three-dimensional grid spanning the system or a certain region of the system determined by the user. The density is assumed to be the reference field (see the REFERENCE keyword). If no scalar fields have been loaded, then the promolecular density is used instead. The reduced density gradient is calculated from the reference fields using the formulas above.

    The output of the NCIPLOT keyword is meant to be visualized with the vmd program. As a result of its operation, NCIPLOT generates the following files:

    • <root>.dat: a text file containing two columns of numbers. The first column is the density multiplied by the sign of the second eigenvalue of the density Hessian. The second is the reduced density gradient. Each row corresponds to one point in the NCIPLOT grid. Plotting this data file directly (for instance, with gnuplot) gives the characteristic $s$ vs. $\rho$ plots seen in the literature.

    • <root>-dens.cube: a cube containing the electron density times the sign of the second Hessian eigenvalue times 100.

    • <root>-grad.cube: a cube containing the reduced density gradient (RDG).

    • <root>.vmd: the VMD script for convenient visualization of the results. To use this file, open vmd from the same directory where the .vmd file resides. Then go to File->Load Visualization State and select the .vmd file.

    The root of all the files generated by the NCIPLOT run can be changed using the ONAME keyword in the NCIPLOT environment (default: the root of the critic2 input file).

    Worked out examples showing how to make NCIPLOTs in molecules, surfaces, and crystals can be found in the examples pages. This may be an easier introduction to the use of the NCIPLOT keyword. The full set of options is described below.

    Changing the Number of NCI Grid Points

    By default, the region represented is the whole unit cell (in a solid) or the whole molecule, with step lengths of 0.1 bohr in each direction. The step lengths or the number of points in each axis can be controlled with the INCREMENTS and NSTEP keywords. The INCREMENTS keyword:

    INCREMENTS x.r y.r z.r

    accepts three arguments x.r, y.r, and z.r, which are interpreted as the step lengths for the NCIPLOT grid in each direction. The units are bohr (crystals) or angstrom (in molecules). The NSTEP keyword:

    NSTEP nx.i ny.i nz.i

    accepts three integers corresponding to the number of steps in each NCIPLOT grid direction.

    Selecting the Region for the NCI Plot

    In medium or large molecules and crystals, representing the NCI domains for the whole structure results in a very busy plot. In those cases, it is recommended that the user focuses on the region of interest and discards all the other domains. The simplest way to do this is to restrict the NCIPLOT to a box encompassing the required region:

    CUBE x0.r y0.r z0.r x1.r y1.r z1.r

    where (x0.r y0.r z0.r) and (x1.r y1.r z1.r) are the coordinates corresponding to the two opposite corners that determine the orthogonal box in which the NCI domains will be calculated. The coordinates are in crystallographic (crystals) or in molecular Cartesian coordinates (molecules, default unit: angstrom). In crystals, a region larger than one periodic cell can also be selected with the CUBE keyword.

    Another option for restricting the NCI plot to a certain region encompassing a number of intermolecular contacts is ues the box that contains a number of fragments of the crystal or molecular structure using:

    CUBE file1.xyz file2.xyz ...

    The box for the NCI grid will contain the fragments given by files file1.xyz, file2.xyz,… in standard xyz format. A small border, given by the RTHRES keyword, is added around the region exactly containing all the atoms in all the xyz files. Note that, unlike in the FRAGMENT keyword, the fragments passed to the CUBE keyword are not used as fragments for the NCI calculation. They are simply used to establish the box for the NCI grid.

    Cutoffs and Options for the Plot

    Some cutoffs are relevant to the visualization of the NCI. The rhocut.r and dimcut.r options of the CUTOFFS keyword:

    CUTOFFS rhocut.r dimcut.r

    determine how the density and reduced density gradient are written to the .dat file (respectively). If, at a given point, the density is above rhocut.r or the reduced density gradient is above dimcut.r, then the point is not written to the .dat file. The defaults are: rhocut.r = 0.2 and dimcut.r = 2.0 for densities from a loaded scalar field and 1.0 for promolecular densities.

    The cutoffs rhoplot.r and dimplot.r in the CUTPLOT keyword:

    CUTPLOT rhoplot.r dimplot.r

    can be used to change how (and whether) the NCI domains are represented. When the density is greater than rhoplot.r, the reduced density gradient at that point in the -grad.cube file is set to a value of 100, effectively eliminating that point from the isosurface plot. In addition, the default color scale represented in RDG isosurfaces ranges from -rhoplot.r to rhoplot.r. The default value is 0.05 (densities from a loaded scalar field) or 0.07 (promolecular). The dimplot.r value is the default isosurface value represented by VMD. Default: 0.5 (loaded scalar field) or 0.3 (promolecular).

    The SRHORANGE keyword:

    SRHORANGE srhomin.r srhomax.r

    is similar to CUTPLOT. When the density times the sign of the second Hessian eigenvalue is greater than srhomax.r or smaller than srhomin.r, the reduced density gradient in the -grad.cube file is set to 100, effectively eliminating the point from the isosurface plot. This is useful when plotting a subset of the peaks found in the .dat file as isosurfaces. By default, no srhomin.r or srhomax.r pruning is done.

    The ONLYNEG keyword:

    ONLYNEG

    tells critic2 to represent only the points where the second eigenvalue of the Hessian is negative.

    In molecular crystals the MOLMOTIF keyword:

    MOLMOTIF

    is used to complete the molecules that lie across unit cell faces by using atoms from neighboring cells. This results in a much more understandable plot, since all molecules are wholly represented.

    Finally, the VOID keyword:

    VOID void.r

    is used to represent only the points where the promolecular density is lower than void.r.

    Checkpoint Files

    By default, an NCI plot will create a checkpoint file with the same name as the root for the run and extension .chk_nci. This checkpoint is useful if several similar test runs need to be made, for instance, using different cutoffs. To prevent critic2 from writing the checkpoint file, use:

    NOCHK

    The NCIPLOT checkpoint file contains the promolecular densities for the whole system and the fragments, the density, and the reduced density gradient. If the checkpoint file exists and is compatible with the current structure and fragments, it is read. Otherwise, it is discarded and a new checkpoint is written. By default, checkpoint files are used.

    NCI Domains Only Between Selected Fragments

    In very large systems, the default NCI plots generated by critic2 contain all NCI domains. This results in a very busy plot and one normally wants to focus on a specific intermolecular interaction. In the current version of NCIplot, it is possible to define molecular fragments in order to focus on a specific interaction between two fragments of the system. This can be done by using the FRAGMENT environment:

    FRAGMENT
 x.r y.r z.r # (in angstrom!)
 ...
ENDFRAGMENT/END
FRAGMENT file.xyz

    Each FRAGMENT block or keyword defines one fragment. Only the intermolecular interactions between fragments are represented (hence, at least two FRAGMENT blocks are necessary). In the FRAGMENT environment, the atomic positions of all the atoms in the fragment are given. The x.r, y.r, z.r are Cartesian coordinates, and the units are angstrom, regardless of the use of the UNITS keyword. This is useful in combination with the WRITE keyword. In the typical usage, WRITE generates an xyz file containing a portion of the crystal or molecule. Then, the user removes all the undesired atoms using a graphical program (like avogadro). Finally, the xyz of the desired fragment are passed back to critic2 using the FRAGMENT keyword. To simplify this process, the fragment can be read directly by FRAGMENT from an xyz file (file.xyz). Note that the fragments need not include all atoms in the crystal or molecule, or may contain several lattice-translated copies of the same atom.

    There are three options that control the behavior of the fragments: RTHRES, RHOPARAM and RHOPARAM2. The RTHRES keyword is:

    RTHRES rthres.r

    When fragments are used, the density and reduced density gradient grids are calculated inside a box encompassing the fragments, with an additional border that is rthres.r thick. This value has units of bohr in crystals, and angstrom in molecules. The default is 2.0 bohr. Note that RTHRES can also be used to control the border around NCI regions defined with the CUBE keyword.

    The RHOPARAM keyword:

    RHOPARAM rhoparam.r

    restricts the representation of NCI domains to those regions where none of the fragments contributes more than rhoparam.r times to the total promolecular density. The default is 0.95 (i.e. a 95%). This keyword is useful to remove intramolecular domains within a fragment.

    The RHOPARAM2 keyword:

    RHOPARAM2 rhoparam2.r

    is used to restrict the representation of NCI domains to the regions where the sum of the promolecular density from all defined fragments is higher than rhoparam2.r times the total promolecular density of the whole system. The default is 0.75 (i.e. 75%). This keyword is useful to remove NCI domains that are unrelated to the selected fragments.

    Tips for Using NCIPLOT

    Some tips for using the NCIPLOT keyword efficiently:

    • If the reference field is a grid, the system is periodic (i.e. the structure was read with CRYSTAL) and no FRAGMENT keyword is used, the dimensions of the grids used in NCIPLOT are the same as the reference field.

    • If the reference field that provides the density for NCIPLOT is a grid field, then it is usually much faster if core-augmentation (ZPSP field option) is not used. The reason is that if the core augmentation is not present, the reduced density gradient and the Hessian components are calculated by Fourier transform, leading to smoother derivatives. However, not using the core augmentation can sometimes result in noisy Hessian components, which may result in alternate blue/red domains (because of a spurious change of sign caused by numerical inaccuracies).

    • Older versions of VMD have problems dealing with non-orthogonal cells. There’s little critic2 can do about this. To work around this problem, use the FRAGMENT or CUBE keywords so an orthogonal NCI region is used.

    Examples

    Updated: