Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support for simulations with multiple small molecule ligands #338

Open
JSLJ23 opened this issue Jun 5, 2024 · 4 comments
Open

Support for simulations with multiple small molecule ligands #338

JSLJ23 opened this issue Jun 5, 2024 · 4 comments

Comments

@JSLJ23
Copy link

JSLJ23 commented Jun 5, 2024 •
edited
Loading

Overview

I am trying to run a simple NPT simulation with multiple ligands from a single SDF file parameterised with GAFF but I am running into ValueErrors for multiple definitions for an atom type.

May I know if there is a proper / better supported way for something like this to be done (i.e. running a simulation with multiple small molecules as rdkit/openff-Molecule objects)? Or if the current approach is sensible, are there any tweak I should make to allow for unique atom naming to circumvent this "Found multiple definitions for atom type" error?

Code

from openff.toolkit.topology import Molecule
from openmm import LangevinMiddleIntegrator
from openmm.app import PME, Modeller, Simulation, StateDataReporter
from openmm.app.dcdreporter import DCDReporter
from openmm.app.forcefield import ForceField
from openmm.openmm import MonteCarloBarostat, Platform
from openmm.unit import (
    angstrom,
    atmosphere,
    femtosecond,
    kelvin,
    kilojoule_per_mole,
    molar,
    nanometers,
    picosecond,
)
from openff.toolkit import Topology
from openmmforcefields.generators import GAFFTemplateGenerator
from rdkit import Chem

sdf_path = "./data/9_mols.sdf"

mols = []
off_mols = []

for mol in Chem.SDMolSupplier(str(sdf_path), removeHs=False):
    mols.append(mol)

# Modeller and force field with the first molecule
off_mol = Molecule.from_rdkit(mols[0], allow_undefined_stereo=True, hydrogens_are_explicit=True)
off_mol.generate_unique_atom_names()
mol_topology = off_mol.to_topology().to_openmm()
mol_positions = off_mol.to_topology().get_positions().to_openmm()
modeller = Modeller(mol_topology, mol_positions)
force_field = ForceField("amber14-all.xml", "amber14/tip3p.xml")
off_mols.append(off_mol)

# All other molecules
for i, mol in enumerate(mols[1:]):
    off_mol = Molecule.from_rdkit(mol, allow_undefined_stereo=True, hydrogens_are_explicit=True)
    off_mol.generate_unique_atom_names()
    mol_topology = off_mol.to_topology().to_openmm()
    mol_positions = off_mol.to_topology().get_positions().to_openmm()
    modeller.add(mol_topology, mol_positions)
    off_mols.append(off_mol)

mol_ff = GAFFTemplateGenerator(molecules=off_mols)
force_field.registerTemplateGenerator(mol_ff.generator)

modeller.addSolvent(force_field, padding=1.0 * nanometers, ionicStrength=0.15 * molar)
system = force_field.createSystem(modeller.topology, nonbondedMethod=PME)

solvated_system_full_path = "./combined_solvated.pdb"
PDBFile.writeFile(modeller.topology, modeller.positions, open(solvated_system_full_path, "w"))

Full traceback:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 1
----> 1 modeller.addSolvent(force_field, padding=1.0 * nanometers, ionicStrength=0.15 * molar)
      2 system = force_field.createSystem(modeller.topology, nonbondedMethod=PME)

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/modeller.py:519](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/modeller.py#line=518), in Modeller.addSolvent(self, forcefield, model, boxSize, boxVectors, padding, numAdded, boxShape, positiveIon, negativeIon, ionicStrength, neutralize, residueTemplates)
    515         raise ValueError('Neither the box size, box vectors, nor padding was specified, and the Topology does not define unit cell dimensions')
    517 # Have the ForceField build a System for the solute from which we can determine van der Waals radii.
--> 519 system = forcefield.createSystem(self.topology, residueTemplates=residueTemplates)
    520 nonbonded = None
    521 for i in range(system.getNumForces()):

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py:1247](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py#line=1246), in ForceField.createSystem(self, topology, nonbondedMethod, nonbondedCutoff, constraints, rigidWater, removeCMMotion, hydrogenMass, residueTemplates, ignoreExternalBonds, switchDistance, flexibleConstraints, drudeMass, **args)
   1243 rigidResidue = [False]*topology.getNumResidues()
   1245 # Find the template matching each residue and assign atom types.
-> 1247 templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds)
   1248 for res in topology.residues():
   1249     if res.name == 'HOH':
   1250         # Determine whether this should be a rigid water.

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py:1452](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py#line=1451), in ForceField._matchAllResiduesToTemplates(self, data, topology, residueTemplates, ignoreExternalBonds, ignoreExtraParticles, recordParameters)
   1449 if matches is None:
   1450     # Try all generators.
   1451     for generator in self._templateGenerators:
-> 1452         if generator(self, res):
   1453             # This generator has registered a new residue template that should match.
   1454             [template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
   1455             if matches is None:
   1456                 # Something went wrong because the generated template does not match the residue signature.

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmmforcefields/generators/template_generators.py:547](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmmforcefields/generators/template_generators.py#line=546), in GAFFTemplateGenerator.generator(self, forcefield, residue)
    544     # Note that we've loaded the GAFF parameters
    545     self._gaff_parameters_loaded[forcefield] = True
--> 547 return super().generator(forcefield, residue)

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmmforcefields/generators/template_generators.py:323](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmmforcefields/generators/template_generators.py#line=322), in SmallMoleculeTemplateGenerator.generator(self, forcefield, residue)
    320         outfile.write(ffxml_contents)
    322 # Add the parameters and residue definition
--> 323 forcefield.loadFile(StringIO(ffxml_contents))
    324 # If a cache is specified, add this molecule
    325 if self._cache is not None:

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py:288](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py#line=287), in ForceField.loadFile(self, files, resname_prefix)
    286     if tree.getroot().find('AtomTypes') is not None:
    287         for type in tree.getroot().find('AtomTypes').findall('Type'):
--> 288             self.registerAtomType(type.attrib)
    290 # Load the residue templates.
    292 for tree in trees:

File [~/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py:439](http://localhost:8888/home/joshua/mambaforge/envs/chem_py311/lib/python3.11/site-packages/openmm/app/forcefield.py#line=438), in ForceField.registerAtomType(self, parameters)
    437 name = parameters['name']
    438 if name in self._atomTypes:
--> 439     raise ValueError('Found multiple definitions for atom type: '+name)
    440 atomClass = parameters['class']
    441 mass = _convertParameterToNumber(parameters['mass'])

ValueError: Found multiple definitions for atom type: c5
@JSLJ23
Copy link
Author

JSLJ23 commented Jun 5, 2024

9_mols.sdf.zip

@peastman peastman transferred this issue from openmm/openmm Jun 5, 2024
@peastman
Copy link
Member

peastman commented Jun 5, 2024

I transferred this issue to openmmforcefields, since that's where the problem is. It seems the template generator is creating multiple atom types with the same name.

@mattwthompson
Copy link
Collaborator

openmm/openmm#4531 might fix this

@peastman
Copy link
Member

peastman commented Jun 5, 2024

Good point. @JSLJ23 can you test with the latest development code for OpenMM and see if it works?

@peastman peastman mentioned this issue Jun 21, 2024
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@peastman @mattwthompson @JSLJ23