Can NCI 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()
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.
 Revealing Noncovalent Interactions, Johnson et al https://doi.org/10.1021/ja100936w