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
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. $\endgroup$compute_gradient
returnsFatal Error: Analytic gradients are not available for this wavefunction.
, but I thinknp.gradient
will suffice. $\endgroup$psi4
to generate a repeatable and automatable code pipeline. $\endgroup$