diff --git a/loop/src/backbone_loop_score.cc b/loop/src/backbone_loop_score.cc
index 5400261e05e1aeb39d4ae7eb09cf468a96e33021..8c43233738d1b954ffadc3382b2c579346329ee4 100644
--- a/loop/src/backbone_loop_score.cc
+++ b/loop/src/backbone_loop_score.cc
@@ -1492,9 +1492,9 @@ void BackboneLoopScorer::SetEnvironment(const ost::mol::EntityHandle& env){
ost::mol::ChainHandleList chain_list = env.GetChainList();
ost::mol::ResidueHandleList res_list;
- if(chain_list.size() > num_chains_){
+ if(chain_list.size() < num_chains_){
std::stringstream ss;
- ss << "More chains than seqres sequences when setting the loop scorers environment! ";
+ ss << "Less chains than seqres sequences when setting the loop scorers environment! ";
ss << "Don't do anything...";
throw promod3::Error(ss.str());
}
@@ -1548,7 +1548,7 @@ void BackboneLoopScorer::SetEnvironment(const ost::mol::EntityHandle& env){
geom::Vec3 n_pos, ca_pos, c_pos, o_pos, cb_pos, old_pos;
uint num, index;
- for(uint i = 0; i < chain_list.size(); ++i){
+ for(uint i = 0; i < num_chains_; ++i){
res_list = chain_list[i].GetResidueList();
for(ost::mol::ResidueHandleList::iterator j = res_list.begin(),e = res_list.end();
diff --git a/modelling/doc/index.rst b/modelling/doc/index.rst
index 3ea3dff228a3cc3f8dbc53492d4865115dc4aca7..b5415e7707b9a0d06a7979db7d37f60d8358af6d 100644
--- a/modelling/doc/index.rst
+++ b/modelling/doc/index.rst
@@ -84,7 +84,13 @@ Modelling Pipeline
list of alignment handles for raw model with multiple chains.
:type aln: :class:`~ost.seq.AlignmentHandle` / :class:`~ost.seq.AlignmentList`
- :param include_ligands: True, if we wish to include ligands in the model.
+ :param include_ligands: True, if we wish to include ligands in the model. This
+ searches for ligands in all OST handles of the views
+ attached to the alignments. Ligands are identified
+ with the "ligand" property in the handle (set by OST
+ based on HET records) or by the chain name '_' (as set
+ in SMTL). All ligands are added to a new chain named
+ '_'.
:type include_ligands: :class:`bool`
:param chain_names: Chains are named by a single chanacter taken from this.
diff --git a/modelling/pymod/_pipeline.py b/modelling/pymod/_pipeline.py
index a09efd0241c10e109bf25b1dca37cfa49db6403d..9790e2ab8f1184ce4321cd2378c25281efb38a41 100644
--- a/modelling/pymod/_pipeline.py
+++ b/modelling/pymod/_pipeline.py
@@ -15,6 +15,56 @@ from ost import mol
from ost.mol import mm
import os
+###############################################################################
+# helper functions
+def _RemoveHydrogens(ent):
+ '''Hydrogen naming depends on the force field used.
+ To be on the safe side, we simply remove all of them.
+ Note that the pipeline ignores hydrogens when building the raw model and
+ when rebuilding sidechains.
+ '''
+ edi = ent.EditXCS()
+ ha = ent.Select('ele=H')
+ for a in ha.atoms:
+ edi.DeleteAtom(a.handle)
+
+def _SetupMmSimulation(mhandle):
+ '''Get mm simulation object for the model.'''
+ # 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)
+ _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!
+
+ # finally set up the simulation
+ try:
+ sim = mm.Simulation(model, 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)
+ ost.LogWarning("Switched to slower reference platform of OpenMM!")
+ return sim
+
+###############################################################################
+
def BuildSidechains(mhandle, merge_distance=4, scorer=None, fragment_db=None,
structure_db=None, torsion_sampler=None):
'''Build sidechains for model.
@@ -121,26 +171,8 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
# ignore LogInfo in stereochemical problems if output up to info done
ignore_stereo_log = (ost.GetVerbosityLevel() == 3)
# setup mm simulation
- settings = mm.Settings()
- settings.integrator = mm.LangevinIntegrator(310, 1, 0.002)
- settings.init_temperature = 0
- settings.forcefield = mm.LoadCHARMMForcefield()
- 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')
- try:
- sim = mm.Simulation(mhandle.model, 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(mhandle.model, settings)
- ost.LogWarning("Switched to slower reference platform of OpenMM!")
+ sim = _SetupMmSimulation(mhandle)
+
# settings to check for stereochemical problems
clashing_distances = mol.alg.DefaultClashingDistances()
bond_stereo_chemical_param = mol.alg.DefaultBondStereoChemicalParams()
@@ -148,7 +180,8 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
for i in range(max_iterations):
# update atoms
- ost.LogInfo("Perform energy minimization (iteration %d)" % (i+1))
+ ost.LogInfo("Perform energy minimization "
+ "(iteration %d, energy: %g)" % (i+1, sim.GetEnergy()))
sim.ApplySD(tolerance = 1.0, max_iterations = max_iter_sd)
sim.ApplyLBFGS(tolerance = 1.0, max_iterations = max_iter_lbfgs)
sim.UpdatePositions()
@@ -171,10 +204,9 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
# checks above would remove bad atoms
if len(temp_ent_stereo_checked.Select("ele!=H").atoms) \
== len(temp_ent.Select("ele!=H").atoms):
+ ost.LogInfo("No more stereo-chemical problems "
+ "-> final energy: %g" % (sim.GetEnergy()))
break
- else:
- ost.LogInfo("Stereo-chemical problems found, need more energy"
- + " minimization (iteration %d)" % (i+1))
# update model
simulation_ent = sim.GetEntity()
diff --git a/modelling/src/model.cc b/modelling/src/model.cc
index a90685cc3823b5a88a3751abe5e11932a86c783d..d53e4d2a0f5464291009e7464ed7d28257f748e5 100644
--- a/modelling/src/model.cc
+++ b/modelling/src/model.cc
@@ -431,17 +431,21 @@ ModellingHandle BuildRawModel(const ost::seq::AlignmentList& aln,
void AddLigands(const ost::seq::AlignmentList& aln, XCSEditor& edi,
EntityHandle& model)
{
-
+ // make sure we handle duplicates
+ std::set<mol::EntityHandle> unique_handles;
std::vector<mol::ResidueViewList> ligand_views;
for(ost::seq::AlignmentList::const_iterator it=aln.begin();
it!=aln.end();
++it) {
mol::EntityHandle handle = it->GetSequence(1).GetAttachedView().GetHandle();
- mol::EntityView ligand_view = handle.Select("ligand=true");
+ // did we already do that?
+ if (unique_handles.find(handle) != unique_handles.end()) continue;
+ mol::EntityView ligand_view = handle.Select("ligand=true or cname='_'");
mol::ResidueViewList ligands = ligand_view.GetResidueList();
if(!ligands.empty()) {
ligand_views.push_back(ligands);
}
+ unique_handles.insert(handle);
}
if (ligand_views.size()>0) {
diff --git a/sidechain/pymod/_reconstruct_sidechains.py b/sidechain/pymod/_reconstruct_sidechains.py
index 9b147775d790ac06447fd60a5dd6bc183feb86fb..c385aac2f632bb716d769501f334b0338494067f 100644
--- a/sidechain/pymod/_reconstruct_sidechains.py
+++ b/sidechain/pymod/_reconstruct_sidechains.py
@@ -23,7 +23,8 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
bbdep = True
residues_with_rotamer_group = list()
rotamer_groups = list()
- prot = ent.Select("peptide=true")
+ # ignore ligand chain and any non-peptides
+ prot = ent.Select("peptide=true and cname!='_'")
incomplete_sidechains = list()
#extract rotamer ids