SCHWED-882: included externally parametrized ligands in energy minimization

modelling/pymod/_pipeline.py

@@ -11,7 +11,7 @@ from _closegaps import *
 from _ring_punches import *
 # external
 import ost
-from ost import mol
+from ost import mol, conop
 from ost.mol import mm
 import os
 
@@ -28,41 +28,139 @@ def _RemoveHydrogens(ent):
     for a in ha.atoms:
       edi.DeleteAtom(a.handle)
 
-def _SetupMmSimulation(mhandle):
-    '''Get mm simulation object for the model.'''
+def _AddHeuristicHydrogens(ent, ff):
+    '''Add hydrogens with mm.HeuristicHydrogenConstructor.'''
+    for res in ent.residues:
+        if not res.IsPeptideLinking():
+            bb = ff.GetBuildingBlock(res.name)
+            edi = ent.EditXCS()
+            h_constructor = mm.HeuristicHydrogenConstructor(bb)
+            h_constructor.ApplyOnResidue(res, edi)
+
+def _GetTopology(ent, settings, force_fields, add_heuristic_hydrogens=False):
+    '''Return topology or None if no topology could be created.
+    Note: if successful, this will update ent (adding hydrogens).
+    Set add_heuristic_hydrogens to True for ligands.
+    '''
+    # try all force fields in order
+    last_exception = None
+    for ff in force_fields:
+        settings.forcefield = ff
+        try:
+            # check if we need to add hydrogens heuristically
+            if add_heuristic_hydrogens:
+                _AddHeuristicHydrogens(ent, ff)
+            # ok now we try...
+            topo = mm.TopologyCreator.Create(ent, settings)
+        except Exception as ex:
+            # keep track of error and continue
+            last_exception = type(ex).__name__ + ": " + ex.message
+            continue
+        else:
+            # all good
+            return topo
+    # if we got here, nothing worked
+    ost.LogError("Could not create mm topology. Last exception: "\
+                 + last_exception)
+    return None
+
+def _AddLigands(ent, top, lig_ent, settings, force_fields):
+    '''Add ligands from lig_ent to topology top and update entity ent.'''
+    # connect them first
+    proc = conop.RuleBasedProcessor(conop.GetDefaultLib())
+    proc.Process(lig_ent)
+    cur_res = lig_ent.residues[0]
+    lig_num = 1
+    while cur_res.IsValid():
+        # setup connected components
+        cur_view = lig_ent.CreateEmptyView()
+        cur_view.AddResidue(cur_res, mol.INCLUDE_ATOMS)
+        cur_res = cur_res.next
+        while cur_res.IsValid() and mol.InSequence(cur_res.prev, cur_res):
+            cur_view.AddResidue(cur_res, mol.INCLUDE_ATOMS)
+            cur_res = cur_res.next
+        # try to add topology with special named chain
+        view_res_str = str([str(r) for r in cur_view.residues])
+        cur_ent = mol.CreateEntityFromView(cur_view, True)
+        edi = cur_ent.EditXCS()
+        edi.RenameChain(cur_ent.chains[0], '_' + str(lig_num))
+        lig_num += 1
+        cur_top = _GetTopology(cur_ent, settings, force_fields, True)
+        if cur_top is None:
+            ost.LogError("Failed to add ligands " + view_res_str + \
+                         " for energy minimization! Skipping...")
+        else:
+            # merge into main topology
+            cur_top.SetFudgeLJ(top.GetFudgeLJ())
+            cur_top.SetFudgeQQ(top.GetFudgeQQ())
+            top.Merge(ent, cur_top, cur_ent)
+
+def _SetupMmSimulation(model, force_fields):
+    '''Get mm simulation object for the model (incl. ligands in chain "_").
+    This tries to generate a topology for the protein and for each connected
+    component in the ligand chain separately by evaluating the force fields in
+    the same order as given. Ligands without force fields are skipped.
+    '''
     # get general settings 
     settings = mm.Settings()
     settings.integrator = mm.LangevinIntegrator(310, 1, 0.002)
     settings.init_temperature = 0
     settings.nonbonded_method = mm.NonbondedMethod.CutoffNonPeriodic
     settings.keep_ff_specific_naming = False
-    # use fast CPU platform by default
-    settings.platform = mm.Platform.CPU
-    if os.getenv('PM3_OPENMM_CPU_THREADS') is None:
-        settings.cpu_properties["CpuThreads"] = "1"
-    else:
-        settings.cpu_properties["CpuThreads"] = os.getenv('PM3_OPENMM_CPU_THREADS')
-    # prepare protein
-    model = mol.CreateEntityFromView(mhandle.model.Select("cname!='_'"), True)
+    
+    # prepare entity with protein
     _RemoveHydrogens(model)
-    # NOTE: both FF work, performance is similar, CHARMM slightly better
-    settings.forcefield = mm.LoadCHARMMForcefield()
-    #settings.forcefield = mm.LoadAMBERForcefield()
-
-    # TODO: add ligands
-    # TODO: check when platform problem appears when having multiple topologies!
+    ent = mol.CreateEntityFromView(model.Select("cname!='_'"), True)
+    top = _GetTopology(ent, settings, force_fields)
+    if top is None:
+        raise RuntimeError("Failed to setup protein for energy minimization!")
+    
+    # prepare ligands: we reprocess them to ensure connectivity
+    lig_ent = mol.CreateEntityFromView(model.Select("cname='_'"), True)
+    if lig_ent.residue_count > 0:
+        _AddLigands(ent, top, lig_ent, settings, force_fields)
 
     # finally set up the simulation
     try:
