Restrain equation is probably incorrect

[UnWaDo]

Bayly's 1993 article on RESP gives hyperbolic restrain as shown below (eq. 10). In resp/espfit.py there is a function restraint, which adds resp_a/np.sqrt(q[i]**2 + resp_b**2) * num_conformers to unrestrained matrix. As far as I understand, it misses the multiplication on charge q[i], so the formula should be q[i]*resp_a/np.sqrt(q[i]**2 + resp_b**2) * num_conformers to match one given in article. Could you, please, justify this: am I misunderstanding something or is there an error in code?

[alenaizan]

Thank you for your question. I believe q_j is not written because it is just the vector q in the A q = B matrix equation. Maybe the notation is not very clear.

You can verify that the code here is correct. The code agrees with the RESP implementation written by Bayly. I just downloaded the source code of Ambertools2022 from here: https://ambermd.org/GetAmber.php#ambertools. You can find the RESP code in the file: AmberTools/src/etc/resp.F. In the file you can find the following lines:

          qwtval(i) = qwt/sqrt(qcal(i)*qcal(i) + 0.01d0)
          a(i,i)= a(i,i) + qwtval(i)

Furthermore, using the script below you can verify that the implementation here agrees with the original RESP implementation:

import psi4
import resp
import numpy as np

mol = psi4.geometry("""
 O    0.00000001   0.06417294   0.00000000
 H   -0.75406643  -0.50923541   0.00000000
 H    0.75406629  -0.50923551   0.00000000
""")
mol.update_geometry()
mol.print_out_in_bohr()

options = {'VDW_SCALE_FACTORS'  : [1.4, 1.6, 1.8, 2.0],
           'VDW_POINT_DENSITY'  : 1.0,
           'RESP_A'             : 0.01,
           'RESP_B'             : 0.1,
           }

# Call for first stage fit
charges1 = resp.resp([mol], options)
print('Electrostatic Potential Charges')
print(charges1[0])
print('Restrained Electrostatic Potential Charges')
print(charges1[1])

grid = np.loadtxt("1_default_grid.dat")
grid /= psi4.constants.bohr2angstroms
esp = np.loadtxt("1_default_grid_esp.dat")

resp_file = open("espot_m1", "w")
resp_file.write("    %i  %i    0          MOLECULE\n" %(mol.natom(), len(esp)))
for i in range(3):
    x = mol.x(i)
    y = mol.y(i)
    z = mol.z(i)
    resp_file.write("%16s%16.7e%16.7e%16.7e\n" %("", x, y, z))

for i in range(len(esp)):
    resp_file.write("%16.7e%16.7e%16.7e%16.7e\n" %(esp[i], grid[i, 0], grid[i, 1], grid[i, 2]))

resp_file.write("\n\n")
resp_file.close()

The script generates the following RESP charges for water: -0.77531135 0.38638091 0.38893044. It also rewrites the grid data and electrostatic potential using the original RESP format. You can run the original RESP code using the following input file (name it input_file):

 &cntrl
  ioutopt=1, iqopt=1, nmol=1, ihfree=1, irstrnt=1, qwt= 0.0100 
 &end 
  1.0
 MOLECULE 
    0    3          
    8    0          
    1    0          
    1    0          

and the following command: resp -i input_file -o out_resp -p punch_resp -t qout_resp -e espot_m1 -s esout_resp

You will get the following charges in the qout_resp file: -0.775311 0.386381 0.388930.

So this verifies that the current implementation is correct.

Let me know if you have other questions.

[UnWaDo]

It is rather strange, because in article there is A_jj = ... + dχ²/dq_j (eq. 13), while eq. 10 gives the form of the derivative (and it contains q_j)

I believe that your implementation is exactly same as Dr. Bayly's, but do you have any ideas, why is it so?

[alenaizan]

I do not know why the equations were written this way. Yes, it is a little confusing.