# Making NCI Plots in molecules, surfaces, and solids

These examples show how to make non-covalent interaction (NCI) plots in molecules and solids. For more details on what NCI plots are, see Johnson et al.’s article and this other article for the implementation in solids, and the manual pages indicated at the end.

## Molecules: phenol dimer

NCI plots can be generated in two ways. If you have the structure and the self-consistent electron density, you can calculate the NCI domains using them. Otherwise, the promolecular density can be used for the generation of the NCI domains, in which case only the structure is required.

A simple NCI plot input file for phenol dimer using the promolecular density is:

molecule phenol_phenol.xyz 0

nciplot
endnciplot


First, the molecular structure of the phenol dimer (contained in an xyz file) is loaded into critic2 with the MOLECULE keyword. Then, the NCI plot is generated by calling the nciplot/endnciplot environment. A number of optional keywords can be used in this environment to tune the behavior of NCIPLOT, a few of which are described in this example. Please consult the NCIPLOT keyword entry in the manual for details.

The 0 after the MOLECULE command instructs critic2 to restrict the molecular cell to the smallest box that encompasses the molecular structure (i.e. a box with zero border). The same effect can be achieved using the MOLCELL keyword:

molcell 0


after loading the structure with MOLECULE. The molecular cell determines the size of the grid in which the NCI plot quantities (density, reduced density gradients, etc.) will be calculated. A small molecular cell has the double advantage of making the calculation cheaper (because you need smaller grids) and also of restricting the NCI plot to the region of interest. The only downside is that one of your NCI domains may be cut off, in which case you will have to increase this number somewhat.

Running the example above generates a number of files. You can view the results using the vmd program. To open the file, use:

vmd -e phenol_phenol.vmd


Rotate the molecule into the position you want and do File -> Render -> Tachyon (internal, in memory) and press “Start rendering”. This will generate the NCI plot file.

In general, it is better if you have run the calculation and have the self-consistent density for the system from which to make the NCI plot. This density can be read in any format understood by critic2. In this example, we can use a wfx wavefunction file from Gaussian:

molecule phenol_phenol.wfx 0
nciplot
endnciplot


Note the only difference is that the structure and density are both read from the wavefunction file, and therefore used in the NCI plot.

## Surfaces: benzene on kaolinite

Now we consider making the NCI plots of a molecule (benzene) on an inorganic surface, kaolinite (001). The equilibrium structure was calculated using Quantum ESPRESSO in two cases: when the benzene molecule is adsorbed to the hydrophobic side (hb) and when it is adsorbed to the hydrophilic side (hl). The structure and self-consistent electron densities are contained in the hbpar1_rho.cube and hlpar1_rho.cube cube files, respectively, generated using QE’s pp.x program.

We consider absorption to the hydrophobic surface first (hbpar1_rho.cube). We could use an NCIplot input file identical to the one in the previous section for molecular systems. However, in surfaces and solids, there are typically too many NCI domains and it is often the case that we want to focus on a particular interaction. In this case, we want only the NCI domains that arise from the interaction between the molecule and the surface. To focus only on these interactions, we create .xyz files containing the relevant fragments (one for benzene and another for the surface) and then pass those to the NCIPLOT environment. That way, only the domains that involve interactions between the two fragments will be plotted.

Therefore, we first need to write the structure contained in the cube files and select the fragments for the NCIplot calculation. To do this, we first write an xyz file containing the surface and the molecule:

crystal hbpar1_rho.cube
write cell.xyz border 1 1 2


This creates an xyz file (cell.xyz) containing the slab plus the molecule. The 1 1 2 instructs critic2 to write a 1x1x2 supercell of the original system, so that a whole slab is contained in the file. The BORDER keyword makes all atoms at the border of the supercell be included in the file. You can now open this file with your favorite program (I use avogadro) and cut away from it the benzene and the surface underneath, and save the two systems to separate xyz files (hbpar1_benzene.xyz and hbpar1_surface.xyz).

Lastly, we generate the NCI plot with the input file:

crystal hbpar1_rho.cube
nciplot
fragment hbpar1_benzene.xyz
fragment hbpar1_surface.xyz
nstep 100 100 100
endnciplot


