SCHWED-882: tested external ff passing and added parameters to BuildFromRawModel

modelling/pymod/_closegaps.py

@@ -1225,11 +1225,12 @@ def CloseLargeDeletions(mhandle, scorer, structure_db, linker_length=8,
         # this is done to not loose all the sidechain information
         n_ter_bb_list = bb_list.Extract(0, frag_start_idx+1)
         t = initial_n_ter_bb_list.GetTransform(n_ter_bb_list)
-        ed = mhandle.model.EditXCS()
+        ed = mhandle.model.EditXCS(ost.mol.BUFFERED_EDIT)
         for r in actual_chain.residues[:frag_start_idx]:
             for a in r.atoms:
                 new_pos = t * ost.geom.Vec4(a.GetPos())
                 ed.SetAtomPos(a, ost.geom.Vec3(new_pos))
+        ed.UpdateICS()
 
         # replace fragment part
         fragment = bb_list.Extract(frag_start_idx, len(bb_list))

modelling/pymod/_pipeline.py

@@ -23,19 +23,21 @@ def _RemoveHydrogens(ent):
     Note that the pipeline ignores hydrogens when building the raw model and
     when rebuilding sidechains.
     '''
-    edi = ent.EditXCS()
+    edi = ent.EditXCS(mol.BUFFERED_EDIT)
     ha = ent.Select('ele=H')
     for a in ha.atoms:
-      edi.DeleteAtom(a.handle)
+        edi.DeleteAtom(a.handle)
+    edi.UpdateICS()
 
 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()
+            edi = ent.EditXCS(mol.BUFFERED_EDIT)
             h_constructor = mm.HeuristicHydrogenConstructor(bb)
             h_constructor.ApplyOnResidue(res, edi)
+            edi.UpdateICS()
 
 def _GetTopology(ent, settings, force_fields, add_heuristic_hydrogens=False):
     '''Return topology or None if no topology could be created.
@@ -43,8 +45,7 @@ def _GetTopology(ent, settings, force_fields, add_heuristic_hydrogens=False):
     Set add_heuristic_hydrogens to True for ligands.
     '''
     # try all force fields in order
-    last_exception = None
-    for ff in force_fields:
+    for i_ff, ff in enumerate(force_fields):
         settings.forcefield = ff
         try:
             # check if we need to add hydrogens heuristically
@@ -53,15 +54,14 @@ def _GetTopology(ent, settings, force_fields, add_heuristic_hydrogens=False):
             # 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
+            # report only for debugging
+            ost.LogVerbose("Could not create mm topology for ff %d. %s" \
+                           % (i_ff, 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):
@@ -80,13 +80,13 @@ def _AddLigands(ent, top, lig_ent, settings, force_fields):
             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:
+            view_res_str = str([str(r) for r in cur_view.residues])
             ost.LogError("Failed to add ligands " + view_res_str + \
                          " for energy minimization! Skipping...")
         else:
@@ -146,7 +146,7 @@ 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()
+    ent_ed = ent.EditXCS(mol.BUFFERED_EDIT)
     chain_names = [ch.name for ch in ent.chains]
     for chain_name in chain_names:
         # all separate ligand chains start with _
@@ -160,6 +160,7 @@ def _GetSimEntity(sim):
                 ent_ed.AppendResidue(lig_ch, res, deep=True)
             # remove old chain
             ent_ed.DeleteChain(cur_ch)
+    ent_ed.UpdateICS()
     return ent
 ###############################################################################
 
@@ -264,19 +265,12 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
     :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`.
+                         CHARMM one (see :meth:`BuildFromRawModel`).
     :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`.
+    :param extra_force_fields: Additional list of force fields to use (see
+                               :meth:`BuildFromRawModel`).
+    :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>`
@@ -398,7 +392,7 @@ def CheckFinalModel(mhandle):
             ost.LogInfo("Stereo-chemical problem in sidechain " + \
                         "of residue " + str(res))
 
-def BuildFromRawModel(mhandle):
+def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     '''Build a model starting with a raw model (see :func:`BuildRawModel`).
 
     This function implements a recommended pipeline to generate complete models
@@ -415,6 +409,21 @@ def BuildFromRawModel(mhandle):
                     alignment.
     :type mhandle:  :class:`ModellingHandle`
 
+    :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: Delivers the model as an |ost_s| entity.
     :rtype: :class:`Entity <ost.mol.EntityHandle>`
     '''
@@ -452,7 +461,8 @@ def BuildFromRawModel(mhandle):
                     structure_db, torsion_sampler)
 
     # minimize energy of final model using molecular mechanics
-    MinimizeModelEnergy(mhandle)
+    MinimizeModelEnergy(mhandle, use_amber_ff=use_amber_ff,
+                        extra_force_fields=extra_force_fields)
 
     # sanity checks
     CheckFinalModel(mhandle)

modelling/tests/CMakeLists.txt

@@ -30,12 +30,14 @@ set(MODELLING_TEST_DATA
   data/cbeta.pdb
   data/compounds.chemlib
   data/del.fasta
+  data/ff_NAD.dat
   data/gly.pdb
   data/hetero-punched.pdb
   data/ins.fasta
   data/sep.fasta
   data/sep.pdb
   data/seq.fasta
+  data/smtl2bl4.1.pdb
   data/ter.fasta
 )
 

modelling/tests/test_pipeline.py

@@ -324,6 +324,39 @@ class PipelineTests(unittest.TestCase):
         self.assertAlmostEqual(results['Ramachandran outliers'], 0.0, places=2)
         self.assertAlmostEqual(results['Rotamer outliers'], 0.0, places=2)
 
+    def testBuildModelWithLigands(self):
+        '''Check that we can run pipeline with a ligand.'''
+        # we get template with 5 ligands and model 1 chain
+        tpl = io.LoadPDB("data/smtl2bl4.1.pdb")
+        model_view = tpl.Select("cname=A")
+        seqres = ''.join([r.one_letter_code for r in model_view.residues])
+        aln = seq.CreateAlignment(seq.CreateSequence('trg', seqres),
+                                  seq.CreateSequence('tpl', seqres))
+        aln.AttachView(1, model_view)
+        mhandle = modelling.BuildRawModel(aln, include_ligands=True)
+        # 2 ligands (NAD, FE2) touching chain A
+        ch = mhandle.model.FindChain('_')
+        self.assertTrue(ch.IsValid())
+        self.assertEqual(ch.residue_count, 2)
+        # model it without extra force fields
+        # -> should see ligands for sidechains and remove them in energy min.
+        # -> we excpect two error logs for skipped ligands
+        log = _FetchLog()
+        ost.PushLogSink(log)
+        log.messages['ERROR'] = list()
+        final_model = modelling.BuildFromRawModel(mhandle)
+        self.assertFalse(final_model.FindChain('_').IsValid())
+        self.assertEqual(len(log.messages['ERROR']), 2)
+        # model it with extra force field
+        # -> only NAD ligand in ff file, should be kept
+        ff = mol.mm.Forcefield.Load("data/ff_NAD.dat")
+        mhandle = modelling.BuildRawModel(aln, include_ligands=True)
+        final_model = modelling.BuildFromRawModel(mhandle, True, [ff])
+        ch = final_model.FindChain('_')
+        self.assertTrue(ch.IsValid())
+        self.assertEqual(ch.residue_count, 1)
+        self.assertEqual(ch.residues[0].name, "NAD")
+
 if __name__ == "__main__":
     from ost import testutils
     testutils.RunTests()