Hello all,
I would like to compute electrostatic atomic charges fitted to the electrostatic potential generated with the density of an excited state. RESP charges seem to be already implemented (Restrained-Electrostatic-Potential), but I want to use an excited-state density matrix instead of the SCF default, to compute the electrostatic potential. How can I pass an excited state density matrix to the RESP solver to evaluate the electrostatic potential on the grid? Thanks a lot
You can change the QM method and basis set with these options here: https://github.com/psi4/psi4numpy/blob/master/One-Electron-Property/Restrained-Electrostatic-Potential/resp_driver.py#L171
However, not sure about PSI4’s capabilities for excited state densities for this.
If you provide your own ESP grid points in a file from another program, then it should work as well.
pinging the code author @alenaizan
I do not know about excited state densities, but as Holger mentioned you can provide your own grid and ESP values at the grid and the program will perform the RESP procedure for you. Check out the “grid” and “esp” options. The “grid” option takes a path to a file which has three columns for the x, y, and z coordinates. The “esp” option takes a path to a file with one column for the ESP calculated at the given grid points.
On another note, you may want to use the updated RESP package at: https://github.com/cdsgroup/resp. This code is better maintained and it has a few examples for using the package.
Thanks a lot for the input. I was able to have the RESP code running, but I still need to calculate the ESP with a different density matrix. Basically, I need to calculate one-electron integrals of the form
Integral_mn {Phi_m(r') Phi_n(r') / |r-r'| dr'}
where Phi_m, and Phi_n are the atomic basis functions, and |r-r’| denotes the modulus of the difference between the electron density position r’ and the observation point r.
In the language of pyscf it would be something like mole.intor('int1e_rinv'), but I am not sure if there is something similar for psi4. Thanks a lot!
It would be easier if you let us know more details what you are trying to do.
I dont know how pyscf exposes it’s integrals, this specific one is not ready-to-use exposed in PSI4.
At this moment what I am trying to do is calculating the potential “by hand”, so as I can then move to a different density. The electronic part of the electrostatic potential can be calculated via:
Vele = Sum_mn {P_mn Int_mn}
where P_mn is the element of the one-electron reduced density matrix between the atomic basis functions Phi_m and Phi_n.
I now only want to calculate the integral part by myself, at this moment only to reproduce an standard SCF calculation.
The integral that are expose to python in psi4 can be found here:
http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper
I think these http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.electric_field would be what you need. In our C++ code they are contracted with the density matrix (e.g. here https://github.com/psi4/psi4/blob/master/psi4/src/psi4/libmints/oeprop.cc#L1105)
Though for integral questions @andysim is the actual expert 
Sorry, didn’t see this until now. It looks like it’s the potential that’s needed, not the field (which is the gradient of the potential). This can be computed by asking for multipoles of order 0, as follows:
import psi4
mol = psi4.geometry("""
O
H 1 0.9
H 1 0.9 2 107
""")
psi4.set_options({'basis' : 'sto-3g'})
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('BASIS'))
mints = psi4.core.MintsHelper(wfn.basisset())
print(mints.ao_multipole_potential(origin=[0.0, 0.0, 0.1], max_k=0)[0].np)
print(mints.ao_multipole_potential(origin=[1.0, 0.0, 0.1], max_k=0)[0].np)
Here each call to ao_multipole_potential returns a list of matrix objects size 1 (hence the need to dereference with the [0]); these can be cast to numpy arrays via .np. The locations to compute the potential are specified in Bohr.