Overview
--------

NCIPLOT is a program that enables the computation and graphical
visualization of inter- and intra-molecular non-covalent interactions
(hydrogen bonds, pi-pi interactions, ...) in systems ranging from
small dimers to large biomolecules.

NCI (Non-Covalent Interactions) is a visualization index based on
the electron density ($\rho$) and the reduced density gradient (RDG, $s$),
defined as:

$s = \frac{1}{2(3\pi^2)^{1/3}}\frac{|\nabla\rho |}{\rho^{4/3}}$

The non-covalent contacts can be identified with the regions of small
reduced density gradient (RDG) at low densities. These regions are
mapped in real-space by plotting an isosurface of $s$ for a low value
of RDG. In addition, the sign of the second eigenvalue of the density
Hessian times the density, is color-mapped onto the isosurfaces,
allowing to characterize both the strength and (un)favourable nature
of these interactions. The references below give a thorough
explanation of the details.

NCIplot reads a gaussian-style wavefunction file (wfn or wfx) or a
molecular geometry from an xyz file and computes the density and the
reduced density gradient on a grid. The calculated values are written
to a gaussian cube file, along with VMD scripts that allow the direct
visualization of the results.

Compilation
-----------

To install nciplot, unpack the contents of the distribution and cd
into the src/ subdirectory:

    tar xjvf nciplot.tar.bz2
    cd nciplot/src/


Then modify the Makefile.inc to suit your compiler and flags. Two
examples (gfortran and ifort) are given in the default
Makefile.inc. For instance, the gfortran flags are:

    FC = gfortran
    FCFLAGS = -O -fopenmp
    LDFLAGS = -fopenmp

Note the -fopenmp flag that enables shared-memory parallelization.
Once you are satisfied with your choice of compiler options, do:

    make mrproper
    make

to build the nciplot executable. To clean the object and module files,

    make clean

To clean objects, modules, and binaries,

    make mrproper

NCIplot requires the NCIPLOT_HOME environment variable to find the
atomic density files contained in the dat/ subdirectory. Set this
variable to the absolute path to NCIplot:

    export NCIPLOT_HOME=/home/xxxx/programs/nciplot

so that the dat directory is ${NCIPLOT_HOME}/dat/ and a link to the
binary is located in ${NCIPLOT_HOME}/bin/. Several tests and examples
can be downloaded from:

http://schooner.chem.dal.ca/downloads/nciplot/

The old PDF manual (outdated) can be found in doc/.

The code has been parallelized for shared-memory architectures using
the OpenMP library. To use this feature, set the OMP_NUM_THREADS
environment variable to the number of cores:

    export OMP_NUM_THREADS=4

Use
---

NCIplot accepts a single input file (typically with nci
extension). Run it as:

    nciplot example.nci

The standard output can be redirected by writing the output file as
the second argument:

    nciplot example.nci example.out

In addition, depending on the input, some files may be created:

- example.dat: a 2-column file containing the reduced density gradient
  (col. 2) vs. density (col. 1) calculated at the points of the 
  grid. This grid is the same as the one used in the files below.
- example-dens.cube : a cube with the electron density.
- example-grad.cube : a cube with the reduced density gradient
  (RDG).
- example-elf.cube : a cube containing the electron localisation
  function (ELF). The generation of this file is activated using the
  ELF keyword.
- example.vmd : the VMD script for convenient visualization of the
  results. 

The simplest possible NCIplot input is:

    1
    molecule.wfn

where the 1 stands for the number of molecules to be read in the input
and molecule.wfn a the wavefunction file generated by gaussian, which
will be used to compute the density and RDG. If no wavefunction file
is available (because, for instance, the molecule is too big) then:

    1
    molecule.xyz

reads the molecular geometry from the xyz file and then uses the
promolecular density for the visualization.

Input syntax
------------

After the number of molecules and the wavefunction (or xyz) file, a
number of keywords can be given to NCIplot to control its behavior
(grid size, which regions to plot, etc.)

Regarding the limits of the grid, by default, the origin is chosen as
the minimum (x,y,z) coordinates of all the molecules in input, and the
upper corner is set to the maximum (x,y,z) triplet. Then a small
quantity (2 bohr, can be changed with the RTHRES keyword) is added in
each direction. The resulting cube encompasses all the input
molecules.

In the following, all lengths are in angstrom and densities in e/bohr^3. 

; RTHRES h.r

: Change the thickness of the border around the default cube. 
  Default: 2 bohr.

; RADIUS x.r y.r z.r h.r 

