OpenFF with RESP or AM1BCC charges

[pradipchm]

Hi,
I am trying to generate gromacs top file with the resp charges with the following code. But the generated top file from this code is not working properly, bonds are breaking and having a deformed structure. Please have a look at my code if anything is wrong.

import numpy
from openff.toolkit.topology import Molecule

from openff.recharge.charges.library import (
    LibraryChargeCollection,
    LibraryChargeGenerator,
)
from openff.recharge.charges.resp import generate_resp_charge_parameter
from openff.recharge.charges.resp.solvers import IterativeSolver
from openff.recharge.esp import ESPSettings
from openff.recharge.esp.psi4 import Psi4ESPGenerator
from openff.recharge.esp.storage import MoleculeESPRecord
from openff.recharge.grids import MSKGridSettings
from openff.recharge.utilities.molecule import extract_conformers

from rdkit import Chem
from openff.interchange.components.interchange import Interchange
from openff.toolkit import ForceField, Molecule, Topology
from openmm import NonbondedForce, unit


import os
os.environ["PSI_SCRATCH"] = "/storage/nas_scr/ps0579/PSIscri"

qc_data_settings = ESPSettings(
    method="HF", basis="6-31G*", grid_settings=MSKGridSettings()
)

smiles = "Cc1cc(c(s1)Nc2ccccc2[N+](=O)[O-])C#N"
rd_mol = Chem.MolFromSmiles(smiles)
rd_mol = Chem.AddHs(rd_mol)

# Convert to OpenFF Molecule
molecule = Molecule.from_rdkit(rd_mol)

molecule.generate_conformers(n_conformers=1)

[input_conformer] = extract_conformers(molecule)
conformer, grid, esp, electric_field = Psi4ESPGenerator.generate(
    molecule, input_conformer, qc_data_settings, minimize=True, n_threads=36
)

qc_data_record = MoleculeESPRecord.from_molecule(
    molecule, conformer, grid, esp, None, qc_data_settings
)


resp_solver = IterativeSolver()
resp_charge_parameter = generate_resp_charge_parameter(
        [qc_data_record], resp_solver
    )

resp_charges = LibraryChargeGenerator.generate(
        molecule, LibraryChargeCollection(parameters=[resp_charge_parameter])
    )

print(f"RESP SMILES         : {resp_charge_parameter.smiles}")
print(f"RESP VALUES (UNIQUE): {resp_charge_parameter.value}")
print("")
print(f"RESP CHARGES        : {numpy.round(resp_charges, 12)}")

molecule.partial_charges = resp_charges.reshape(-1) * unit.elementary_charge

forcefield = ForceField("openff_unconstrained-2.2.0.offxml")
topology = Topology.from_pdb(
    "ROY.pdb",
    unique_molecules=[molecule],
)

system = forcefield.create_openmm_system(topology, charge_from_molecules=[molecule])

# Inspect the openmm.System object to see which charges
# actually ended up on each particle
for force in system.getForces():
    if type(force) == NonbondedForce:
        for i in range(force.getNumParticles()):
            print(force.getParticleParameters(i)[0])


interchange = Interchange.from_smirnoff(force_field=forcefield, topology=topology, charge_from_molecules=[molecule])
interchange.to_top("system_RESP_HF.top")
interchange.to_gro("system_RESP_HF.gro")

I have used this one for am1bcc which looks working fine.

from openff.toolkit import ForceField, Molecule, Topology
molecule = Molecule.from_smiles("Cc1cc(c(s1)Nc2ccccc2[N+](=O)[O-])C#N")
molecule.assign_partial_charges(partial_charge_method="am1bcc")


# Load the topology from a PDB file and `Molecule` objects
topology = Topology.from_pdb(
    "ROY.pdb",
    unique_molecules=[molecule],
)

topology
forcefield = ForceField("openff-2.2.0.offxml")
forcefield
omm_system = forcefield.create_openmm_system(topology)
omm_system
from openff.interchange import Interchange

