Interface to the PyFraME package for Polarizable Embedding

[JonasProg]

Hi,

I would like to write an interface in PySCF to the PyFraME package for embedding calculations, starting with polarizable embedding QM/MM. I want to enable SCF and geometry optimizations with PE (and possibly Hessian calculations).

Regarding the gradients necessary for the geometry optimizations (and second-order derivatives for the Hessian), I would like to know which intor to use to get geometric derivatives of the following one-electron electrostatic-potential integral and its potential derivatives:

t^{[\alpha]}_{\mu\nu} = \int \chi_\mu(r; R_n) T^{[\alpha]}(R_s, r) \chi_\nu(r; R_n) dr

where

T^{[\alpha]}(R_s, r) = \partial^{\alpha}_{r} \frac{1}{|r - R_s|}

Here $R_s$ is the position of a MM site, $R_n$ is the nuclear position, and $\alpha$ is a multi-index associated with the three cartesian coordinates.

The method is described in: https://doi.org/10.1021/acs.jctc.0c01323

[Original Issue https://github.com/pyscf/pyscf/issues/2214]

[Update 04.06.24] I've got the SCF and PE-SCF gradients for the geometry optimizations running now. But I would need some pointers to get the integrals for Hessian calculations.

[MatthewRHermes]

I'm not sure I'm interpreting the calculus notation correctly ($\partial^\alpha_r = \frac{\partial r}{\partial \alpha}$?) but it looks to me like

mol.set_rinv_orig (Rs)
TRs = mol.intor ('int1e_iprinv', comp=3)

from gto.moleintor and libcint

[MatthewRHermes]

To get higher-order derivatives, I think you can just put more ip in front of rinv:

TRs2 = mol.intor ('int1e_ipiprinv', comp=9) # second derivatives
TRs3 = mol.intor ('int1e_ipipiprinv', comp=27) # third derivatives
TRs4 = mol.intor ('int1e_ipipipiprinv', comp=81) # fourth derivatives
# according to my installation's cint_funcs.h file, this is the highest-order derivative available