ESP charges#

Since there is no unique definition for partial charges and no corresponding physical observable, they can be assigned in several ways, such as being derived from the quantum mechanical electrostatic potential

\[\begin{equation*} V(\boldsymbol{r}) = \sum_{I} \frac{Z_I e}{4\pi\varepsilon_0 |\boldsymbol{r}-\mathrm{\textbf{R}}_I|} - e \sum_{\mu,\nu} D_{\mu\nu} \int \frac{ \phi_\mu^*(\boldsymbol{r}')\phi_\nu(\boldsymbol{r}') }{ 4\pi\varepsilon_0 |\boldsymbol{r}-\boldsymbol{r}'| } d^3\boldsymbol{r}' \end{equation*}\]

that can be replaced with a potential caused by the partial charges:

\[\begin{equation*} \widetilde{V}(\boldsymbol{r}) = \sum_{I} \frac{ q_I }{ 4\pi\varepsilon_0 |\boldsymbol{r}-\textbf{R}_I| } \end{equation*}\]

The Merz–Kollman scheme [BMK90, SK84] minimizes the squared norm difference between these two quantities 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 – the grid points are distributed on successive layers of scaled van der Waals surfaces. This measure is referred to as the figure-of-merit

\[\begin{equation*} \chi_{\mathrm{esp}}^2 = \sum_a \bigl(V(\boldsymbol{r}_a) - \widetilde{V}(\boldsymbol{r}_a)\bigl)^2 \end{equation*}\]

The resulting electrostatic potential (ESP) charges are obtained by solving the equation

\[\begin{equation*} \mathrm{\textbf{A}} \, \mathrm{\textbf{q}} = \mathrm{\textbf{b}} \end{equation*}\]

where

\[\begin{equation*} A_{JI} = \frac{1}{4\pi\varepsilon_0} \sum_{a} \frac{1}{r_{aI}r_{aJ}} \end{equation*}\]

and

\[\begin{equation*} b_J = \sum_{a} \frac{V_a}{r_{aJ}} \end{equation*}\]
import py3Dmol as p3d

methanol_xyz = """
  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
"""

v = p3d.view(width=400, height=200)

v.addModel("6\n" + methanol_xyz, "xyz")
v.setStyle({'stick':{}})
v.zoomTo()
v.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_str(methanol_xyz, units = 'angstrom')
basis = vlx.MolecularBasis.read(molecule, "def2-SVP")

print(basis.get_string(molecule))

A calculation of ESP charges is performed by invoking the RespChargesDriver driver class.

esp_drv = vlx.RespChargesDriver()

esp_drv.print_keywords()
                     ================================================================================                     
                       @resp charges                                                                                      
                     --------------------------------------------------------------------------------                     
                       number_layers       integer        number of layers of scaled vdW surfaces                         
                       density             float          density of points in each layer                                 
                       restrain_hydrogen   boolean        restrain hydrogen atoms                                         
                       weak_restraint      float          strength of restraint in 1st RESP stage                         
                       strong_restraint    float          strength of restraint in 2nd RESP stage                         
                       max_iter            integer        maximum iterations in RESP fit                                  
                       threshold           float          convergence threshold of RESP fit                               
                       xyz_file            string         xyz file containing the conformers                              
                       net_charge          float          net charge of the molecule                                      
                       multiplicity        integer        spin multiplicity of the molecule                               
                       weights             sequence       weight factors of conformers                                    
                       energies            sequence       energies of conformers                                          
                       temperature         float          temperature for Boltzmann weight factor                         
                       fitting_points      multi-line     points on which ESP charges are computed                        
                     ================================================================================                     

Let us modify a few of the default settings and calculate the ESP charges of methanol with use of the compute method.

esp_drv.update_settings({
    'number_layers': 5,
    'density': 10.0,
})

esp_charges = esp_drv.compute(molecule, basis, 'esp')
print("Atom   Charge (a.u.)")
print("-"*20)

for label, charge in zip(molecule.get_labels(), esp_charges):
    print(f"{label :s} {charge : 18.6f}")

print("-"*20)
print(f"Total: {esp_charges.sum() : 13.6f}")
Atom   Charge (a.u.)
--------------------
H           0.038025
C           0.246240
H          -0.018417
H          -0.018484
O          -0.655265
H           0.407901
--------------------
Total:      0.000000