interchange = Interchange.from_smirnoff(
    force_field=forcefield,
    topology=topology,
    charge_from_molecules=[molecule],
)
interchange
# Export GROMACS files.
interchange.to_top("system.top")
interchange.to_gro("system.gro")

Also I am bit confused about constrained vs unconstrained version of openff, anyway I am using LINCS while running gromacs.

[mattwthompson]

How are you running the GROMACS simulations (+ what files are you using)? You're using the constrained version of Sage¹ in that snippet, which includes constraints between all hydrogens and the heavy atoms they're associated with, so you need to also include constraints = h-bonds or similar in the .mdp file² for this to make sense. I think that without that setting, there's nothing keeping hydrogens attached.

Alternatively you could use the un-constrained version by loading ForceField("openff_unconstrained-2.2.0.offxml").

Footnotes

1. Search "constrain" here: https://docs.openforcefield.org/en/latest/faq.html#smirnoff-force-fields ↩

2. This may change in a future release, see https://github.com/openforcefield/openff-interchange/issues/1180 ↩

[mattwthompson]

No, unfortunately you're on your own for generating RESP charges, if that's what you're trying to use here, since OpenFF provides AM1-BCC and not RESP charges

[pradipchm]

Thanks. As it correctly gets the charges for the other molecule. So while calculating openff recharge if you can print psi4 input and output it would be nice to understand the error. I also put a similar issue related to openff recharge, while calculating openff for large molecule it shows an error about optimization failure.

[mattwthompson]

ESP calculations with large molecules are non-trivial, your best bet in figuring out Psi4 issues is the Psi4 documentation and/or forums at https://forum.psicode.org/

[pradipchm]

Thanks. I will look into it. What I am trying to say here, while running openff recharge, if you could generate the psi4 input and output file (now I can't see exact psi4 input and output), that will be helpful to figure out any issue related to psi4 calculations. Hope you got my point.

[pradipchm]

Just a small update. I separately calculated RESP using Psi4, which looks correct. It might be a reasonable idea to have two separate steps: charge calculations followed by generating the topology.

from psikit import Psikit
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from IPython.display import SVG
from rdkit.Chem.Draw import rdMolDraw2D
from collections import defaultdict
import pandas as pd
from rdkit.Chem.rdchem import AtomPDBResidueInfo



pk = Psikit(threads=36, memory=20)

molecule='Cc1cc(c(s1)Nc2ccccc2[N+](=O)[O-])C#N'
pk.read_from_smiles(molecule)
pk.optimize(basis_sets= "HF/6-31g(d)", maxiter=200)
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])


mol = Chem.MolFromSmiles(molecule)
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol, AllChem.ETKDG())

element_counts = defaultdict(int)

for atom in mol.GetAtoms():
    symbol = atom.GetSymbol()
    element_counts[symbol] += 1
    atom_name = f"{symbol}{element_counts[symbol]:<2}"  # Left-align to ensure 4-character width
    pdb_info = AtomPDBResidueInfo(atom_name)
    atom.SetMonomerInfo(pdb_info)

atom_data = []

for i, atom in enumerate(mol.GetAtoms()):
    atom_name = atom.GetPDBResidueInfo().GetName().strip()
    charge = respCharges[i]
    atom_data.append([atom_name, charge])


charges_df = pd.DataFrame(atom_data, columns=['Atom', 'RESP_Charge'])
charges_df.to_csv('molecule_charges.csv', sep='\t', index=False, float_format='%.6f')

Chem.MolToPDBFile(mol, 'molecule.pdb')

the output looks like

Electrostatic potential parameters

Grid information (see 1_default_grid.dat in Angstrom)
    Number of grid points:            1343

Quantum electrostatic potential (see 1_default_grid_esp.dat)