We first load the surface structure from the cube file, then the self-consistent electron density, which becomes the reference field and therefore will be used by NCIPLOT to generate the NCI domains. The xyz files containing the fragments written in the previous step are passed to the NCIPLOT environment, which ensures that only the NCI domains contained between those fragments are plotted. Lastly, the default grid used for the NCI plots (0.1 bohr between points in each direction) is too large for this example. Therefore, we use the NSTEP option to limit the grids used for the plot to 100 points in each direction. In a production calculation, you should use a finer grid (perhaps the default grid) to improve the quality of the plot and avoid artifacts.

For the case of benzene on the hydrophilic surface, the procedure to make the NCI plots is entirely analogous. First, we create an xyz file containing the surface:

crystal hlpar1_rho.cube
write cell.xyz border 1 1 2


Then, we cut the benzene and surface fragments into their respective xyz files (hlpar1_benzene.xyz and hlpar1_surface.xyz) and finally we generate the NCI plot with:

crystal hlpar1_rho.cube
nciplot
fragment hlpar1_benzene.xyz
fragment hlpar1_surface.xyz
nstep 100 100 100
endnciplot


As in the molecular case, the promolecular density can be used instead of the calculated electron density if the latter is not available. To do this, simply remove the LOAD command. This requires only the surface structure.

## Solids: thymine molecular crystal

Same as in the surface case, NCI plots of a periodic crystal are very busy because of the large number of NCI domains and atoms. We need to use the FRAGMENT keyword to focus the NCI plot on the interactions we are interested in. In this example we examine the thymine crystal, in which there are two dominant interactions: the intra-layer hydrogen bonds and the inter-layer stacking interactions.

Before calculating the NCI plots, we need to generate the fragment structures. The equilibrium structure and electron density were calculated with Quantum ESPRESSO, and both are contained in the thymine_rho.cube file. To create the fragments, we use the following input:

crystal thymine_rho.cube
write cell.xyz molmotif


This creates an xyz file (cell.xyz) containing the atoms in the unit cell. The MOLMOTIF option makes critic2 include all atoms from neighboring cells necessary to complete all molecular structures. If a larger unit cell is needed to define the fragments, we could do:

crystal thymine_rho.cube
write cell.xyz molmotif 2 2 2


to perform the same operation on a 2x2x2 supercell. You can also write all the monomers in the unit cell with:

write thymine_monomers.xyz molmotif nmer 1


and all the dimers with:

write thymine_pairs.xyz molmotif nmer 2


Using avogadro, we now cut the fragments of interest. To examine the intra-layer hydrogen bond, we select three adjacent molecules interacting via hydrogen bonds within the same layer (thymine_intra*.xyz). For the stacked layer interaction, we pick a three-molecule motif from one layer (thymine_inter1.xyz) and a single molecule, exactly above this triangular motif, from an adjacent layer (thymine_inter1.xyz).

Now we generate the NCI plot for the inter-layer stacking interaction with:

crystal thymine_rho.cube

nciplot
fragment thymine_inter1.xyz
fragment thymine_inter2.xyz
nstep 100 100 100
endnciplot


The two fragments we generated before (the triangular motif and the single molecule on top of it) are passed to the NCIPLOT environment using the FRAGMENT option. This restricts the NCI domain to only those that appear between the selected fragments. The NSTEP option reduces the number of points in the calculated density and reduced density gradient grids. For production calculations, you should use the defaults (remove the NSTEP) or increase these numbers; otherwise oscillations appear, as shown in the NCI plots below.

For the intra-layer interactions, we pass to the FRAGMENT option the three fragments corresponding to the three hydrogen-bonded molecules within the same layer:

crystal thymine_rho.cube

nciplot
fragment thymine_intra1.xyz
fragment thymine_intra2.xyz
fragment thymine_intra3.xyz
nstep 175 175 75
endnciplot


As before, we use NSTEP to restrict the number of points in the calculated grids. You should use the defaults, or at least a finer grid, for production plots.

As before, the promolecular density can be used instead of the calculated electron density if the latter is not available. To do this, simply remove the LOAD commands. In this case, only the structure of the molecular solid is needed.

## Example files package

Files: example_12_01.tar.xz.

Updated: