diff --git a/CHANGELOG b/CHANGELOG
index e83fc8600d486f04ba3a8f6cd0de7f3468855d4d..24a36a9d5d65a22fa5a7a8c4227e63608e3a4936 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -5,6 +5,14 @@
Changelog
================================================================================
+
+Release 3.x.x
+--------------------------------------------------------------------------------
+
+* Check for and report non-planar rings in the modelling result. This check is
+ also incorporated in energy minimization, potentially increasing the number
+ of minimization steps.
+
Release 3.2.1
--------------------------------------------------------------------------------
diff --git a/modelling/doc/model_checking.rst b/modelling/doc/model_checking.rst
index ab5b5eb9eec534e44b685fb09571c30ce77ed2a7..400f3737f7389d157300c6bce938acf6abc68353 100644
--- a/modelling/doc/model_checking.rst
+++ b/modelling/doc/model_checking.rst
@@ -35,6 +35,15 @@ Detecting Ring Punches
.. autofunction:: FilterCandidatesWithSC
+
+Detecting Non-Planar Rings
+--------------------------------------------------------------------------------
+
+.. autofunction:: GetNonPlanarRings
+
+.. autofunction:: HasNonPlanarRings
+
+
Model Checking With MolProbity
--------------------------------------------------------------------------------
diff --git a/modelling/pymod/CMakeLists.txt b/modelling/pymod/CMakeLists.txt
index 58e2665c98f6d6ac1e81901a956547ddaa3a2836..c58f17c9c50a9433ad97832afc00e1677541dbcc 100644
--- a/modelling/pymod/CMakeLists.txt
+++ b/modelling/pymod/CMakeLists.txt
@@ -27,6 +27,7 @@ set(MODELLING_PYMOD
_mhandle_helper.py
_alignment_fiddling.py
_raw_model.py
+ _planar_rings.py
)
pymod(NAME modelling
diff --git a/modelling/pymod/__init__.py b/modelling/pymod/__init__.py
index 6cf9be125390e77c9307fdff5c346c56a2f540c0..a0eb62b66c12c597a411d19967d500d3f75918b5 100644
--- a/modelling/pymod/__init__.py
+++ b/modelling/pymod/__init__.py
@@ -27,3 +27,4 @@ from ._fragger_handle import *
from ._monte_carlo import *
from ._mhandle_helper import *
from ._raw_model import *
+from ._planar_rings import *
diff --git a/modelling/pymod/_pipeline.py b/modelling/pymod/_pipeline.py
index a4b5e8709578a0a1b79831442c2d8c93e22dc0e9..8b79f214120d3c6bf4d345cdd223f2e7f5f33bac 100644
--- a/modelling/pymod/_pipeline.py
+++ b/modelling/pymod/_pipeline.py
@@ -25,6 +25,7 @@ from ._modelling import *
from ._reconstruct_sidechains import *
from ._closegaps import *
from ._ring_punches import *
+from ._planar_rings import *
from ._mhandle_helper import *
# external
import ost
@@ -363,8 +364,9 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
ost.PopVerbosityLevel()
# checks above would remove bad atoms
cur_energy = sim.GetEnergy()
- if len(temp_ent_stereo_checked.Select("ele!=H").atoms) \
- == len(temp_ent.Select("ele!=H").atoms):
+ n_checked = len(temp_ent_stereo_checked.Select("ele!=H").atoms)
+ n_tot = len(temp_ent.Select("ele!=H").atoms)
+ if n_checked == n_tot and HasNonPlanarRings(temp_ent) == False:
ost.LogInfo("No more stereo-chemical problems "
"-> final energy: %g" % cur_energy)
break
@@ -396,6 +398,7 @@ def CheckFinalModel(mhandle):
- Remaining stereo-chemical problems after energy minimization. The affected
residues will have the boolean property "stereo_chemical_problem_backbone"
set to True, if the problem affects backbone atoms.
+ - Residues that contain non-planar rings.
:param mhandle: Modelling handle for which to perform checks.
:type mhandle: :class:`ModellingHandle`
@@ -477,6 +480,13 @@ def CheckFinalModel(mhandle):
ost.LogInfo(msg)
AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res])
+ # report non-planar rings
+ non_planar_rings = GetNonPlanarRings(mhandle.model)
+ for res in non_planar_rings:
+ msg = "Ring of " + str(res) + " is non-planar!"
+ ost.LogWarning(msg)
+ AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res])
+
def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list(),
model_termini=False):
'''Build a model starting with a raw model (see :func:`BuildRawModel`).
diff --git a/modelling/pymod/_planar_rings.py b/modelling/pymod/_planar_rings.py
new file mode 100644
index 0000000000000000000000000000000000000000..dcb4c28771200f0092be647dcfb7daf5d8a28d49
--- /dev/null
+++ b/modelling/pymod/_planar_rings.py
@@ -0,0 +1,113 @@
+# Copyright (c) 2013-2020, SIB - Swiss Institute of Bioinformatics and
+# Biozentrum - University of Basel
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+'''Helper functions to deal with non-planar rings'''
+import ost
+from ost import geom
+
+from promod3 import core
+
+def _GetMaxDist(res):
+ """Get max distance of any atom from ring plane
+
+ None is returned if res.one_letter_code is not in ['Y', 'F', 'W', 'H']
+ or when not all expected ring atoms are present
+
+ :param res: Residue for which you want to get the max dist to ring plane
+ :type res: :class:`~ost.mol.ResidueHandle` or :class:`~ost.mol.ResidueView`
+ :returns: Max distance of any ring atom to optimal ring plane or None if
+ residue contains no ring or not all expected ring atoms are
+ present.
+ """
+ prof = core.StaticRuntimeProfiler.StartScoped("planar_rings::GetMaxDist")
+ if res.one_letter_code in "YF":
+ atom_names = ["CG", "CD1", "CD2", "CE1", "CE2", "CZ"]
+ elif res.one_letter_code == 'W':
+ atom_names = ["CG", "CD1", "NE1", "CD2", "CE2", "CE3", "CZ2", "CZ3",
+ "CH2"]
+ elif res.one_letter_code == 'H':
+ atom_names = ["CG", "CD2", "ND1", "CE1", "NE2"]
+ else:
+ return None
+ positions = geom.Vec3List()
+ for atom_name in atom_names:
+ a = res.FindAtom(atom_name)
+ if a.IsValid():
+ positions.append(a.GetPos())
+ else:
+ return None
+ origin=positions.center
+ normal=positions.principal_axes.GetRow(0)
+ plane = geom.Plane(origin, normal)
+ return max([abs(geom.Distance(plane, pos)) for pos in positions])
+
+def GetNonPlanarRings(ent, max_dist_thresh = 0.1):
+ '''Get list of residues with rings that non-planar in the given structure.
+
+ Only residues with res.one_letter_code in ['Y', 'F', 'W', 'H'] that contain
+ all expected ring atoms are considered.
+
+ Planarity is defined by the maximum distance of any ring atom to the optimal
+ plane describing the ring. This plane is constructed by defining a point on
+ the plane, here the geometric center of the ring atoms and a normal. We
+ construct an orthogonal basis [e1, e2, e3]. e1 points in direction with
+ lowest variance of ring atoms and e3 in direction with highest variance.
+ e1 is thus the plane normal. Internally this is done using a singular value
+ decomposition algorithm.
+
+ To put the default threshold of 0.1 in perspective: if you calculate these
+ distances on the same non-redundant set of experimental structures as
+ |project| used to derive its rotamer libraries, 99.9 % of the residues are
+ within 0.065 (HIS), 0.075 (TRP), 0.057 (TYR), 0.060 (PHE).
+
+ :param ent: Structure for which to detect non-planar rings.
+ :type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
+ :param max_dist_thresh: A residue that contains a ring is considered
+ non-planar if the max distance of any ring-atom to
+ the optimal ring plane is above this threshold.
+ :type max_dist_thresh: :class:`float`
+
+ :return: :class:`list` of residues (:class:`~ost.mol.ResidueHandle`/
+ :class:`~ost.mol.ResidueView`) which have a non-planar ring.
+ '''
+ prof_name = 'ring_punches::GetNonPlanarRings'
+ prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
+ dists = [_GetMaxDist(res) for res in ent.residues]
+ mdt = max_dist_thresh
+ return [r for d, r in zip(dists, ent.residues) if d is not None and d > mdt]
+
+def HasNonPlanarRings(ent, max_dist_thresh = 0.1):
+ '''Check if any ring is non-planar in the given structure.
+
+ Calls :func:`GetNonPlanarRings` with given parameters and returns if any
+ residue is considered non-planar.
+
+ :param ent: Structure for which to detect non-planar rings
+ :type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
+ :param max_dist_thresh: A residue that contains a ring is considered
+ non-planar if the max distance of any ring-atom to
+ the optimal ring plane is above this threshold.
+ :type max_dist_thresh: :class:`float`
+
+
+ :return: True, if any ring is non-planar
+ :rtype: :class:`bool`
+ '''
+ prof_name = 'ring_punches::HasNonPlanarRings'
+ prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
+ return len(GetNonPlanarRings(ent, max_dist_thresh)) > 0
+
+# these methods will be exported into module
+__all__ = ('GetNonPlanarRings', 'HasNonPlanarRings')