Support for simulations with multiple small molecule ligands

[JSLJ23]

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