: Sets the grid limits to (x-h,y-h,z-h) and (x+h,y+h,z+h). 

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

: Sets the grid limits to (x0,y0,z0) and (x1,y1,z1). 

; LIGAND nfile.i h.r

: Takes the molecule number nfile.i and calculates (x0,y0,z0) as the
  minimum of all the (x,y,z) atomic coordinates; and (x1,y1,z1) as the
  maximum. Then adds a region of widht h.r to those values:
  (x0-h,y0-h,z0-h) and (x0+h,y0+h,z0+h). The final cube will contain
  the molecule number nfile.i and a small region of width h.r around
  it. This is useful in the study of ligands embedded in larger
  molecules. The use of this keyword also deactivates writing the
  'intramolecular' points to the dat file (see INTERMOLECULAR).

; ATCUBE

     ATCUBE
       file1.i at1.i at2.i at3.i ...
       file2.i at4.i ...
     END | ENDATCUBE

: Calculate the cube that encompasses atoms from different molecules:
  atoms at1.i, at2.i,... from file file1.i; at4.i from file2.i,
  etc. Then add a region of width 2 bohr around the cube.

; EXC xfun.s cfun.s

: Write the exchange-correlation energy density to a cube file
  (-xc.cube). The first string selects the exchange part and the
  second one is the correlation part. The exchange contribution may be
  one of:

  + s, slater, lda : homogeneous electron gas exchange.
  + b88 : Becke88 exchange.
  + pbe : PBE exchange.
  + `-` : none.
  
: And the correlation part:

  + pw | lda : homogeneous electron gas correlation, in the Perdew-Wang
               parametrization. 
  + pbe : Perdew-Burke-Ernzerhof (96) correlation.
  + b88 : Becke88.
  + `-` : none.

The number of points along each side of the grid is set using:

; INCREMENTS x.r y.r z.r

: x.r, y.r and z.r are the step lengths.  The number of points
  is calculated as (xend - xini) / xinc, where xini and xend are the
  limits above.

The INTERMOLECULAR keyword can be used to avoid writing the
(density,rdg) points corresponding to intramolecular interactions to
the dat file. Specifically, a point is not represented if the density
comes mostly (by default, more than 95%) from only one of the
molecules in the input. The fraction can be controlled using the
RHOCUT2 keyword:

; RHOCUT2 f.r g.r

with f.r = 0.95 and g.r = 0.75 by default (see right below for an
explanation of g.r).

By default, the fragments to which INTERMOLECULAR applies are the
molecules defined in input. This behavior can be modified using the
**FRAGMENT** keyword:

    FRAGMENT
      file1.i at1.i at2.i at3.i ...
      file2.i at4.i ...
    END | ENDFRAGMENT
 
Several FRAGMENT environments can appear in the same input, each
defining a different fragment. The same atom should appear only in one
of the defined fragments. If the combined list of all the atoms inside
a FRAGMENT environment is not exhaustive, then the remaining atoms are
not considered. Inside a FRAGMENT environment, atoms can be defined by
choosing the molecule (ifile1.i) and then the list of atoms within
that molecule (at1.i at2.i etc.). If the FRAGMENT keyword is used,
then the INTERMOLECULAR keyword is automatically activated. Note that
if the EXC keyword is used, then the criterion applies also to
that cube. That is, if the considered point is intramolecular (more
than 95% of the density coming from only one fragment), then the value
of the exc is set to a very high value to remove the point from the
isosurface representation. 

In the case when the are atoms not assigned to any fragment, then the
second RHOCUT2 parameter (g.r) applies. A point is written to the dat
and cube files if: 

$\rho_i < f.r * (\sum \rho_i)$

where i runs over fragments and:

$(\sum \rho_i) > g.r * \rho_{total}$

This way, only the points close to the selected fragments are
considered for plotting.

A discussion on the density and rdg cutoffs follow. A point is
represented in the dat file if the density is below rhocut.r
(def. 0.2), the rdg is below dimcut.r (default: 2.0 if any molecule
comes in wfn format, 1.0 otherwise) and the point is not considered
intramolecular. These two parameters can be controlled using the
keyword:

; CUTOFFS rhocut.r dimcut.r

A point is represented in the rdg cube file (-grad.cube) if the
density is below rhoplot.r and the point is not considered
intramolecular. The rhoplot.r value is set by the keyword CUTPLOT: 

; CUTPLOT rhoplot.r 

