RESP charges - eChem book

The restrained electrostatic potential (RESP) charge model Bayly et al. (1993)Cieplak et al. (1995) represents an improvement to the Merz–Kollman (MK) scheme as the ESP figure-of-merit is rather insensitive to variations in charges of atoms buried inside the molecule.

Dependence of figure-of-merit, \chi^2_\mathrm{esp}, with respect to variations in atomic charges. Four separate atoms are here considered. — Figure 1: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

\chi_{\mathrm{rstr}}^2 = \alpha \sum_{A=1}^M \big[ (q_A^2 + \beta^2)^{1/2} - \beta \big]

(1)

so that the diagonal matrix elements of the $A$ -matrix in MK scheme become equal to

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

(2)

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.

Let us determine the RESP charges for methanol at the Hartree–Fock/6-31G(d) level of theory. We will require the partial charges of the hydrogen atoms in the methyl group to be identical.

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

import py3Dmol as p3d

viewer = p3d.view(width=300, height=200)

viewer.addModel(methanol_xyz, "xyz")
viewer.setStyle({"stick": {}})

viewer.show()

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 the RESP charges with the compute method of the RespChargesDriver class.

resp_drv = vlx.RespChargesDriver()

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

resp_charges = resp_drv.compute(molecule, basis, scf_results)

                                                                                                                          
                                                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 9 iterations.                                       
                                                                                                                          
                                        No. | Atom |  Constraints | Charges (a.u.)                                        
                                       --------------------------------------------                                       
                                          1     H                      0.075643                                           
                                          2     C                      0.117258                                           
                                          3     H                      0.013902                                           
                                          4     H                      0.013047                                           
                                          5     O                     -0.639004                                           
                                          6     H                      0.419154                                           
                                       --------------------------------------------                                       
                                               Total Charge  :  0.000000                                                  
                                                                                                                          
                                                       Fit Quality                                                        
                                                      -------------                                                       
                                       Relative Root-Mean-Square Error  :  0.139861                                       
                                                                                                                          
                                                                                                                          
                                                     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 4 iterations.                                       
                                                                                                                          
                                    No. | Atom | Frozen | Constraints | Charges (a.u.)                                    
                                   ----------------------------------------------------                                   
                                      1     H       No                     0.033747                                       
                                      2     C       No                     0.118610                                       
                                      3     H       No         1           0.033747                                       
                                      4     H       No         1           0.033747                                       
                                      5     O      Yes                    -0.639004                                       
                                      6     H      Yes                     0.419154                                       
                                   ----------------------------------------------------                                   
                                               Total Charge  :  0.000000                                                  
                                                                                                                          
                                                       Fit Quality                                                        
                                                      -------------                                                       
                                       Relative Root-Mean-Square Error  :  0.203249                                       
                                                                                                                          
                                       Reference:                                                                         
                                       J. Phys. Chem. 1993, 97, 10269-10280.

Third, we print out the results.

print("Atom     RESP charge")
print(20 * "-")

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

print(20 * "-")

print(f"Total: {resp_charges.sum() : 13.6f}")

Atom     RESP charge
--------------------
H           0.033747
C           0.118610
H           0.033747
H           0.033747
O          -0.639004
H           0.419154
--------------------
Total:      0.000000

Compared to the ESP results, it is noted that the partial charge of the buried carbon atom is here lower (0.12 a.u.) as to be expected from the design of the RESP scheme.

References¶

Bayly, C. I., Cieplak, P., Cornell, W., & Kollman, P. A. (1993). A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem., 97, 10269–10280. 10.1021/j100142a004
Cieplak, P., Cornell, W. D., Bayly, C., & Kollman, P. A. (1995). Application of the multimolecule and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins. J. Comput. Chem., 16, 1357–1377. 10.1002/jcc.540161106