RESP charges#

The restrained electrostatic potential (RESP) charge model [BCCK93, CCBK95] represents an improvement to the Merz–Kollman scheme as the ESP figure-of-merit is rather insensitive to variations in charges of atoms buried inside the molecule.

../../_images/chi_square.svg

Fig. 34 Dependence of figure-of-merit, \(\chi^2_\mathrm{esp}\), with respect to variations in atomic charges. Four separate atoms are here considered.#

To avoid unphysically high magnitudes of the charges of interior atoms, a hyperbolic penalty function is added

\[\begin{equation*} \chi_{\mathrm{rstr}}^2 = \alpha \sum_I \bigl((q_I^2+\beta^2)^{1/2}-\beta\bigl) \end{equation*}\]

so that the diagonal matrix elements of the \(A\)-matrix in Merz–Kollman scheme become equal to

\[\begin{equation*} A_{JJ} = \frac{1}{4\pi\varepsilon_0} \sum_{a} \frac{1}{r_{aJ}^2} + \alpha \, (q_J^2+\beta^2)^{-1/2} \end{equation*}\]

with a dependency on the partial charge. Consequently, RESP charges are obtained by solving the matrix equation iteratively until the charges and Lagrange multipliers become self-consistent. In addition to that, the RESP charge model allows for the introduction of constraints on charges of equivalent atoms due to symmetry operations or bond rotations.

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 RESP charges is performed by invoking the RespChargesDriver driver class.

resp_drv = vlx.RespChargesDriver()

resp_drv.update_settings({
    'equal_charges': '1 = 3, 1 = 4'
})

resp_charges = resp_drv.compute(molecule, basis, 'resp')
                                  *** Warning: Recommended basis set 6-31G* is not used!                                  
                                                                                                                          
                                                                                                                          
                                            Self Consistent Field Driver Setup                                            
                                           ====================================                                           
                                                                                                                          
                   Wave Function Model             : Spin-Restricted Hartree-Fock                                         
                   Initial Guess Model             : Superposition of Atomic Densities                                    
                   Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace                     
                   Max. Number of Iterations       : 50                                                                   
                   Max. Number of Error Vectors    : 10                                                                   
                   Convergence Threshold           : 1.0e-06                                                              
                   ERI Screening Scheme            : Cauchy Schwarz + Density                                             
                   ERI Screening Mode              : Dynamic                                                              
                   ERI Screening Threshold         : 1.0e-12                                                              
                   Linear Dependence Threshold     : 1.0e-06                                                              
                                                                                                                          
* Info * Nuclear repulsion energy: 40.7905269056 a.u.                                                                     
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * SAD initial guess computed in 0.00 sec.                                                                          
                                                                                                                          
* Info * Starting Reduced Basis SCF calculation...                                                                        
* Info * ...done. SCF energy in reduced basis set: -114.888819617018 a.u. Time: 0.10 sec.                                 
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
                                                                                                                          
               Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change               
               --------------------------------------------------------------------------------------------               
                  1      -114.952364853481    0.0000000000      0.12293596      0.01016005      0.00000000                
                  2      -114.953993923559   -0.0016290701      0.02896702      0.00286098      0.04180420                
                  3      -114.954077690538   -0.0000837670      0.00649658      0.00067164      0.00868558                
                  4      -114.954084213326   -0.0000065228      0.00229285      0.00021230      0.00267146                
                  5      -114.954084968419   -0.0000007551      0.00061792      0.00006954      0.00116667                
                  6      -114.954085013328   -0.0000000449      0.00017914      0.00002270      0.00022484                
                  7      -114.954085018145   -0.0000000048      0.00002955      0.00000300      0.00008349                
                  8      -114.954085018328   -0.0000000002      0.00000431      0.00000034      0.00001603                
                  9      -114.954085018334   -0.0000000000      0.00000100      0.00000008      0.00000309                
                 10      -114.954085018334   -0.0000000000      0.00000017      0.00000002      0.00000078                
                                                                                                                          
