RESP charges calculation using Psikit (Psi4 + RdKit)

The importance of RESP charges and how to calculate them with open source tools

Melchor Sanchez-Martinez
13 September 2021

Restrained electrostatic potential atomic partial (RESP) charges calculation is nowadays a common practice to prepare ligands for MD simulations. It is a highly useful and widely used method of assigning partial charges to molecules for simulations. Molecular Dynamics, MD, force-fields (FFs) are usually prepared to be used with proteins, to simulate proteins. Some of the most popular FFs are Amber, Charmm or Gromos that can be used within different simulation programs like Gromcas, NAMD or OpenMM, that are prepared to work with different FFs, or within more speciffic tools like Amber or Gromos, that were developed to work with an specific family of FFs althougth there are ways to incorporate other FFs. There are also specific force-fields prepared to simulate DNA, like the parmbsc01, and even organic molecules, like GAFF/GAFF2 or the newly developed Parsley.

Although there are specific parameterized FFs to simulate organic molecules and there are different tools such as antechamber, ACPYPE or ATB that allow parameterizing organic molecules to be simulated with different FFs, the resulting atomic charges are not always accurate / correct. Some of the tools like antechamber, open source tool from AmberTools, recommend you to use Quantum Mechanics-based (QM) methods to obtain more accurate partial charges. Within these methods there are AM1-BCC or RESP charges. AM1-BCC generates charges that emulate the electrostatic potential (ESP) HF/6-31G* of a molecule (more information can be found in the original article describing the method). The so-called ESP charges are derived from fitting a classical Coulomb model to the molecular electrostatic potentials QM. Although simulation methods using ESP charges are used to reproduce the geometries of hydrogen-bonded complexes, the strength of these interactions is overestimated. To overcome this problem, it is possible to use a restriction function during the adjustment of the partial charges to the electrostatic potentials, which attenuate the obtained charges. These are the so-called RESP charges, which in principle will produce more accurate results.

RESP charges calculation using Python, the Psi4 library. Code available in GitHub; Described in the original publication

I’m not sure if these kind of methods were first implemented for AMBER or not, but since they have been implemented in antechamber for a long time, they are regularly used for AMBER simulations, or Amber FF based simulations, of protein-ligand complexes from years ago. Considering that it is claimed that AmberFFs are those that most correctly reproduce the behavior of proteins and that they are, if not the most, one of the most widely used FFs, the use of RESP charges is quite useful and common. However, its use within antechamber is linked to Gaussian. And to use Gaussian you need a license; a fact that has limited its use. Fortunately, today, there are several options for calculating these charges with open source tools.

Actually RESP charges can be calculated with different open source software tools like CP2K or psi4 (in this case not directly in psi4 but through the resp library from cdsgroup) and even web servers like q4mdforcefields.org. I especially like psi4 because it is open source and a Python library so you can easily manage QM calculations within Python scripts. This possibility has allowed the emergence of libraries that merge Psi4 with other popular cheminformatics tools such as RdKit, obtaining Psikit. In this blog post you can find a clear example of how to use Psikit for Calculate RESP charges, also in the Psikit repository you can find an example. There are also proposed improved versions of the RESP charges, RESP2 and W-RESP. RESP2 is described in a recent paper and the code is available in github. However, it is based on OpenEye-toolkit which is licensed, it is not open source. The idea behind RESP2 does not seem too complicated, so maybe it could be implemented soon as an open source tool. W-RESP, also described in a recent paper, is based on Gaussian RESP software. They implement a new restriction function on RESP. However in the biorxiv paper (not checked in the JCTC paper that appeared online few months ago) there is not code provided and the new fitting function is not described. So actually the use of RESP charges seems to be the best option.

metformin='CN(C)C(=N)NC(=N)N'
pk.read_from_smiles(metformin)
pk.optimize('HF/6-31g(d)')
respCharges=pk.calc_resp_charges()
res=defaultdict(list)
for i, atom in enumerate(pk.mol.GetAtoms()):
    res['SYMBOL'].append(atom.GetSymbol())
    res['RESP'].append(respCharges[i])
molecule_df = pd.DataFrame(res)

pk.read_from_molfile('metformin.mol')
pk.optimize('HF/6-31g(d)')
respCharges=pk.calc_resp_charges()
res=defaultdict(list)
for i, atom in enumerate(pk.mol.GetAtoms()):
    res['SYMBOL'].append(atom.GetSymbol())
    res['RESP'].append(respCharges[i])
molecule_df = pd.DataFrame(res)

Example of RESP charges calculation starting from Smiles or a MOL file using Psikit. Code inspired by Iwatobipen’s blog post

Now we have more clear why RESP charges are interesting and how we can use it, but what can be the “real” applications of RESP charges? So for instance CADD projects. In a drug discovery project, the use of RESP charges could be important, for example, to perform MD simulations of docking complexes with the aim of incorporating flexibility after rigid docking to test the stability of the complex and obtain a more accurate description of the binding mode due to the incorporation of the so-called induced fit events (more information on post-processing of docking poses with MD can be found here or here) and / or to perform more accurate free (binding) energy calculations of the complex. In practice, what is done is to optimize the geometry of the docked molecule, calculate the RESP charges of the optimized conformation and add them together with the Cartesian coordinates, xyz, from the docking calculation, to a new SDF/MOL/MOL2 file. Then we have several options to use this file but one that I like it is to input this file to antechamber obtaining a topology file incorporating RESP charges. An interesting option would be to make antechamber recognize psi4 RESP calculation output, this will simplify the file manipulation, but by now I have not try it. Maybe it’s already done and I am not aware of it. If it is the case, please let me know!

So, let’s try Psikit! Then add the RESP charges to you protein-ligand simulations and check if your results improve! Here you can find an example that uses Psikit to calculate RESP charges that then are reflected into a parameterized file coming from antechamber that finally is converted to Gromacs topology.