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

Question

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

Answer 1

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)
hess=np.array([np.gradient(grad[i],spacing) for i in range(2)])
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.

Answer 2

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.