-        sim = mm.Simulation(model, settings)
+        # use fast CPU platform by default
+        settings.platform = mm.Platform.CPU
+        num_cpu_threads = os.getenv('PM3_OPENMM_CPU_THREADS')
+        if num_cpu_threads is None:
+            settings.cpu_properties["CpuThreads"] = "1"
+        else:
+            settings.cpu_properties["CpuThreads"] = num_cpu_threads
+        
+        # TODO: check when platform problem appears when having multiple 
+        #       topologies! -> should be here
+
+        sim = mm.Simulation(top, ent, settings)
     except RuntimeError as ex:
         ost.LogWarning("OpenMM failed to initialize with error: %s" % str(ex))
         # switch to "mm.Platform.Reference" as fallback
         settings.platform = mm.Platform.Reference
-        sim = mm.Simulation(model, settings)
+        sim = mm.Simulation(top, ent, settings)
         ost.LogWarning("Switched to slower reference platform of OpenMM!")
     return sim
 
+def _GetSimEntity(sim):
+    '''Get Entity from mm sim and reverse ligand chain naming.'''
+    ent = sim.GetEntity().Copy()
+    # merge ligand chains into _
+    ent_ed = ent.EditXCS()
+    chain_names = [ch.name for ch in ent.chains]
+    for chain_name in chain_names:
+        # all separate ligand chains start with _
+        if chain_name[0] == '_':
+            # add to chain _
+            if not ent.FindChain('_').IsValid():
+                ent_ed.InsertChain('_')
+            lig_ch = ent.FindChain('_')
+            cur_ch = ent.FindChain(chain_name)
+            for res in cur_ch.residues:
+                ent_ed.AppendResidue(lig_ch, res, deep=True)
+            # remove old chain
+            ent_ed.DeleteChain(cur_ch)
+    return ent
 ###############################################################################
 
 def BuildSidechains(mhandle, merge_distance=4, scorer=None, fragment_db=None,
@@ -134,7 +232,8 @@ def BuildSidechains(mhandle, merge_distance=4, scorer=None, fragment_db=None,
             
 
 def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
-                        max_iter_lbfgs=10):
+                        max_iter_lbfgs=10, use_amber_ff=False,
+                        extra_force_fields=list()):
     '''Minimize energy of final model using molecular mechanics.
 
     Uses :mod:`ost.mol.mm` to perform energy minimization.
@@ -164,14 +263,36 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
     :param max_iter_lbfgs: Max. number of iterations within LBFGS method
     :type max_iter_lbfgs:  :class:`int`
 
+    :param use_amber_ff: if True, use the AMBER force field instead of the def.
+                         CHARMM one (see :meth:`ost.mol.mm.LoadAMBERForcefield`
+                         and :meth:`ost.mol.mm.LoadCHARMMForcefield`).
+                         Both do a similarly good job without ligands (CHARMM
+                         slightly better), but you will want to be consistent
+                         with the optional force fields in `extra_force_fields`.
+    :type use_amber_ff:  :class:`bool`
+
+    :param extra_force_fields: Additional list of force fields to use if a 
+                               (ligand) residue cannot be parametrized with the
+                               default force field. The force fields are tried
+                               in the order as given and ligands without an
+                               existing parametrization are skipped.
+    :type extra_force_fields:  :class:`list` of :class:`ost.mol.Forcefield`.
+
     :return: The model including all oxygens as used in the minimizer.
     :rtype:  :class:`Entity <ost.mol.EntityHandle>`
     '''
     ost.LogInfo("Minimize energy.")
     # ignore LogInfo in stereochemical problems if output up to info done
     ignore_stereo_log = (ost.GetVerbosityLevel() == 3)
+
+    # setup force fields
+    if use_amber_ff:
+        force_fields = [mm.LoadAMBERForcefield()]
+    else:
+        force_fields = [mm.LoadCHARMMForcefield()]
+    force_fields.extend(extra_force_fields)
     # setup mm simulation
-    sim = _SetupMmSimulation(mhandle)
+    sim = _SetupMmSimulation(mhandle.model, force_fields)
     
     # settings to check for stereochemical problems
     clashing_distances = mol.alg.DefaultClashingDistances()
@@ -213,7 +334,7 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
     mhandle.model = mol.CreateEntityFromView(\
                             simulation_ent.Select("peptide=true and ele!=H"),
                             True)
-    # return model with oxygens
+    # return model with hydrogens
     return mol.CreateEntityFromView(simulation_ent.Select("peptide=true"), True)
 
 def CheckFinalModel(mhandle):

modelling/src/model.cc

@@ -479,16 +479,11 @@ void AddLigands(const ost::seq::AlignmentList& aln, XCSEditor& edi,
           continue;
         }
 
+        // deep copy of residue into ligand chain
         mol::ResidueHandle dst_res=edi.AppendResidue(ligand_chain,
-                                                     src_ligand.GetName());
+                                                     src_ligand.GetHandle(),
+                                                     true);
         dst_res.SetIsLigand(true);
-
-        for (AtomViewList::iterator j=atoms.begin(), 
-             e2=atoms.end(); j!=e2; ++j) {
-               edi.InsertAtom(dst_res, j->GetName(), 
-                              j->GetPos(), j->GetElement(), 
-                              j->GetOccupancy(), j->GetBFactor());
-        }
       }
     }
   }