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

\[\begin{equation*} V(\mathbf{r}) = \sum_{A=1}^M \frac{Z_A e}{4\pi\varepsilon_0 |\mathbf{r}-\mathrm{\textbf{R}}_A|} - e \sum_{\alpha,\beta} D_{\alpha\beta} \int \frac{ \phi_\alpha^*(\mathbf{r}')\phi_\beta(\mathbf{r}') }{ 4\pi\varepsilon_0 |\mathbf{r}-\mathbf{r}'| } d^3\mathbf{r}' \end{equation*}\]

with a classical potential associated with a set of atomic partial charges known as the electrostatic potential (ESP) charges

\[\begin{equation*} \tilde{V}(\mathbf{r}) = \sum_{A=1}^M \frac{ q_A }{ 4\pi\varepsilon_0 |\mathbf{r}-\textbf{R}_A| } \end{equation*}\]

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

\[ L(\{q_A\},\lambda) = \chi_{\mathrm{esp}}^2 - 2 \lambda \Big( \sum_{A=1}^M q_A -Q \Big) \]

where

\[\begin{equation*} \chi_{\mathrm{esp}}^2 = \sum_a \big[ V(\mathbf{r}_a) - \tilde{V}(\mathbf{r}_a) \big]^2 \end{equation*}\]

Finding the stationary point of this Lagrangian results in the following equation from which the ESP charges are obtained

\[\begin{equation*} \begin{pmatrix} A_{11} & \cdots & A_{1M} & 1 \\ \vdots & \ddots & \vdots & 1 \\ A_{M1} & \cdots & A_{MM} & 1 \\ 1 & \cdots & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} q_1\\ \vdots \\ q_M \\ \lambda \\ \end{pmatrix} = \begin{pmatrix} b_1\\ \vdots \\ b_M \\ Q \\ \end{pmatrix} \end{equation*}\]

where

\[\begin{equation*} A_{IJ} = \frac{1}{4\pi\varepsilon_0} \sum_{a} \frac{1}{r_{aI}r_{aJ}} ; \qquad r_{aI} = |\mathbf{r}_a - \mathbf{R}_I | \end{equation*}\]

and

\[\begin{equation*} b_I = \sum_{a} \frac{V(\mathbf{r}_a)}{r_{aI}} \end{equation*}\]

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.

Hide 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