When the density is greater than rhoplot.r, the rdg in the -grad.cube
file is set to 100d0, effectively eliminating the point from the
isosurface plot. The default value is 0.05 (if at leamax one file is
wfn) or 0.07 (otherwise). 

The rdg iso-value represented in the VMD script is set by the ISORDG
keyword: 

; ISORDG isordg.r

The default is 0.5 (any molecule is a wfn) or 0.3 (all molecules are
xyz). 

A point is represented in the exc cube file (-xc.cube) if it is not
considered intramolecular. All points are represented in the elf cube
file (-elf.cube).

The number of files on output can be controlled using the keywords:

; OUTPUT n.i

: n.i can be 1, 2 or 3. The files written in each case are:

  + 1: only dat.
  + 2: only cubes and vmd script.
  + 3: all.

: The default is 3.

; ELF

: Generates a -elf.cube file containing the ELF scalar field. Only wfn
  files contribute to ELF.

; EDR d.r

: Generate the electron delocalization range function EDR(r;d) in -edr-d.cube 
: Janesko, Scalmani, Frisch, J Chem Phys 2014, 141, 144104
: d.r is a distance d in bohr. 

; EDRDMAX ainc.r amax.r anum.i 

: Generate D(r), the distance d maximizing EDR(r;d) at point r, in -D.cube 
: Janesko, Wiberg, Scalmani, Frisch, J Chem Theory Comput 2016, 141, 144104 
: Generated by interpolating EDR(r;d.i) at several d.i values 
: Generates an even-tempered set of anum.i different exponents a=1/d^2, starting at 
: max value amax.r (bohr^(-2)), with increment between values ainc.r . 

Finally, there are some miscellaneous keywords:

; ONAME root.s

: By default, the root of the input filename is used to
  generate the output file names. With ONAME, it is possible to change
  this root.

; DGRID

: There are two sets of promolecular electron densities that can be
  used in NCIplot: exponential fits and radial grids. The exponential
  fits are available up to Z=18 (Ar, first 3 rows of the periodic
  table) while the radial grids are available up to Pu (Z=94). The
  default behavior is using the exponential fits unless the molecule
  contains some atom with Z>18 or there are charged atoms
  (cations). Using DGRID, only radial grids are used to calculate
  promolecular densities. 

For more information, please refer to the manual in the doc/
subdirectory. 

References and citation
-----------------------

The basic references for nciplot are:

* J. Contreras-Garcia, E. Johnson, S. Keinan, R. Chaudret, J-P
  Piquemal, D. Beratan, W. Yang, J. Chem. Theor. Comp. **7**, 625
  (2011). (http://dx.doi.org/10.1021/ct100641a)
* E R. Johnson, S. Keinan, P. Mori-Sanchez, J. Contreras-Garcia, A
  J. Cohen, and W. Yang, J. Am. Chem. Soc., **132**, 6498 (2010). 
  (http://dx.doi.org/10.1021/ja100936w)

Some other works using NCI:

* J. Contreras-Garcia, E. R. Johnson, W. Yang, J. Phys. Chem. A
  **115**, 12983 (2011). (http://dx.doi.org/10.1021/jp204278k)
* J. Contreras-Garcia, M. Calatayud, J-P Piquemal, J. M. Recio,
  Comp. Theo. Chem. **998**, 193 (2012).
  (http://dx.doi.org/10.1016/j.comptc.2012.07.043)
* N. Gillet, R. Chaudret, J. Contreras-Garcia, W. Yang, B. Silvi,
  J.-P.Piquemal, J. Chem. Theor. Comp. **8**, 3993 (2012). 
  (http://dx.doi.org/10.1021/ct300234g)

NCI for the solid state:

* A. Otero-de-la-Roza, J. Contreras-Garcia, E. R. Johnson,
  Phys. Chem. Chem. Phys. **14**, 12165 (2012).
  (http://dx.doi.org/10.1039/C2CP41395G)
* E. R. Johnson, A. Otero-de-la-Roza, J. Chem. Theory Comput. **8**,
  5124 (2012). (http://dx.doi.org/10.1021/ct3006375)

NCI for experimental densities:

* G. Saleh, C. Gatti, L. Lo Presti, J. Contreras-Garcia,
  Chem. Eur. J. **18**, 15523 (2012). 
  (http://dx.doi.org/10.1002/chem.201201290) 

Copyright notice
----------------

Copyright (c) 2013-2015 Alberto Otero de la Roza,
Julia Conteras-Garcia, Erin R. Johnson, and Weitao Yang

nciplot is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.