ESP charges#
Let us consider a molecular system with \(M\) nuclei and an electron charge distribution described by the one-particle density matrix \(D\).
Our goal is to represent the quantum mechanical electrostatic potential
with a classical potential associated with a set of atomic partial charges known as the electrostatic potential (ESP) charges
There is of course no unique way of doing so. The Merz–Kollman (MK) scheme [BMK90, SK84] minimizes the squared norm difference between these two potentials evaluated on a set of grid points in the solvent-accessible region of the molecule with respect to variations in the partial charges and a constraint of a conservation of the total molecular charge, \(Q\). The grid points are distributed on successive layers of scaled van der Waals surfaces at positions \(\mathbf{r}_a\).
To achieve this minimization, a Lagrangian is defined
where
Finding the stationary point of this Lagrangian results in the following equation from which the ESP charges are obtained
where
and
The MK scheme is implemented in VeloxChem and we will use this implementation to determine the ESP charges of methanol at the Hartree–Fock/6-31G(d) level of theory.
Show code cell source
import py3Dmol as p3d
methanol_xyz = """6
H 1.2001 0.0363 0.8431
C 0.7031 0.0083 -0.1305
H 0.9877 0.8943 -0.7114
H 1.0155 -0.8918 -0.6742
O -0.6582 -0.0067 0.1730
H -1.1326 -0.0311 -0.6482
"""
viewer = p3d.view(width=300, height=200)
viewer.addModel(methanol_xyz, "xyz")
viewer.setStyle({"stick": {}})
viewer.show()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
import veloxchem as vlx
molecule = vlx.Molecule.read_xyz_string(methanol_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31g*")
A calculation of ESP charges is performed with the compute
method of the RespChargesDriver
class.
esp_drv = vlx.RespChargesDriver()
esp_charges = esp_drv.compute(molecule, basis, "esp")
print("Atom ESP charge")
print(20 * "-")
for label, esp_charge in zip(molecule.get_labels(), esp_charges):
print(f"{label :s} {esp_charge : 18.6f}")
print(20 * "-")
print(f"Total: {esp_charges.sum() : 13.6f}")
Atom ESP charge
--------------------
H 0.056273
C 0.196542
H -0.007483
H -0.008463
O -0.657469
H 0.420600
--------------------
Total: -0.000000