Constraints
    Charge constraints
        Total charge of   0.02985592 on the set   2
        Total charge of  -0.26162170 on the set   3
        Total charge of   0.02887014 on the set   4
        Total charge of  -0.03515440 on the set   5
        Total charge of  -0.01518896 on the set   6
        Total charge of  -0.17469074 on the set   7
        Total charge of   0.07790033 on the set   8
        Total charge of  -0.14164940 on the set   9
        Total charge of  -0.11735707 on the set  10
        Total charge of  -0.17514554 on the set  11
        Total charge of  -0.13400900 on the set  12
        Total charge of  -0.00490757 on the set  13
        Total charge of   0.68585148 on the set  14
        Total charge of  -0.46801572 on the set  15
        Total charge of  -0.41592542 on the set  16
        Total charge of   0.34951689 on the set  17
        Total charge of  -0.45731755 on the set  18
        Total charge of   0.20767427 on the set  22
        Total charge of   0.30359798 on the set  23
        Total charge of   0.13487416 on the set  24
        Total charge of   0.16410663 on the set  25
        Total charge of   0.15021367 on the set  26
        Total charge of   0.17399682 on the set  27
    Equality constraints
        Equal charges on atoms
                                 19  20  21

Restraint
    Hyperbolic restraint to a charge of zero
    Hydrogen atoms are not restrained
    resp_a:                           0.0010
    resp_b:                           0.1000

Fit

Electrostatic Potential Charges
   Center  Symbol       ESP      RESP
      1    C      -0.17435181 -0.15955093
      2    C       0.02985592  0.02985592
      3    C      -0.26162170 -0.26162170
      4    C       0.02887014  0.02887014
      5    C      -0.03515440 -0.03515440
      6    S      -0.01518896 -0.01518896
      7    N      -0.17469074 -0.17469074
      8    C       0.07790033  0.07790033
      9    C      -0.14164940 -0.14164940
     10    C      -0.11735707 -0.11735707
     11    C      -0.17514554 -0.17514554
     12    C      -0.13400900 -0.13400900
     13    C      -0.00490757 -0.00490757
     14    N       0.68585148  0.68585148
     15    O      -0.46801572 -0.46801572
     16    O      -0.41592542 -0.41592542
     17    C       0.34951689  0.34951689
     18    N      -0.45731755 -0.45731755
     19    H       0.08962552  0.08469190
     20    H       0.08962552  0.08469190
     21    H       0.08962552  0.08469190
     22    H       0.20767427  0.20767427
     23    H       0.30359798  0.30359798
     24    H       0.13487416  0.13487416
     25    H       0.16410663  0.16410663
     26    H       0.15021367  0.15021367
     27    H       0.17399682  0.17399682
Total Charge:      0.00000000  0.00000000

[j-wags]

I'm not confident (and I don't have time to try and reproduce this, apologies) but this might be something having to do with the different way the mols are created in the two code paths. In the AM1BCC path it's

molecule = Molecule.from_smiles("Cc1cc(c(s1)Nc2ccccc2[N+](=O)[O-])C#N")

but in the RESP path it's

smiles = "Cc1cc(c(s1)Nc2ccccc2[N+](=O)[O-])C#N"
rd_mol = Chem.MolFromSmiles(smiles)
rd_mol = Chem.AddHs(rd_mol)

# Convert to OpenFF Molecule
molecule = Molecule.from_rdkit(rd_mol)

There are lots of "gotchas" surrounding molecule sanitization that might cause these to produce different numbers of Hs (or other subtly different things about the mols), and the Molecule.from_smiles route is the one we'd recommend (since then we handle the hydrogen addition/sanitization internally in OpenFF code in a consistent way). So maybe try using Molecule.from_smiles in the RESP code path and see if that helps?

[pradipchm]

I have used Molecule.from_smiles instead of rdmol, but the result does not change and still there are some wrong bonds.

[j-wags]

Darn, I guess that's not the problem. Thanks for checking this out!