# Can the NCI (noncovalent interaction) index be generated starting with psi4?

Can NCI[1] index plots be generated starting from psi4? I'd be happy to just generate the NCI index values and mesh the isosurface myself.

With a long-term goal to visualise various medicinal-chemistry interactions with quantum chemical detail, it seems like the NCI index is the simplest starting point since it relies on just the electron density. Electron density can be generated using psi4 by outputting a cube file, or by querying the wave function at various points - which is how a cube file would be generated anyway.

How can one go from density to NCI index value? The second ingredient of NCI index is the reduced gradient, which might be obtained empirically with np.gradient on the density. The final ingredient is the second eigenvalue of the Hessian of the electron density at each evaluated point. However psi4 tells me Hessians aren't available for any functional or basis set I've tried (same thing happens with dimers):

import psi4
h2o = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
""")

nrg, wfn = psi4.energy('scf/cc-pVDZ', return_wfn=True)
wfn.compute_hessian()


output:

RuntimeError:
Fatal Error: Analytic Hessians are not available for this wavefunction.
Error occurred in file: /Users/github/builds/conda-builds/psi4-multiout_1557977521159/work/psi4/src/psi4/libmints/wavefunction.h on line: 321
The most recent 5 function calls were:


Alternatively the Hessian just describes the second derivative, so would calling np.gradient twice on the density be equivalent? Still, the Hessian at a single grid point is just three numbers - i.e. a 1d array, which has no eigenvalue.

Any tips here appreciated.

[1] Revealing Noncovalent Interactions, Johnson et al https://doi.org/10.1021/ja100936w

• Have you tried wfn.compute_gradient to see if it works? They may have only defined the basic interface for computing the Hessian, but may not have actually implemented it yet. Also, you want the gradient/hessian of the density, which is distinct from the gradient/hessian of the wavefunction. You may just want to generate the cube file of the density and then use numpy to manipulate that. Note that the Hessian at each point would be 9 values (a 3x3 array) and you would have diagonalize it at each point to get the 3 eigenvalues.
– Tyberius
Jun 9, 2021 at 15:55
• Thanks @Tyberius that clears up a lot. compute_gradient returns Fatal Error: Analytic gradients are not available for this wavefunction., but I think np.gradient will suffice. Jun 9, 2021 at 19:26
• I dont have time to write a full answer, but you might look into MultiWfn Jun 9, 2021 at 19:50
• Thanks @TristanMaxson , but I'm aiming to use psi4 to generate a repeatable and automatable code pipeline. Jun 9, 2021 at 19:57
• MultiWfn can actually be used in an automated way and can do NCI analysis (I have automated this before myself actually) Jun 10, 2021 at 1:12

Not a full answer of an alternative program, but just collecting the conclusions from the comments along with some additional suggestions.

It appears that compute_gradient and compute_hessian aren't implemented yet for wavefunction objects in Psi4. Even still, what you are interested in is the density and its gradient/hessian rather than those of the wavefunction itself.

Psi4 can generate a cube file of the density, which can be done by changing the last several lines of your input to:

psi4.set_options({"cubeprop_tasks": ["density"]})

E, wfn = psi4.energy('scf', return_wfn=True)
psi4.cubeprop(wfn)


Other cube file options, such as the spacing or overall range of points to include, are specified in the same way.

You should be able to readily extract the values to an (N,N,N) numpy array once you know the Cube format. From here, you can calculate the gradient and hessian using

grad=np.gradient(density,spacing)
h_eigs=[]
for i in range(N):
for j in range(N):
for k in range(N):
h=hess[:,:,i,j,k]
h_eigs.append(np.eig(h)[0])


which will give you a (3,N,N,N) array and a (3,3,N,N,N) array respectively. The code for the Hessian looks different because np.gradient assumes all axes of the array are variables of your functions, so just doing np.gradient(grad) would give you a (4,3,N,N,N) array.

The last step is getting the eigenvalues of the hessian matrix at each point in space. There is probably a better way to arrange this within numpy, but this gets across the concept that you need to iterate through every grid point and collect the eigenvalues.

• cheers! I've got a function to parse cube format files into a NumPy grid already so all that's left is to try it! I'll try and figure out how to broadcast the eig across the grid points too. Jun 9, 2021 at 21:39
• works perfectly. I'll post a gist below describing with the code Jun 10, 2021 at 8:19

If you call wfn.compute_hessian(), Psi4 tries to compute the nuclear Hessian which is indeed not implemented for most methods.

Calling np.gradient is not advisable, as it merely computes finite-difference approximations of derivatives instead of the true derivative. This is especially going to be a problem for cube files, since the resolution is poor from the beginning. A much better alternative is to compute the matrix of second derivatives analytically; this code exists in Psi4 as it is necessary for e.g. the calculation of DFT forces with GGA functionals.

The electron density is given in the basis set as $$n({\bf r}) = \sum_{\mu \nu} P_{\mu\nu} \chi_\mu({\bf r}) \chi_\nu({\bf r})$$ where $$\chi_\mu$$ is the $$\mu$$:th basis function and $${\mathbf P}$$ is the density matrix.

The $$j$$:th component of the density gradient is $$\partial_j n({\bf r}) = \sum_{\mu \nu} P_{\mu\nu} \left[ \chi_\mu({\bf r}) \{\partial_j \chi_\nu({\bf r})\} + \{\partial_j \chi_\mu({\bf r})\} \chi_\nu({\bf r})\right]$$ The $$i,j$$ element of the density Hessian is then $$\partial_i \partial_j n({\bf r}) = \sum_{\mu \nu} P_{\mu\nu} \left[ 2 \{\partial_i \chi_\mu({\bf r})\} \{\partial_j \chi_\nu({\bf r})\} + \{\partial_i \partial_j \chi_\mu({\bf r})\} \chi_\nu({\bf r}) + \{\partial_i \partial_j \chi_\nu({\bf r})\} \chi_\mu({\bf r}) \right]$$

You can get the alpha and beta density matrices with wfn.Da() and wfn.Db(); I'm not sure how to evaluate the basis functions, their gradients and laplacians.

In PySCF you can access the density matrix as dm = mf.make_rdm1() and get the values of the basis functions with the numint.eval_ao() function.

• Thanks - I think this might be a better solution longer term (although the finite difference approach does work fine for 0.1 angstrom spacing). It might just take me some time to figure out how pyscf works and how to connect it with those formulas :) Jun 12, 2021 at 3:25