ESP charges

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 at positions \(\mathbf{r}_a\) are distributed on successive layers of scaled van der Waals surfaces, and the default is to introduce four layers associated with scaling factors of 1.4, 1.6, 1.8, and 2.0.

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. Added to Lagrangian presented above, it also offers a means to add constraints to require charges of equivalent atoms to be equal.

We will use this implementation to determine the ESP charges of methanol at the Hartree–Fock/6-31G(d) level of theory. The hydrogens of the methyl group will be required to have identical charge.

Hide code cell source 
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
"""

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

import veloxchem as vlx

First, we determine the reference state of the system.

molecule = vlx.Molecule.read_xyz_string(methanol_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G*", ostream=None)

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_results = scf_drv.compute(molecule, basis)

Second, we calculate ESP charges with the compute method of the RespChargesDriver class.

esp_drv = vlx.RespChargesDriver()

esp_drv.equal_charges = "1=3, 1=4"

esp_charges = esp_drv.compute(molecule, basis, scf_results, "esp")
Hide code cell output 
                                                                                                                          
                                                RESP Charges Driver Setup                                                 
                                               ===========================                                                
                                                                                                                          
                                         Number of Conformers         :  1                                                
                                         Number of Layers             :  4                                                
                                         Points per Square Angstrom   :  1.0                                              
                                         Total Number of Grid Points  :  420                                              
                                                                                                                          
                                                 Merz-Kollman ESP Charges                                                 
                                                --------------------------                                                
                                                                                                                          
                                              No.   Atom      Charge (a.u.)                                               
                                             -------------------------------                                              
                                                1     H         0.023220                                                  
                                                2     C         0.148458                                                  
                                                3     H         0.023220                                                  
                                                4     H         0.023220                                                  
                                                5     O        -0.594785                                                  
                                                6     H         0.376669                                                  
                                             -------------------------------                                              
                                               Total Charge  :  0.000000                                                  
                                                                                                                          
                                                       Fit Quality                                                        
                                                      -------------                                                       
                                       Relative Root-Mean-Square Error  :  0.180034                                       

                                                                                                                          
                                       Reference:                                                                         
                                       J. Comput. Chem. 1984, 5, 129-145.

Third, we print out the results.

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.023220
C           0.148458
H           0.023220
H           0.023220
O          -0.594785
H           0.376669
--------------------
Total:      0.000000