* Info * Checkpoint written to file: veloxchem_esp_2022-04-14T19.52.59.scf.h5                                             
                                                                                                                          
* Info * SCF tensors written to file: veloxchem_esp_2022-04-14T19.52.59.scf.tensors.h5                                    
                                                                                                                          
               *** SCF converged in 10 iterations. Time: 1.34 sec.                                                        
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :     -114.9540850183 a.u.                                              
               Electronic Energy                  :     -155.7446119240 a.u.                                              
               Nuclear Repulsion Energy           :       40.7905269056 a.u.                                              
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000001716 a.u.                                              
                                                                                                                          
                                                                                                                          
               Ground State Information                                                                                   
               ------------------------                                                                                   
               Charge of Molecule            :  0.0                                                                       
               Multiplicity (2S+1)           :  1.0                                                                       
               Magnetic Quantum Number (M_S) :  0.0                                                                       
                                                                                                                          
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.69455 a.u.                                                                  
               (   2 C   1p+1:     0.16) (   2 C   1p0 :    -0.24) (   3 H   1s  :     0.15)                              
               (   4 H   1s  :     0.15) (   5 O   3s  :    -0.18) (   5 O   1p+1:    -0.25)                              
               (   5 O   1p0 :    -0.33) (   5 O   2p0 :    -0.17) (   6 H   1s  :     0.30)                              
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.62208 a.u.                                                                  
               (   2 C   1p-1:    -0.37) (   2 C   2p-1:    -0.15) (   3 H   1s  :    -0.24)                              
               (   4 H   1s  :     0.24) (   5 O   1p-1:    -0.35) (   5 O   2p-1:    -0.24)                              
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.59265 a.u.                                                                  
               (   1 H   1s  :     0.27) (   2 C   1p+1:     0.33) (   2 C   1p0 :     0.18)                              
               (   5 O   1p+1:    -0.37) (   5 O   1p0 :     0.21) (   5 O   2p+1:    -0.24)                              
               (   5 O   2p0 :     0.15)                                                                                  
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.49859 a.u.                                                                  
               (   1 H   1s  :     0.29) (   2 C   1p0 :     0.32) (   3 H   1s  :    -0.17)                              
               (   4 H   1s  :    -0.17) (   5 O   3s  :    -0.21) (   5 O   1p0 :    -0.38)                              
               (   5 O   2p0 :    -0.27) (   6 H   1s  :     0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.44480 a.u.                                                                  
               (   2 C   1p-1:    -0.24) (   3 H   1s  :    -0.24) (   4 H   1s  :     0.24)                              
               (   5 O   1p-1:     0.53) (   5 O   2p-1:     0.44)                                                        
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.17822 a.u.                                                                  
               (   1 H   2s  :     0.77) (   2 C   3s  :    -2.21) (   2 C   2p+1:    -0.41)                              
               (   2 C   2p0 :     0.20) (   3 H   2s  :     1.07) (   4 H   2s  :     1.07)                              
               (   5 O   3s  :    -0.51) (   5 O   2p+1:     0.19) (   5 O   2p0 :     0.22)                              
               (   6 H   2s  :     0.96)                                                                                  
                                                                                                                          
               Molecular Orbital No.  11:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.22492 a.u.                                                                  
               (   1 H   2s  :    -1.25) (   2 C   3s  :     2.06) (   2 C   2p+1:     0.19)                              
               (   2 C   2p0 :     0.20) (   3 H   2s  :    -0.84) (   4 H   2s  :    -0.84)                              
               (   5 O   3s  :    -0.90) (   5 O   1p0 :     0.15) (   5 O   2p0 :     0.31)                              
               (   6 H   2s  :     1.14)                                                                                  
                                                                                                                          
               Molecular Orbital No.  12:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.26581 a.u.                                                                  
               (   1 H   2s  :     2.00) (   2 C   3s  :     0.22) (   2 C   1p0 :    -0.29)                              
               (   2 C   2p+1:    -0.41) (   2 C   2p0 :    -1.30) (   3 H   2s  :    -1.07)                              
               (   4 H   2s  :    -1.03) (   5 O   3s  :    -0.26) (   5 O   2p0 :     0.21)                              
                                                                                                                          
               Molecular Orbital No.  13:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.26603 a.u.                                                                  
               (   2 C   1p-1:    -0.33) (   2 C   2p-1:    -1.36) (   3 H   2s  :     1.81)                              
               (   4 H   2s  :    -1.83) (   5 O   2p-1:     0.19)                                                        
                                                                                                                          
               Molecular Orbital No.  14:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.36619 a.u.                                                                  
               (   1 H   2s  :     0.41) (   2 C   3s  :    -1.24) (   2 C   1p+1:     0.16)                              
               (   2 C   2p+1:     1.63) (   2 C   2p0 :    -0.57) (   3 H   2s  :    -0.22)                              
               (   4 H   2s  :    -0.22) (   5 O   3s  :     1.05) (   5 O   1p+1:     0.27)                              
               (   5 O   2p+1:     0.87) (   6 H   2s  :     0.89)                                                        
                                                                                                                          
                                                                                                                          
                                                RESP Charges Driver Setup                                                 
                                               ===========================                                                
                                                                                                                          
                                         Number of Conformers         :  1                                                
                                         Number of Layers             :  4                                                
                                         Points per Square Angstrom   :  1.0                                              
                                         Total Number of Grid Points  :  420                                              
                                                                                                                          
                                                                                                                          
                                                     First Stage Fit                                                      
                                                    -----------------                                                     
                                                                                                                          
                                         Restraint Strength           :  0.0005                                           
                                         Restrained Hydrogens         :  No                                               
                                         Max. Number of Iterations    :  50                                               
                                         Convergence Threshold (a.u.) :  1e-06                                            
                                                                                                                          
                                      *** Charge fitting converged in 10 iterations.                                      
                                                                                                                          
                                        No. | Atom |  Constraints | Charges (a.u.)                                        
                                       --------------------------------------------                                       
                                          1     H                      0.060752                                           
                                          2     C                      0.159466                                           
                                          3     H                      0.004218                                           
                                          4     H                      0.003584                                           
                                          5     O                     -0.633671                                           
                                          6     H                      0.405651                                           
                                       --------------------------------------------                                       
                                               Total Charge  :  0.000000                                                  
                                                                                                                          
                                                       Fit Quality                                                        
                                                      -------------                                                       
                                       Relative Root-Mean-Square Error  :  0.148570                                       
                                                                                                                          
                                                                                                                          
                                                     Second Stage Fit                                                     
                                                    ------------------                                                    
                                                                                                                          
                                         Restraint Strength           :  0.001                                            
                                         Restrained Hydrogens         :  No                                               
                                         Max. Number of Iterations    :  50                                               
                                         Convergence Threshold (a.u.) :  1e-06                                            
                                                                                                                          
                                      *** Charge fitting converged in 3 iterations.                                       
                                                                                                                          
                                    No. | Atom | Frozen | Constraints | Charges (a.u.)                                    
                                   ----------------------------------------------------                                   
                                      1     H       No                     0.022789                                       
                                      2     C       No                     0.159653                                       
                                      3     H       No         1           0.022789                                       
                                      4     H       No         1           0.022789                                       
                                      5     O      Yes                    -0.633671                                       
                                      6     H      Yes                     0.405651                                       
                                   ----------------------------------------------------                                   
                                               Total Charge  :  0.000000                                                  
                                                                                                                          
                                                       Fit Quality                                                        
                                                      -------------                                                       
                                       Relative Root-Mean-Square Error  :  0.202017                                       
                                                                                                                          
                                       Reference:                                                                         
                                       J. Phys. Chem. 1993, 97, 10269-10280.                                              
                                                                                                                          
print("Atom   Charge (a.u.)")
print("-"*20)

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

print("-"*20)
print(f"Total: {resp_charges.sum() : 13.6f}")
Atom   Charge (a.u.)
--------------------
H           0.022789
C           0.159653
H           0.022789
H           0.022789
O          -0.633671
H           0.405651
--------------------
Total:      0.000000