diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 8553f53977e515837cb5e3ddf38e5fe01e643d1f..a74f28d4d27440806392edbeba1ee36255b67351 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,4 +1,4 @@
-Changes in Release 1.
+Changes in Release <RELEASE NUMBER>
--------------------------------------------------------------------------------
* nonstandard C++ module was moved from ost.conop to ost.mol.alg. This implies
@@ -6,7 +6,19 @@ Changes in Release 1.
CopyNonConserved that were previousely imported from ost.conop are now
to be imported from ost.mol.alg.
* Added Molck API to the ost.mol.alg module.
+ * Removed habit of changing secondary structure of entities when loading
+ from mmCIF PDB files. Before, OST would turn secondary structure 'EEH'
+ into 'ECH' to make it look nicer in DNG. Now, 'EEH' stays 'EEH'.
+Changes in Release 1.7.1
+--------------------------------------------------------------------------------
+
+ * Fixed an issue that could cause the star format parser (mmCIF, chemical
+ components dictionary) to enter an infinite loop
+ * Chemical components dictionary was extended by new chemical classes
+ introduced by PDB
+ * Fixed unit tests
+ * Improved documentation
Changes in Release 1.7
--------------------------------------------------------------------------------
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2f582069ba52a1b79e281c8866d6f591c5f1942e..2cf4c283a4cbfd59833a4d76837da48948dbe9f1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,7 +7,7 @@ project(OpenStructure CXX C)
set (CMAKE_EXPORT_COMPILE_COMMANDS 1)
set (OST_VERSION_MAJOR 1)
set (OST_VERSION_MINOR 7)
-set (OST_VERSION_PATCH 0)
+set (OST_VERSION_PATCH 1)
set (OST_VERSION_STRING ${OST_VERSION_MAJOR}.${OST_VERSION_MINOR}.${OST_VERSION_PATCH} )
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake_support)
include(OST)
diff --git a/modules/io/doc/mmcif.rst b/modules/io/doc/mmcif.rst
index 730ca1b8e97ed554287158184c9d9eb20283da37..23dc32f461a08cf1633362e60bb3b9c54ecb09af 100644
--- a/modules/io/doc/mmcif.rst
+++ b/modules/io/doc/mmcif.rst
@@ -638,7 +638,7 @@ of the annotation available.
.. method:: GetChainIntervalList()
- See :attr:`chainintervalls`
+ See :attr:`chainintervals`
.. method:: GetOperations()
diff --git a/modules/io/src/mol/mmcif_reader.cc b/modules/io/src/mol/mmcif_reader.cc
index 27a7c9bc8406308409525862a7b1d906417d2c36..905e28442bab8bd68300426bdeb99df8c6d4cb0e 100644
--- a/modules/io/src/mol/mmcif_reader.cc
+++ b/modules/io/src/mol/mmcif_reader.cc
@@ -1543,17 +1543,7 @@ void MMCifReader::AssignSecStructure(mol::EntityHandle ent)
continue;
}
mol::SecStructure alpha(mol::SecStructure::ALPHA_HELIX);
- // some PDB files contain helix/strand entries that are adjacent to each
- // other. To avoid visual artifacts, we effectively shorten the first of
- // the two secondary structure segments to insert one residue of coil
- // conformation.
- mol::ResNum start = i->start, end = i->end;
- if (helix_list_.end() != i+1 && // unit test
- (*(i+1)).start.GetNum() <= end.GetNum()+1 &&
- (*(i+1)).end.GetNum() > end.GetNum()) {
- end = mol::ResNum((*(i+1)).start.GetNum()-2);
- }
- chain.AssignSecondaryStructure(alpha, start, end);
+ chain.AssignSecondaryStructure(alpha, i->start, i->end);
}
for (MMCifHSVector::const_iterator i=strand_list_.begin(),
@@ -1565,14 +1555,7 @@ void MMCifReader::AssignSecStructure(mol::EntityHandle ent)
continue;
}
mol::SecStructure extended(mol::SecStructure::EXTENDED);
- mol::ResNum start = i->start, end = i->end;
- // see comment for helix assignment
- if (strand_list_.end() != i+1 && // unit test
- (*(i+1)).start.GetNum() <= end.GetNum()+1 &&
- (*(i+1)).end.GetNum() > end.GetNum()) {
- end=mol::ResNum((*(i+1)).start.GetNum()-2);
- }
- chain.AssignSecondaryStructure(extended, start, end);
+ chain.AssignSecondaryStructure(extended, i->start, i->end);
}
}
diff --git a/modules/io/src/mol/pdb_reader.cc b/modules/io/src/mol/pdb_reader.cc
index a203234ecc275f82bceb87b75781ad1e626fb100..ffa04e6ecc0144a342bd477c84ab35b570f78e5a 100644
--- a/modules/io/src/mol/pdb_reader.cc
+++ b/modules/io/src/mol/pdb_reader.cc
@@ -520,16 +520,7 @@ void PDBReader::AssignSecStructure(mol::EntityHandle ent)
continue;
}
mol::SecStructure alpha(mol::SecStructure::ALPHA_HELIX);
- // some PDB files contain helix/strand entries that are adjacent to each
- // other. To avoid visual artifacts, we effectively shorten the first of the
- // two secondary structure segments to insert one residue of coil
- // conformation.
- mol::ResNum start=i->start, end=i->end;
- if (helix_list_.end()!=i+1 && (*(i+1)).start.GetNum()<=end.GetNum()+1 &&
- (*(i+1)).end.GetNum()>end.GetNum()) {
- end=mol::ResNum((*(i+1)).start.GetNum()-2);
- }
- chain.AssignSecondaryStructure(alpha, start, end);
+ chain.AssignSecondaryStructure(alpha, i->start, i->end);
}
for (HSList::const_iterator i=strand_list_.begin(), e=strand_list_.end();
@@ -540,13 +531,7 @@ void PDBReader::AssignSecStructure(mol::EntityHandle ent)
continue;
}
mol::SecStructure extended(mol::SecStructure::EXTENDED);
- mol::ResNum start=i->start, end=i->end;
- // see comment for helix assignment
- if (strand_list_.end()!=i+1 && (*(i+1)).start.GetNum()<=end.GetNum()+1 &&
- (*(i+1)).end.GetNum()>end.GetNum()) {
- end=mol::ResNum((*(i+1)).start.GetNum()-2);
- }
- chain.AssignSecondaryStructure(extended, start, end);
+ chain.AssignSecondaryStructure(extended, i->start, i->end);
}
}
diff --git a/modules/io/src/mol/star_parser.cc b/modules/io/src/mol/star_parser.cc
index 80043259b4df5fdf34ee81de6bf471696d03bc5a..214eb626c9f743a1cbbb218bede0ad7f795b2f73 100644
--- a/modules/io/src/mol/star_parser.cc
+++ b/modules/io/src/mol/star_parser.cc
@@ -550,11 +550,15 @@ void StarParser::Parse()
case 'd':
if (tline.length()>=5 && StringRef("data_", 5)==tline.substr(0, 5)) {
this->ParseData();
+ } else {
+ throw IOException("Missing 'data_' control structure");
}
break;
case 'g':
if (tline.length()>=7 && StringRef("global_", 7)==tline.substr(0, 7)) {
this->ParseGlobal();
+ } else {
+ throw IOException("Missing 'global_' control structure");
}
break;
case '#':
diff --git a/modules/io/tests/test_star_parser.cc b/modules/io/tests/test_star_parser.cc
index 6ff217b748c123770f789a2168d540be3f1ef94d..c4e40fdbda137d74293282a592884fb5aa7e8305 100644
--- a/modules/io/tests/test_star_parser.cc
+++ b/modules/io/tests/test_star_parser.cc
@@ -389,6 +389,20 @@ BOOST_AUTO_TEST_CASE(star_items_as_row)
BOOST_TEST_MESSAGE(" done.");
}
+BOOST_AUTO_TEST_CASE(star_broken_data)
+{
+ std::ifstream s("testfiles/broken_data.cif");
+ LoopTestParser star_p(s);
+ BOOST_CHECK_THROW(star_p.Parse(), IOException);
+}
+
+BOOST_AUTO_TEST_CASE(star_broken_global)
+{
+ std::ifstream s("testfiles/broken_global.cif");
+ LoopTestParser star_p(s);
+ BOOST_CHECK_THROW(star_p.Parse(), IOException);
+}
+
BOOST_AUTO_TEST_CASE(star_missing_data)
{
std::ifstream s("testfiles/missing_data.cif");
diff --git a/modules/io/tests/testfiles/broken_data.cif b/modules/io/tests/testfiles/broken_data.cif
new file mode 100644
index 0000000000000000000000000000000000000000..2ae242255c2bd8c6e952ffee6d1815ee01f6335a
--- /dev/null
+++ b/modules/io/tests/testfiles/broken_data.cif
@@ -0,0 +1 @@
+d lines starting with a 'd' kicked the parser into an infinite loop.
diff --git a/modules/io/tests/testfiles/broken_global.cif b/modules/io/tests/testfiles/broken_global.cif
new file mode 100644
index 0000000000000000000000000000000000000000..0f340e0befe0077229b53b121979181ec3c890fc
--- /dev/null
+++ b/modules/io/tests/testfiles/broken_global.cif
@@ -0,0 +1 @@
+g lines starting with a 'g' kicked the parser into an infinite loop.
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 1816100c024249600998e0aae690dc8d9e25e072..d9ba6a6a60588da14ff02de9793aa46d899bd69a 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -957,6 +957,102 @@ Algorithms on Structures
:class:`~ost.mol.EntityHandle`
+
+.. class:: FindMemParam
+
+ Result object for the membrane detection algorithm described below
+
+ .. attribute:: axis
+
+ initial search axis from which optimal membrane slab could be found
+
+ .. attribute:: tilt_axis
+
+ Axis around which we tilt the membrane starting from the initial axis
+
+ .. attribute:: tilt
+
+ Angle to tilt around tilt axis
+
+ .. attribute:: angle
+
+ After the tilt operation we perform a rotation around the initial axis
+ with this angle to get the final membrane axis
+
+ .. attribute:: membrane_axis
+
+ The result of applying the tilt and rotation procedure described above.
+ The membrane_axis is orthogonal to the membrane plane and has unit length.
+
+ .. attribute:: pos
+
+ Real number that describes the membrane center point. To get the actual
+ position you can do: pos * membrane_axis
+
+ .. attribute:: width
+
+ Total width of the membrane in A
+
+ .. attribute:: energy
+
+ Pseudo energy of the implicit solvation model
+
+ .. attribute:: membrane_representation
+
+ Dummy atoms that represent the membrane. This entity is only valid if
+ the according flag has been set to True when calling FindMembrane.
+
+
+.. method:: FindMembrane(ent, assign_membrane_representation=True, fast=False)
+
+ Estimates the optimal membrane position of a protein by using an implicit
+ solvation model. The original algorithm and the used energy function are
+ described in: Lomize AL, Pogozheva ID, Lomize MA, Mosberg HI (2006)
+ Positioning of proteins in membranes: A computational approach.
+
+ There are some modifications in this implementation and the procedure is
+ as follows:
+
+ * Initial axis are constructed that build the starting point for initial
+ parameter grid searches.
+
+ * For every axis, the protein is rotated so that the axis builds the z-axis
+
+ * In order to exclude internal hydrophilic pores, only the outermost atoms
+ with respect the the z-axis enter an initial grid search
+ * The width and position of the membrane is optimized for different
+ combinations of tilt and rotation angles (further described in
+ :class:`FindMemParam`). The top 20 parametrizations
+ (only top parametrization if *fast* is True) are stored for further
+ processing.
+
+ * The 20 best membrane parametrizations from the initial grid search
+ (only the best if *fast* is set to True) enter a final
+ minimization step using a Levenberg-Marquardt minimizer.
+
+
+ :param ent: Entity of a transmembrane protein, you'll get weird
+ results if this is not the case. The energy term
+ of the result is typically a good indicator whether
+ *ent* is an actual transmembrane protein.
+ :type ent: :class:`ost.mol.EntityHandle` / :class:`ost.mol.EntityView`
+
+ :param assign_membrane_representation: Whether to construct a membrane
+ representation using dummy atoms
+
+ :type assign_membrane_representation: :class:`bool`
+
+ :param fast: If set to false, the 20 best results of the initial grid
+ search undergo a Levenberg-Marquardt minimization and
+ the parametrization with optimal minimized energy is
+ returned.
+ If set to yes, only the best result of the initial grid
+ search is selected and returned after
+ Levenberg-Marquardt minimization.
+
+ :returns: The results object
+ :rtype: :class:`ost.mol.alg.FindMemParam`
+
.. _traj-analysis:
Trajectory Analysis
diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index 5519c6adf2741e213eba44cb82fbaf559f4c19e8..f4d1fdb71b68519b303d5df13dea34a2d529c230 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -9,6 +9,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
export_sec_structure.cc
export_non_standard.cc
export_molck.cc
+ export_membrane.cc
)
set(OST_MOL_ALG_PYMOD_MODULES
diff --git a/modules/mol/alg/pymod/export_membrane.cc b/modules/mol/alg/pymod/export_membrane.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fe8ebe52de9e43c2a98499718f1951015794453f
--- /dev/null
+++ b/modules/mol/alg/pymod/export_membrane.cc
@@ -0,0 +1,59 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+
+
+#include <boost/python.hpp>
+
+#include <ost/mol/alg/find_membrane.hh>
+
+using namespace boost::python;
+
+namespace{
+ ost::mol::alg::FindMemParam FindMembraneView(ost::mol::EntityView& v,
+ bool assign_membrane_representation,
+ bool fast) {
+ return ost::mol::alg::FindMembrane(v, assign_membrane_representation, fast);
+ }
+
+ ost::mol::alg::FindMemParam FindMembraneHandle(ost::mol::EntityHandle& h,
+ bool assign_membrane_representation,
+ bool fast) {
+ return ost::mol::alg::FindMembrane(h, assign_membrane_representation, fast);
+ }
+}
+
+void export_find_membrane() {
+
+ class_<ost::mol::alg::FindMemParam>("FindMemParam", no_init)
+ .def_readonly("tilt", &ost::mol::alg::FindMemParam::tilt)
+ .def_readonly("angle", &ost::mol::alg::FindMemParam::angle)
+ .def_readonly("width", &ost::mol::alg::FindMemParam::width)
+ .def_readonly("pos", &ost::mol::alg::FindMemParam::pos)
+ .def_readonly("energy", &ost::mol::alg::FindMemParam::energy)
+ .def_readonly("axis", &ost::mol::alg::FindMemParam::axis)
+ .def_readonly("tilt_axis", &ost::mol::alg::FindMemParam::tilt_axis)
+ .def_readonly("membrane_axis", &ost::mol::alg::FindMemParam::GetMembraneAxis)
+ .def_readonly("membrane_representation", &ost::mol::alg::FindMemParam::membrane_representation)
+ ;
+
+ def("FindMembrane", FindMembraneView, (arg("ent"), arg("assign_membrane_representation")=true,
+ arg("fast")=false));
+ def("FindMembrane", FindMembraneHandle, (arg("ent"), arg("assign_membrane_representation")=true,
+ arg("fast")=false));
+}
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index f32e3256f656586953ffa654c6617cb0377781cb..1642c2cafe124f4eb189d2ddb8d4d308458e0309 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -44,6 +44,7 @@ void export_Molck();
void export_contact_overlap();
void export_accessibility();
void export_sec_struct();
+void export_find_membrane();
#if OST_IMG_ENABLED
void export_entity_to_density();
#endif
@@ -117,6 +118,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
export_contact_overlap();
export_accessibility();
export_sec_struct();
+ export_find_membrane();
#if OST_IMG_ENABLED
export_entity_to_density();
#endif
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index ef2ae22ae10ccc9172285b7d79e96b298c7d6fec..42c571d79444add71b47a8da9916a106a2e1caf1 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -22,6 +22,7 @@ set(OST_MOL_ALG_HEADERS
sec_struct.hh
nonstandard.hh
molck.hh
+ find_membrane.hh
)
set(OST_MOL_ALG_SOURCES
@@ -47,6 +48,7 @@ set(OST_MOL_ALG_SOURCES
sec_struct.cc
nonstandard.cc
molck.cc
+ find_membrane.cc
)
set(MOL_ALG_DEPS ost_mol ost_seq)
diff --git a/modules/mol/alg/src/accessibility.cc b/modules/mol/alg/src/accessibility.cc
index a1becdd7f83c3de0b939e22788d6bf4bb3f22926..276b53ba2256802e9c9afc52c2f4b49e308e796f 100644
--- a/modules/mol/alg/src/accessibility.cc
+++ b/modules/mol/alg/src/accessibility.cc
@@ -203,6 +203,11 @@ Real GetAtomAccessibilityNACCESS(Real x_pos, Real y_pos, Real z_pos,
Real a = x_pos - x[i];
Real b = y_pos - y[i];
Real c = a*a + b*b;
+
+ if(c == Real(0.0)) {
+ return 0.0;
+ }
+
dx[i] = a;
dy[i] = b;
dsqr[i] = c;
diff --git a/modules/mol/alg/src/find_membrane.cc b/modules/mol/alg/src/find_membrane.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f48fb1149815f65e162e2f39a2c961dff494a0a1
--- /dev/null
+++ b/modules/mol/alg/src/find_membrane.cc
@@ -0,0 +1,1018 @@
+#include <ost/mol/alg/find_membrane.hh>
+#include <ost/mol/alg/accessibility.hh>
+#include <ost/geom/vecmat3_op.hh>
+#include <ost/message.hh>
+
+#include <limits>
+#include <exception>
+#include <list>
+#include <cmath>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+
+namespace{
+
+// Copyright notice of the Levenberg Marquardt minimizer we use...
+
+// Copyright (c) 2007, 2008, 2009 libmv authors.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to
+// deal in the Software without restriction, including without limitation the
+// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+// sell copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+// IN THE SOFTWARE.
+//
+// A simple implementation of levenberg marquardt.
+//
+// [1] K. Madsen, H. Nielsen, O. Tingleoff. Methods for Non-linear Least
+// Squares Problems.
+// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
+//
+// TODO(keir): Cite the Lourakis' dogleg paper.
+
+template<typename Function,
+ typename Jacobian,
+ typename Solver = Eigen::JacobiSVD<
+ Eigen::Matrix<typename Function::FMatrixType::RealScalar,
+ Function::XMatrixType::RowsAtCompileTime,
+ Function::XMatrixType::RowsAtCompileTime> > >
+class LevenbergMarquardt {
+ public:
+ typedef typename Function::XMatrixType::RealScalar Scalar;
+ typedef typename Function::FMatrixType FVec;
+ typedef typename Function::XMatrixType Parameters;
+ typedef Eigen::Matrix<typename Function::FMatrixType::RealScalar,
+ Function::FMatrixType::RowsAtCompileTime,
+ Function::XMatrixType::RowsAtCompileTime> JMatrixType;
+ typedef Eigen::Matrix<typename JMatrixType::RealScalar,
+ JMatrixType::ColsAtCompileTime,
+ JMatrixType::ColsAtCompileTime> AMatrixType;
+
+ // TODO(keir): Some of these knobs can be derived from each other and
+ // removed, instead of requiring the user to set them.
+ enum Status {
+ RUNNING,
+ GRADIENT_TOO_SMALL, // eps > max(J'*f(x))
+ RELATIVE_STEP_SIZE_TOO_SMALL, // eps > ||dx|| / ||x||
+ ERROR_TOO_SMALL, // eps > ||f(x)||
+ HIT_MAX_ITERATIONS,
+ };
+
+ LevenbergMarquardt(const Function &f)
+ : f_(f), df_(f) {}
+
+ struct SolverParameters {
+ SolverParameters()
+ : gradient_threshold(1e-20),
+ relative_step_threshold(1e-20),
+ error_threshold(1e-16),
+ initial_scale_factor(1e-3),
+ max_iterations(200) {}
+ Scalar gradient_threshold; // eps > max(J'*f(x))
+ Scalar relative_step_threshold; // eps > ||dx|| / ||x||
+ Scalar error_threshold; // eps > ||f(x)||
+ Scalar initial_scale_factor; // Initial u for solving normal equations.
+ int max_iterations; // Maximum number of solver iterations.
+ };
+
+ struct Results {
+ Scalar error_magnitude; // ||f(x)||
+ Scalar gradient_magnitude; // ||J'f(x)||
+ int iterations;
+ Status status;
+ };
+
+ Status Update(const Parameters &x, const SolverParameters ¶ms,
+ JMatrixType *J, AMatrixType *A, FVec *error, Parameters *g) {
+ *J = df_(x);
+ *A = (*J).transpose() * (*J);
+ *error = -f_(x);
+ *g = (*J).transpose() * *error;
+ if (g->array().abs().maxCoeff() < params.gradient_threshold) {
+ return GRADIENT_TOO_SMALL;
+ } else if (error->norm() < params.error_threshold) {
+ return ERROR_TOO_SMALL;
+ }
+ return RUNNING;
+ }
+
+ Results minimize(Parameters *x_and_min) {
+ SolverParameters params;
+ return minimize(params, x_and_min);
+ }
+
+ Results minimize(const SolverParameters ¶ms, Parameters *x_and_min) {
+ Parameters &x = *x_and_min;
+ JMatrixType J;
+ AMatrixType A;
+ FVec error;
+ Parameters g;
+
+ Results results;
+ results.status = Update(x, params, &J, &A, &error, &g);
+
+ Scalar u = Scalar(params.initial_scale_factor*A.diagonal().maxCoeff());
+ Scalar v = 2;
+
+ Parameters dx, x_new;
+ int i;
+ for (i = 0; results.status == RUNNING && i < params.max_iterations; ++i) {
+// LOG(INFO) << "iteration: " << i;
+// LOG(INFO) << "||f(x)||: " << f_(x).norm();
+// LOG(INFO) << "max(g): " << g.array().abs().maxCoeff();
+// LOG(INFO) << "u: " << u;
+// LOG(INFO) << "v: " << v;
+ AMatrixType A_augmented = A + u*AMatrixType::Identity(J.cols(), J.cols());
+ Solver solver(A_augmented, Eigen::ComputeThinU | Eigen::ComputeThinV);
+ dx = solver.solve(g);
+ if (dx.norm() <= params.relative_step_threshold * x.norm()) {
+ results.status = RELATIVE_STEP_SIZE_TOO_SMALL;
+ break;
+ }
+
+ x_new = x + dx;
+ // Rho is the ratio of the actual reduction in error to the reduction
+ // in error that would be obtained if the problem was linear.
+ // See [1] for details.
+ Scalar rho((error.squaredNorm() - f_(x_new).squaredNorm())
+ / dx.dot(u*dx + g));
+ if (rho > 0) {
+ // Accept the Gauss-Newton step because the linear model fits well.
+ x = x_new;
+ results.status = Update(x, params, &J, &A, &error, &g);
+ Scalar tmp = Scalar(2*rho-1);
+ u = u*std::max(Scalar(1/3.), 1 - (tmp*tmp*tmp));
+ v = 2;
+ continue;
+ }
+
+ // Reject the update because either the normal equations failed to solve
+ // or the local linear model was not good (rho < 0). Instead, increase
+ // to move closer to gradient descent.
+ u *= v;
+ v *= 2;
+ }
+ if (results.status == RUNNING) {
+ results.status = HIT_MAX_ITERATIONS;
+ }
+ results.error_magnitude = error.norm();
+ results.gradient_magnitude = g.norm();
+ results.iterations = i;
+ return results;
+ }
+
+ private:
+ const Function &f_;
+ Jacobian df_;
+};
+
+geom::Mat3 RotationAroundAxis(geom::Vec3 axis, Real angle) {
+
+ Real aa, ab, ac, ba, bb, bc, ca, cb, cc, one_m_cos, cos_ang, sin_ang;
+
+ cos_ang = std::cos(angle);
+ sin_ang = std::sin(angle);
+ one_m_cos = 1-cos_ang;
+
+ aa = cos_ang+axis[0]*axis[0]*one_m_cos;
+ ab = axis[0]*axis[1]*one_m_cos-axis[2]*sin_ang;
+ ac = axis[0]*axis[2]*one_m_cos+axis[1]*sin_ang;
+
+ ba = axis[1]*axis[0]*one_m_cos+axis[2]*sin_ang;
+ bb = cos_ang+axis[1]*axis[1]*one_m_cos;
+ bc = axis[1]*axis[2]*one_m_cos-axis[0]*sin_ang;
+
+ ca = axis[2]*axis[0]*one_m_cos-axis[1]*sin_ang;
+ cb = axis[2]*axis[1]*one_m_cos+axis[0]*sin_ang;
+ cc = cos_ang+axis[2]*axis[2]*one_m_cos;
+
+ geom::Mat3 result(aa, ab, ac, ba, bb, bc, ca, cb, cc);
+ return result;
+}
+
+geom::Vec3 RotateAroundAxis(geom::Vec3 point, geom::Vec3 axis, Real angle) {
+
+ geom::Mat3 rot = RotationAroundAxis(axis, angle);
+ geom::Vec3 result = rot*point;
+ return result;
+}
+
+
+// Levenberg Marquardt specific objects for the membrane finding algorithm
+
+struct EnergyF {
+ EnergyF(const std::vector<geom::Vec3>& p, const std::vector<Real>& t_e,
+ Real l, Real o, const geom::Vec3& ax, const geom::Vec3& tilt_ax):
+ positions(p), transfer_energies(t_e), axis(ax),
+ tilt_axis(tilt_ax), lambda(l), offset(o) { }
+
+ typedef Eigen::Matrix<Real, 4, 1> XMatrixType;
+ typedef Eigen::Matrix<Real, 1, 1> FMatrixType;
+
+ Eigen::Matrix<Real,1,1> operator()(const Eigen::Matrix<Real, 4, 1>& x) const {
+
+ FMatrixType result;
+ result(0,0) = 0.0;
+ geom::Vec3 tilted_axis = axis;
+
+ tilted_axis = RotateAroundAxis(tilted_axis, tilt_axis, x(0, 0));
+ tilted_axis = RotateAroundAxis(tilted_axis, axis, x(1, 0));
+
+ Real pos_on_axis;
+ Real distance_to_center;
+ Real half_width = Real(0.5) * x(2, 0);
+ Real exponent;
+ Real one_over_lambda = Real(1.0) / lambda;
+
+ int n = transfer_energies.size();
+
+ for(int i = 0; i < n; ++i) {
+ pos_on_axis = geom::Dot(tilted_axis, positions[i]);
+ distance_to_center = std::abs(x(3, 0) - pos_on_axis);
+ exponent = (distance_to_center - half_width) * one_over_lambda;
+ result(0, 0) += (1.0 / (1.0 + std::exp(exponent))) * transfer_energies[i];
+ }
+
+ result(0, 0) += offset;
+
+ return result;
+ }
+
+ std::vector<geom::Vec3> positions;
+ std::vector<Real> transfer_energies;
+ geom::Vec3 axis;
+ geom::Vec3 tilt_axis;
+ Real lambda;
+
+ // define an offset parameter...
+ // The levenberg-marquardt algorithm is designed to minimize an error
+ // function, aims to converge towards zero. The offset parameter gets added
+ // at the end of the energy calculation to get a positive result =>
+ // forces algorithm to minimize
+ Real offset;
+};
+
+
+struct EnergyDF {
+
+ EnergyDF(const EnergyF& f): function(f),d_tilt(0.02),d_angle(0.02),
+ d_width(0.4),d_pos(0.4) { }
+
+ Eigen::Matrix<Real,1,4> operator()(const Eigen::Matrix<Real, 4, 1>& x) const {
+
+ Eigen::Matrix<Real,1,4> result;
+ Eigen::Matrix<Real, 4, 1> parameter1 = x;
+ Eigen::Matrix<Real, 4, 1> parameter2 = x;
+
+ parameter1(0,0)+=d_tilt;
+ parameter2(0,0)-=d_tilt;
+ result(0,0) = (function(parameter1)(0, 0) -
+ function(parameter2)(0, 0)) / (2*d_tilt);
+
+ parameter1=x;
+ parameter2=x;
+ parameter1(1,0)+=d_angle;
+ parameter2(1,0)-=d_angle;
+ result(0,1) = (function(parameter1)(0,0) -
+ function(parameter2)(0,0)) / (2*d_angle);
+
+ parameter1=x;
+ parameter2=x;
+ parameter1(2,0)+=d_width;
+ parameter2(2,0)-=d_width;
+ result(0,2) = (function(parameter1)(0,0) -
+ function(parameter2)(0,0)) / (2*d_width);
+
+ parameter1=x;
+ parameter2=x;
+ parameter1(3,0)+=d_pos;
+ parameter2(3,0)-=d_pos;
+ result(0,3) = (function(parameter1)(0,0) -
+ function(parameter2)(0,0)) / (2*d_pos);
+ return result;
+ }
+
+ EnergyF function;
+ Real d_tilt;
+ Real d_angle;
+ Real d_width;
+ Real d_pos;
+};
+
+
+void GetRanges(const std::vector<geom::Vec3>& atom_positions,
+ Real& min_x, Real& max_x, Real& min_y, Real& max_y,
+ Real& min_z, Real& max_z) {
+
+ min_x = std::numeric_limits<Real>::max();
+ max_x = -min_x;
+
+ min_y = std::numeric_limits<Real>::max();
+ max_y = -min_y;
+
+ min_z = std::numeric_limits<Real>::max();
+ max_z = -min_z;
+
+ for(uint i = 0; i < atom_positions.size(); ++i) {
+ const geom::Vec3& pos = atom_positions[i];
+ min_x = std::min(min_x, pos[0]);
+ max_x = std::max(max_x, pos[0]);
+ min_y = std::min(min_y, pos[1]);
+ max_y = std::max(max_y, pos[1]);
+ min_z = std::min(min_z, pos[2]);
+ max_z = std::max(max_z, pos[2]);
+ }
+}
+
+
+void FloodLevel(char* data, int x_start, int y_start,
+ int x_extent, int y_extent,
+ int orig_value, int dest_value) {
+
+ //http://lodev.org/cgtutor/floodfill.html
+ if(orig_value != data[x_start*y_extent + y_start]) {
+ return;
+ }
+
+ std::vector<std::pair<int,int> > queue;
+ queue.push_back(std::make_pair(x_start, y_start));
+
+ int y1,y,x;
+ bool spanLeft, spanRight;
+ std::pair<int,int> actual_position;
+
+ while(!queue.empty()){
+
+ actual_position = queue.back();
+ queue.pop_back();
+ x = actual_position.first;
+ y = actual_position.second;
+
+ y1 = y;
+
+ while(y1 >= 0 && data[x*y_extent + y1] == orig_value) {
+ --y1;
+ }
+
+ y1++;
+ spanLeft = spanRight = 0;
+
+ while(y1 < y_extent &&
+ data[x*y_extent + y1] == orig_value ) {
+
+ data[x*y_extent + y1] = dest_value;
+
+ if(!spanLeft && x > 0 &&
+ data[(x-1)*y_extent + y1] == orig_value) {
+ queue.push_back(std::make_pair(x-1,y1));
+ spanLeft = 1;
+ }
+ else if(spanLeft && x > 0 &&
+ data[(x-1)*y_extent + y1] != orig_value) {
+ spanLeft = 0;
+ }
+
+ if(!spanRight && x < x_extent - 1 &&
+ data[(x+1)*y_extent + y1] == orig_value) {
+ queue.push_back(std::make_pair(x+1,y1));
+ spanRight = 1;
+ }
+ else if(spanRight && x < x_extent - 1 &&
+ data[(x+1)*y_extent + y1] != orig_value) {
+ spanRight = 0;
+ }
+ ++y1;
+ }
+ }
+}
+
+
+void GetExposedAtoms(const std::vector<geom::Vec3>& atom_positions,
+ const std::vector<Real>& transfer_energies,
+ std::vector<geom::Vec3>& exposed_atom_positions,
+ std::vector<Real>& exposed_transfer_energies) {
+
+ // sum of approx. vdw radius of the present heavy atoms (1.8)
+ // plus 1.4 for water.
+ Real radius = 3.2;
+ Real one_over_radius = Real(1.0) / radius;
+
+ // lets setup a grid in which we place the atoms
+ Real min_x, max_x, min_y, max_y, min_z, max_z;
+ GetRanges(atom_positions, min_x, max_x, min_y, max_y, min_z, max_z);
+
+ // we guarantee that the thing is properly solvated in the x-y plane and add
+ // some space around this is also necessary to avoid overflow checks in
+ // different places
+ min_x -= Real(2.1) * radius;
+ min_y -= Real(2.1) * radius;
+ min_z -= Real(2.1) * radius;
+ max_x += Real(2.1) * radius;
+ max_y += Real(2.1) * radius;
+ max_z += Real(2.1) * radius;
+
+ int num_xbins = std::ceil((max_x - min_x) * one_over_radius);
+ int num_ybins = std::ceil((max_y - min_y) * one_over_radius);
+ int num_zbins = std::ceil((max_z - min_z) * one_over_radius);
+ int num_bins = num_xbins * num_ybins * num_zbins;
+ char* grid = new char[num_bins];
+ memset(grid, 0, num_bins);
+
+ for(uint i = 0; i < atom_positions.size(); ++i) {
+
+ const geom::Vec3& pos = atom_positions[i];
+ int x_bin = (pos[0] - min_x) * one_over_radius;
+ int y_bin = (pos[1] - min_y) * one_over_radius;
+ int z_bin = (pos[2] - min_z) * one_over_radius;
+
+ // we're really crude here and simply set all 27 cubes with central
+ // cube defined by x_bin, y_bin and z_bin to one
+ for(int z = z_bin - 1; z <= z_bin + 1; ++z) {
+ for(int x = x_bin - 1; x <= x_bin + 1; ++x) {
+ for(int y = y_bin - 1; y <= y_bin + 1; ++y) {
+ grid[z*num_xbins*num_ybins + x*num_ybins + y] = 1;
+ }
+ }
+ }
+ }
+
+ // lets call flood fill for every layer along the z-axis from every
+ // corner in the x-y plane.
+ for(int z = 0; z < num_zbins; ++z) {
+ char* level = &grid[z*num_xbins*num_ybins];
+ FloodLevel(level, 0, 0, num_xbins, num_ybins, 0, 2);
+ FloodLevel(level, 0, num_ybins - 1, num_xbins, num_ybins, 0, 2);
+ FloodLevel(level, num_xbins - 1, 0, num_xbins, num_ybins, 0, 2);
+ FloodLevel(level, num_xbins - 1, num_ybins - 1, num_xbins, num_ybins, 0, 2);
+ }
+
+ // every cube in every layer that has currently value 1 that has a city-block
+ // distance below 3 to any cube with value 2 is considered to be in contact
+ // with the outer surface... lets set them to a value of three
+ for(int z = 0; z < num_zbins; ++z) {
+ char* level = &grid[z*num_xbins*num_ybins];
+ for(int x = 0; x < num_xbins; ++x) {
+ for(int y = 0; y < num_ybins; ++y) {
+ if(level[x*num_ybins + y] == 1) {
+ int x_from = std::max(0, x - 3);
+ int x_to = std::min(num_xbins-1, x + 3);
+ int y_from = std::max(0, y - 3);
+ int y_to = std::min(num_ybins-1, y + 3);
+ bool exposed = false;
+ for(int i = x_from; i <= x_to && !exposed; ++i) {
+ for(int j = y_from; j <= y_to; ++j) {
+ if(level[i*num_ybins + j] == 2) {
+ level[x*num_ybins + y] = 3;
+ exposed = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // all positions that lie in a cube with value 3 are considered to be exposed...
+ exposed_atom_positions.clear();
+ exposed_transfer_energies.clear();
+ for(uint i = 0; i < atom_positions.size(); ++i) {
+ const geom::Vec3& pos = atom_positions[i];
+ int x_bin = (pos[0] - min_x) * one_over_radius;
+ int y_bin = (pos[1] - min_y) * one_over_radius;
+ int z_bin = (pos[2] - min_z) * one_over_radius;
+ if(grid[z_bin*num_xbins*num_ybins + x_bin*num_ybins + y_bin] == 3) {
+ exposed_atom_positions.push_back(pos);
+ exposed_transfer_energies.push_back(transfer_energies[i]);
+ }
+ }
+
+ // cleanup
+ delete [] grid;
+}
+
+
+void ScanAxis(const std::vector<geom::Vec3>& atom_positions,
+ const std::vector<Real>& transfer_energies,
+ const geom::Vec3& axis,
+ Real& best_width, Real& best_center, Real& best_energy) {
+
+ int n_pos = atom_positions.size();
+
+ geom::Vec3 normalized_axis = geom::Normalize(axis);
+ std::vector<Real> pos_on_axis(n_pos);
+
+ for(int i = 0; i < n_pos; ++i) {
+ pos_on_axis[i] = geom::Dot(atom_positions[i], normalized_axis);
+ }
+
+ Real min_pos = pos_on_axis[0];
+ Real max_pos = min_pos;
+
+ for(int i = 1; i < n_pos; ++i) {
+ min_pos = std::min(min_pos, pos_on_axis[i]);
+ max_pos = std::max(max_pos, pos_on_axis[i]);
+ }
+
+ min_pos = std::floor(min_pos);
+ max_pos = std::ceil(max_pos);
+
+ int full_width = int(max_pos - min_pos) + 1;
+
+ //energies representing the energy profile along the axis
+ std::vector<Real> mapped_energies(full_width, 0.0);
+
+ for(int i = 0; i < n_pos; ++i) {
+ mapped_energies[int(pos_on_axis[i] - min_pos)] += transfer_energies[i];
+ }
+
+ best_width = 0;
+ best_center = 0;
+ best_energy = std::numeric_limits<Real>::max();
+
+ for(int window_width = 10; window_width <= 40; ++window_width) {
+
+ if(window_width > full_width) {
+ break;
+ }
+
+ Real energy=0.0;
+ for(int i = 0; i < window_width; ++i) {
+ energy += mapped_energies[i];
+ }
+
+ if(energy < best_energy) {
+ best_width = window_width;
+ best_center = min_pos + Real(0.5) * window_width + Real(0.5);
+ best_energy = energy;
+ }
+
+ for(int pos = 1; pos < full_width - window_width; ++pos) {
+ energy -= mapped_energies[pos-1];
+ energy += mapped_energies[pos+window_width-1];
+ if(energy < best_energy){
+ best_width = window_width;
+ best_center = min_pos + pos + Real(0.5) * window_width + Real(0.5);
+ best_energy = energy;
+ }
+ }
+ }
+}
+
+
+struct LMInput {
+ ost::mol::alg::FindMemParam mem_param;
+ geom::Transform initial_transform;
+ std::vector<geom::Vec3> exposed_atom_positions;
+ std::vector<Real> exposed_transfer_energies;
+};
+
+
+void SampleZ(const std::vector<geom::Vec3>& atom_pos,
+ const std::vector<Real>& transfer_energies,
+ const geom::Transform& initial_transform,
+ int n_solutions, std::list<LMInput>& top_solutions) {
+
+
+ std::vector<geom::Vec3> transformed_atom_pos(atom_pos.size());
+ for(uint at_idx = 0; at_idx < atom_pos.size(); ++at_idx) {
+ transformed_atom_pos[at_idx] = initial_transform.Apply(atom_pos[at_idx]);
+ }
+
+ std::vector<geom::Vec3> exposed_atom_positions;
+ std::vector<Real> exposed_transfer_energies;
+ GetExposedAtoms(transformed_atom_pos, transfer_energies,
+ exposed_atom_positions, exposed_transfer_energies);
+
+ std::vector<Real> tilt_angles;
+ std::vector<Real> rotation_angles;
+ for(int tilt_deg = 0; tilt_deg <= 45; tilt_deg += 5) {
+ if(tilt_deg == 0) {
+ tilt_angles.push_back(0.0);
+ rotation_angles.push_back(0.0);
+ }
+ else {
+ Real tilt_angle = Real(tilt_deg) / Real(180.) * Real(M_PI);
+ for(int angle_deg = 0; angle_deg < 360; angle_deg += 5) {
+ tilt_angles.push_back(tilt_angle);
+ rotation_angles.push_back(Real(angle_deg) / Real(180.) * Real(M_PI));
+ }
+ }
+ }
+
+ geom::Vec3 normalized_axis(0.0,0.0,1.0);
+ geom::Vec3 tilt_axis(1.0,0.0,0.0);
+
+ for(uint i = 0; i < tilt_angles.size(); ++i) {
+
+ Real tilt_angle = tilt_angles[i];
+ Real rotation_angle = rotation_angles[i];
+
+ geom::Vec3 tilted_axis = RotateAroundAxis(normalized_axis, tilt_axis,
+ tilt_angle);
+ geom::Vec3 scan_axis = RotateAroundAxis(tilted_axis, normalized_axis,
+ rotation_angle);
+
+ Real actual_width, actual_center, actual_energy;
+ ScanAxis(exposed_atom_positions, exposed_transfer_energies, scan_axis,
+ actual_width, actual_center, actual_energy);
+
+ if(static_cast<int>(top_solutions.size()) >= n_solutions &&
+ actual_energy > top_solutions.back().mem_param.energy) {
+ continue;
+ }
+
+ LMInput lm_input;
+ lm_input.mem_param.axis = normalized_axis;
+ lm_input.mem_param.tilt_axis = tilt_axis;
+ lm_input.mem_param.tilt = tilt_angle;
+ lm_input.mem_param.angle = rotation_angle;
+ lm_input.mem_param.width = actual_width;
+ lm_input.mem_param.pos = actual_center;
+ lm_input.mem_param.energy = actual_energy;
+ lm_input.initial_transform = initial_transform;
+ lm_input.exposed_atom_positions = exposed_atom_positions;
+ lm_input.exposed_transfer_energies = exposed_transfer_energies;
+
+ if(top_solutions.empty()) {
+ top_solutions.push_back(lm_input);
+ }
+ else {
+ bool added = false;
+ for(std::list<LMInput>::iterator sol_it = top_solutions.begin();
+ sol_it != top_solutions.end(); ++sol_it) {
+ if(sol_it->mem_param.energy > lm_input.mem_param.energy) {
+ top_solutions.insert(sol_it, lm_input);
+ added = true;
+ break;
+ }
+ }
+
+ if(!added) {
+ top_solutions.push_back(lm_input);
+ }
+
+ while(static_cast<int>(top_solutions.size()) > n_solutions) {
+ top_solutions.pop_back();
+ }
+ }
+ }
+}
+
+
+ost::mol::alg::FindMemParam GetFinalSolution(const std::list<LMInput>& top_solutions,
+ Real lambda) {
+
+ Real best_energy = std::numeric_limits<Real>::max();
+ std::list<LMInput>::const_iterator best_sol_it = top_solutions.begin();
+ Eigen::Matrix<Real, 4, 1> lm_parameters;
+ Eigen::Matrix<Real, 4, 1> best_lm_parameters;
+ LevenbergMarquardt<EnergyF, EnergyDF>::Results lm_result;
+
+ for(std::list<LMInput>::const_iterator sol_it = top_solutions.begin();
+ sol_it != top_solutions.end(); ++sol_it) {
+
+ Real offset = std::max(Real(20000.), std::abs(sol_it->mem_param.energy * 2));
+
+ EnergyF en_f(sol_it->exposed_atom_positions,
+ sol_it->exposed_transfer_energies,
+ lambda, offset,
+ sol_it->mem_param.axis,
+ sol_it->mem_param.tilt_axis);
+
+ lm_parameters(0,0) = sol_it->mem_param.tilt;
+ lm_parameters(1,0) = sol_it->mem_param.angle;
+ lm_parameters(2,0) = sol_it->mem_param.width;
+ lm_parameters(3,0) = sol_it->mem_param.pos;
+
+ LevenbergMarquardt<EnergyF,EnergyDF> lm(en_f);
+ lm_result = lm.minimize(&lm_parameters);
+
+ Real minimized_energy = en_f(lm_parameters)(0, 0) - en_f.offset;
+
+ if(minimized_energy < best_energy) {
+ best_energy = minimized_energy;
+ best_sol_it = sol_it;
+ best_lm_parameters = lm_parameters;
+ }
+ }
+
+ ost::mol::alg::FindMemParam mem_param = best_sol_it->mem_param;
+ mem_param.energy = best_energy;
+ mem_param.tilt = best_lm_parameters(0,0);
+ mem_param.angle = best_lm_parameters(1,0);
+ mem_param.width = best_lm_parameters(2,0);
+ mem_param.pos = best_lm_parameters(3,0);
+
+ // the solution is still relative to the initial transform that has
+ // been applied when calling the SampleZ funtion!
+ geom::Transform t = best_sol_it->initial_transform;
+ mem_param.tilt_axis = t.ApplyInverse(mem_param.tilt_axis);
+ mem_param.axis = t.ApplyInverse(mem_param.axis);
+
+ return mem_param;
+}
+
+
+ost::mol::EntityHandle CreateMembraneRepresentation(
+ const std::vector<geom::Vec3>& atom_positions,
+ const ost::mol::alg::FindMemParam& param,
+ Real membrane_margin = 15,
+ Real delta = 2.0) {
+
+ // let's first construct two planes defining the membrane
+ geom::Vec3 membrane_axis = param.GetMembraneAxis();
+ geom::Vec3 one = param.pos * membrane_axis +
+ membrane_axis * param.width / 2;
+ geom::Vec3 two = param.pos * membrane_axis -
+ membrane_axis * param.width / 2;
+ geom::Plane plane_one = geom::Plane(one, membrane_axis);
+ geom::Plane plane_two = geom::Plane(two, membrane_axis);
+
+ // let's find all positions that are somehow close to those planes
+ geom::Vec3List close_pos;
+ geom::Vec3List close_pos_one;
+ geom::Vec3List close_pos_two;
+
+ for(uint i = 0; i < atom_positions.size(); ++i) {
+
+ Real d1 = geom::Distance(plane_one, atom_positions[i]);
+ Real d2 = geom::Distance(plane_two, atom_positions[i]);
+
+ if(d1 < Real(3.)) {
+ close_pos_one.push_back(atom_positions[i]);
+ }
+
+ if(d2 < Real(3.)) {
+ close_pos_two.push_back(atom_positions[i]);
+ }
+
+ if(d1 < Real(3.) || d2 < Real(3.)) {
+ close_pos.push_back(atom_positions[i]);
+ }
+ }
+
+ // the geometric center of the close pos vector in combination with the
+ // membrane axis define the central "line" of the disks that will represent
+ // the membrane
+ geom::Vec3 center_pos = close_pos.GetCenter();
+ geom::Line3 center_line = geom::Line3(center_pos, center_pos + membrane_axis);
+
+ // the final radius of the "disks" is based on the maximal distance of any
+ // position in close_pos to the center_line plus the membrane_margin
+
+ Real max_d_to_center_line = 0;
+ for(uint i = 0; i < close_pos.size(); ++i) {
+ Real d = geom::Distance(center_line, close_pos[i]);
+ max_d_to_center_line = std::max(max_d_to_center_line, d);
+ }
+
+ Real disk_radius = max_d_to_center_line + membrane_margin;
+ int num_sampling_points = (Real(2.) * disk_radius) / delta;
+
+ // reassign the top and bottom positions, that have been only arbitrary
+ // points on the membrane planes
+ one = geom::IntersectionPoint(center_line, plane_one);
+ two = geom::IntersectionPoint(center_line, plane_two);
+
+ // find a pair of perpendicular vectors, that are on the plane
+ geom::Vec3 arbitrary_vec(1.0, 0.0, 0.0);
+ if(geom::Angle(membrane_axis, arbitrary_vec) < 0.1) {
+ // parallel is not cool in this case
+ arbitrary_vec = geom::Vec3(0.0, 1.0, 0.0);
+ }
+ geom::Vec3 plane_x = geom::Normalize(geom::Cross(membrane_axis, arbitrary_vec));
+ geom::Vec3 plane_y = geom::Normalize(geom::Cross(membrane_axis, plane_x));
+
+ // final representing positions come in here
+ std::vector<geom::Vec3> final_pos;
+
+ // do plane one
+ geom::Vec3 origin = one - delta * num_sampling_points * 0.5 * plane_x -
+ delta * num_sampling_points * 0.5 * plane_y;
+
+ for(int i = 0; i < num_sampling_points; ++i) {
+ for(int j = 0; j < num_sampling_points; ++j) {
+ geom::Vec3 pos = origin + i*delta*plane_x + j*delta*plane_y;
+ if(geom::Distance(pos, one) < disk_radius) {
+ bool far_far_away = true;
+ // this is slow...
+ for(uint k = 0; k < close_pos_one.size(); ++k) {
+ if(geom::Length2(pos - close_pos_one[k]) < Real(16.)) {
+ far_far_away = false;
+ break;
+ }
+ }
+ if(far_far_away) {
+ final_pos.push_back(pos);
+ }
+ }
+ }
+ }
+
+ // do plane two
+ origin = two - delta * num_sampling_points * 0.5 * plane_x -
+ delta * num_sampling_points * 0.5 * plane_y;
+
+ for(int i = 0; i < num_sampling_points; ++i) {
+ for(int j = 0; j < num_sampling_points; ++j) {
+ geom::Vec3 pos = origin + i*delta*plane_x + j*delta*plane_y;
+ if(geom::Distance(pos, two) < disk_radius) {
+ bool far_far_away = true;
+ // this is slow...
+ for(uint k = 0; k < close_pos_two.size(); ++k) {
+ if(geom::Length2(pos - close_pos_two[k]) < Real(16.)) {
+ far_far_away = false;
+ break;
+ }
+ }
+ if(far_far_away) {
+ final_pos.push_back(pos);
+ }
+ }
+ }
+ }
+
+ // create hacky entity that contains membrane representing positions and
+ // return
+ ost::mol::EntityHandle membrane_ent = ost::mol::CreateEntity();
+ ost::mol::XCSEditor ed = membrane_ent.EditXCS();
+
+ ost::mol::ChainHandle chain = ed.InsertChain("M");
+ ost::mol::ResidueHandle res = ed.AppendResidue(chain, "MEM");
+ String atom_names = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+ uint atom_name_idx = 0;
+ uint atom_name_secondary_idx = 0;
+
+ for(uint i = 0; i < final_pos.size(); ++i) {
+
+ if(atom_name_secondary_idx == atom_names.size()) {
+ ++atom_name_idx;
+ atom_name_secondary_idx = 0;
+ }
+ if(atom_name_idx == atom_names.size()) {
+ res = ed.AppendResidue(chain, "MEM");
+ atom_name_idx = 0;
+ atom_name_secondary_idx = 0;
+ }
+
+ String atom_name = "--";
+ atom_name[0] = atom_names[atom_name_idx];
+ atom_name[1] = atom_names[atom_name_secondary_idx];
+
+ ed.InsertAtom(res, atom_name, final_pos[i]);
+ ++atom_name_secondary_idx;
+ }
+
+ return membrane_ent;
+}
+
+
+} // anon namespace
+
+
+namespace ost{ namespace mol{ namespace alg{
+
+geom::Vec3 FindMemParam::GetMembraneAxis() const {
+
+ geom::Vec3 result = RotateAroundAxis(axis,tilt_axis,tilt);
+ result = RotateAroundAxis(result,axis,angle);
+ return result;
+}
+
+
+FindMemParam FindMembrane(ost::mol::EntityView& ent,
+ bool assign_membrane_representation, bool fast) {
+
+ ost::mol::EntityView peptide_view = ent.Select("peptide=true and ele!=H");
+ Accessibility(peptide_view);
+
+ std::vector<geom::Vec3> atom_pos;
+ std::vector<Real> transfer_energies;
+
+ atom_pos.reserve(peptide_view.GetAtomCount());
+ transfer_energies.reserve(peptide_view.GetAtomCount());
+
+ ost::mol::AtomViewList atoms = peptide_view.GetAtomList();
+ String stupid_string("S_N_O_C");
+
+ for(ost::mol::AtomViewList::iterator it = atoms.begin();
+ it != atoms.end(); ++it) {
+
+ if(!it->HasProp("asaAtom")) {
+ continue;
+ }
+
+ String element = it->GetElement();
+ if(stupid_string.find(element) == std::string::npos) {
+ continue;
+ }
+
+ Real asa = it->GetFloatProp("asaAtom");
+ atom_pos.push_back(it->GetPos());
+
+ if(element == "S") {
+ transfer_energies.push_back(asa * Real(10.0));
+ }
+ else if(element == "N") {
+ transfer_energies.push_back(asa * Real(53.0));
+ }
+ else if(element == "O") {
+ transfer_energies.push_back(asa * Real(57.0));
+ }
+ else if(element == "C") {
+ // check whether we find a double bond to distinguish between
+ // hibridization states
+ bool assigned_energy = false;
+ ost::mol::BondHandleList bond_list = it->GetBondList();
+ for(ost::mol::BondHandleList::iterator bond_it = bond_list.begin();
+ bond_it != bond_list.end(); ++bond_it){
+ unsigned char bond_order = bond_it->GetBondOrder();
+ if(bond_order > '1'){
+ transfer_energies.push_back(asa * Real(-19.0));
+ assigned_energy = true;
+ break;;
+ }
+ }
+ if(!assigned_energy) {
+ transfer_energies.push_back(asa * Real(-22.6));
+ }
+ }
+ }
+
+ if(atom_pos.size() < 10) {
+ throw ost::Error("Cannot detect membrane with such a low number "
+ "of heavy atoms!");
+ }
+
+ // we always optimizer along the z-axis.
+ // We therefore have to transform the positions. We use a rotation
+ // around the z-axis with subsequent rotation around the x-axis for this task
+ std::vector<geom::Transform> transformations;
+ int n_euler_angles = 3;
+ int n_transformations = n_euler_angles * n_euler_angles * n_euler_angles;
+ std::vector<Real> euler_angles(n_euler_angles);
+ euler_angles[0] = 0.0;
+ euler_angles[1] = M_PI/3;
+ euler_angles[2] = 2*M_PI/3;
+
+ for(int i = 0; i < n_euler_angles; ++i) {
+ for(int j = 0; j < n_euler_angles; ++j) {
+ for(int k = 0; k < n_euler_angles; ++k) {
+ geom::Mat3 rot_matrix = geom::EulerTransformation(euler_angles[i],
+ euler_angles[j],
+ euler_angles[k]);
+ geom::Transform transform;
+ transform.SetRot(rot_matrix);
+ transformations.push_back(transform);
+ }
+ }
+ }
+
+ // lets use the generated transforms to search for initial solutions that can
+ // then be fed into a final minimization...
+ std::list<LMInput> top_solutions;
+ int n_initial_solutions = fast ? 1 : 20;
+
+ for(int transformation_idx = 0; transformation_idx < n_transformations;
+ ++transformation_idx) {
+ SampleZ(atom_pos, transfer_energies, transformations[transformation_idx],
+ n_initial_solutions, top_solutions);
+ }
+
+ // Perform the final minimization and return the best solution.
+ // please note, that the returned solution is transformed back in order
+ // to match the initial atom positions
+ FindMemParam final_solution = GetFinalSolution(top_solutions, 0.9);
+
+ if(assign_membrane_representation) {
+ final_solution.membrane_representation = CreateMembraneRepresentation(
+ atom_pos,
+ final_solution);
+ }
+
+ return final_solution;
+}
+
+
+FindMemParam FindMembrane(ost::mol::EntityHandle& ent,
+ bool assign_membrane_representation,
+ bool fast) {
+
+ ost::mol::EntityView ent_view = ent.CreateFullView();
+ return FindMembrane(ent_view, assign_membrane_representation, fast);
+}
+
+}}} // ns
diff --git a/modules/mol/alg/src/find_membrane.hh b/modules/mol/alg/src/find_membrane.hh
new file mode 100644
index 0000000000000000000000000000000000000000..962ea68fc09680304c2da095f78900eff5483225
--- /dev/null
+++ b/modules/mol/alg/src/find_membrane.hh
@@ -0,0 +1,31 @@
+#include <ost/mol/mol.hh>
+
+#include <ost/geom/geom.hh>
+#include <ost/io/binary_data_source.hh>
+#include <ost/io/binary_data_sink.hh>
+
+namespace ost { namespace mol{ namespace alg{
+
+struct FindMemParam{
+ FindMemParam() { }
+
+ geom::Vec3 GetMembraneAxis() const;
+ geom::Vec3 axis;
+ geom::Vec3 tilt_axis;
+ Real tilt;
+ Real angle;
+ Real width;
+ Real pos;
+ Real energy;
+ ost::mol::EntityHandle membrane_representation;
+};
+
+FindMemParam FindMembrane(ost::mol::EntityHandle& ent,
+ bool assign_membrane_representation,
+ bool fast);
+
+FindMemParam FindMembrane(ost::mol::EntityView& ent,
+ bool assign_membrane_representation,
+ bool fast);
+
+}}} // ns
diff --git a/modules/mol/base/doc/entity.rst b/modules/mol/base/doc/entity.rst
index 653625756a1b38c696a50da10840078fb4d74c58..3b0aa5c9fdcd51cdb563922b9fb384f8ec046736 100644
--- a/modules/mol/base/doc/entity.rst
+++ b/modules/mol/base/doc/entity.rst
@@ -118,7 +118,13 @@ The Handle Classes
an enabled ``USE_NUMPY`` flag (see :ref:`here <cmake-flags>` for details).
:type: :class:`numpy.array`
-
+
+ .. attribute:: valid
+
+ Validity of handle.
+
+ :type: bool
+
.. method:: GetName()
:returns: Name associated to this entity.
@@ -327,6 +333,10 @@ The Handle Classes
:type radius: float
:returns: :class:`AtomHandleList` (list of :class:`AtomHandle`)
+
+ .. method:: IsValid()
+
+ See :attr:`valid`
.. class:: ChainHandle
@@ -424,6 +434,12 @@ The Handle Classes
:meth:`GetCenterOfAtoms`.
:type: :class:`~ost.geom.Vec3`
+
+ .. attribute:: valid
+
+ Validity of handle.
+
+ :type: bool
.. method:: FindResidue(res_num)
@@ -469,6 +485,10 @@ The Handle Classes
See :attr:`description`
+ .. method:: IsValid()
+
+ See :attr:`valid`
+
.. class:: ResidueHandle
The residue is either used to represent complete molecules or building blocks
@@ -620,6 +640,43 @@ The Handle Classes
Residue index (starting at 0) within chain.
+ .. attribute:: central_atom
+
+ Central atom used for rendering traces. For peptides, this is usually
+ the CA atom. For nucleotides, this is usually the P atom.
+
+ :type: :class:`AtomHandle`
+
+ .. attribute:: central_normal
+
+ Normal computed for :attr:`central_atom`. Only defined for peptides and
+ nucleotides if all required atoms available. Otherwise, the (1,0,0) vector
+ is returned.
+
+ :type: :class:`~ost.geom.Vec3`
+
+ .. attribute:: valid
+
+ Validity of handle.
+
+ :type: bool
+
+ .. attribute:: next
+
+ Residue after this one in the same chain. Invalid handle returned if there
+ is no next residue. Residues are ordered as in :attr:`ChainHandle.residues`
+ independently on whether they are connected or not (see
+ :func:`InSequence` to check for connected residues).
+
+ :type: :class:`ResidueHandle`
+
+ .. attribute:: prev
+
+ Residue before this one in the same chain. Otherwise same behaviour as
+ :attr:`next`.
+
+ :type: :class:`ResidueHandle`
+
.. method:: FindAtom(atom_name)
Get atom by atom name. See also :attr:`atoms`
@@ -666,7 +723,20 @@ The Handle Classes
.. method:: GetIndex()
See :attr:`index`
+
+ .. method:: GetCentralAtom()
+ SetCentralAtom()
+
+ See :attr:`central_atom`
+
+ .. method:: GetCentralNormal()
+
+ See :attr:`central_normal`
+
+ .. method:: IsValid()
+ See :attr:`valid`
+
.. class:: AtomHandle
@@ -732,6 +802,7 @@ The Handle Classes
The atom's occupancy in the range 0 to 1. Read/write. Also available as
:meth:`GetOccupancy`, :meth:`SetOccupancy`.
+
:type: float
.. attribute:: b_factor
@@ -782,6 +853,12 @@ The Handle Classes
:type: int
+ .. attribute:: valid
+
+ Validity of handle.
+
+ :type: bool
+
.. method:: FindBondToAtom(other_atom)
Finds and returns the bond formed between this atom and `other_atom`. If no
@@ -899,8 +976,7 @@ The Handle Classes
.. method:: IsValid()
See :attr:`valid`
-
- :rtype: bool
+
The View Classes
--------------------------------------------------------------------------------
@@ -981,6 +1057,12 @@ The View Classes
:type: :class:`EntityHandle`
+ .. attribute:: valid
+
+ Validity of view.
+
+ :type: bool
+
.. method:: GetName()
:returns: :func:`~EntityHandle.GetName` of entity :attr:`handle`.
@@ -1199,7 +1281,8 @@ The View Classes
.. method:: GetBondCount()
Get number of bonds
- :rtype: int
+
+ :rtype: :class:`int`
.. method:: GetBondList()
@@ -1241,6 +1324,10 @@ The View Classes
See :attr:`atom_count`
+ .. method:: IsValid()
+
+ See :attr:`valid`
+
.. class:: ChainView
A view representation of a :class:`ChainHandle`. Mostly, the same
@@ -1350,7 +1437,7 @@ The View Classes
.. attribute:: valid
- Validity of handle.
+ Validity of view.
:type: bool
@@ -1752,6 +1839,28 @@ Other Entity-Related Functions
:returns: :class:`EntityHandle`
+.. function:: InSequence(res, res_next)
+
+ :return: True, if both *res* and *res_next* are :attr:`~ResidueHandle.valid`,
+ *res_next* is the residue following *res* (see
+ :attr:`ResidueHandle.next`), both residues are linking (i.e.
+ :attr:`~ChemClass.IsPeptideLinking` or
+ :attr:`~ChemClass.IsNucleotideLinking`) and are connected by an
+ appropriate bond.
+ :rtype: :class:`bool`
+ :param res: First residue to check.
+ :type res: :class:`ResidueHandle`
+ :param res_next: Second residue to check.
+ :type res_next: :class:`ResidueHandle`
+
+.. function:: BondExists(atom_a, atom_b)
+
+ :return: True, if *atom_a* and *atom_b* are connected by a bond.
+ :rtype: :class:`bool`
+ :param atom_a: First atom to check.
+ :type atom_a: :class:`AtomHandle`
+ :param atom_b: Second atom to check.
+ :type atom_b: :class:`AtomHandle`
Residue Numbering
--------------------------------------------------------------------------------
diff --git a/modules/mol/base/src/impl/residue_impl.cc b/modules/mol/base/src/impl/residue_impl.cc
index fe1e7f60b0448d79e661fabc2b4cbe2630cd58ae..014f4069f5cae65e808d25151db5cee4ef10f213 100644
--- a/modules/mol/base/src/impl/residue_impl.cc
+++ b/modules/mol/base/src/impl/residue_impl.cc
@@ -212,12 +212,12 @@ AtomImplPtr ResidueImpl::GetCentralAtom() const
for (AtomImplList::const_iterator it=atom_list_.begin();
it!=atom_list_.end();++it) {
if((*it)->Name()=="P") return *it;
- }
+ }
} else if (chem_class_.IsPeptideLinking()) {
for (AtomImplList::const_iterator it=atom_list_.begin();
it!=atom_list_.end();++it) {
if((*it)->Name()=="CA") return *it;
- }
+ }
}
return AtomImplPtr();
@@ -266,18 +266,21 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const
geom::Vec3 nrvo(1,0,0);
if (chem_class_.IsPeptideLinking()) {
AtomImplPtr a1 = FindAtom("C");
- AtomImplPtr a2 = FindAtom("O");
+ AtomImplPtr a2 = FindAtom("O");
if(a1 && a2) {
nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos());
} else {
a1 = FindAtom("CB");
- a2 = FindAtom("CA");
+ a2 = FindAtom("CA");
if(a1 && a2) {
nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos());
} else {
- geom::Vec3 v0=GetCentralAtom()->TransformedPos();
- nrvo=geom::Cross(geom::Normalize(v0),
- geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1])));
+ AtomImplPtr a0 = GetCentralAtom();
+ if (a0) {
+ geom::Vec3 v0 = a0->TransformedPos();
+ nrvo = geom::Cross(geom::Normalize(v0),
+ geom::Normalize(geom::Vec3(-v0[2], v0[0], v0[1])));
+ }
LOG_VERBOSE("warning: could not find atoms for proper central normal calculation");
}
}
@@ -288,9 +291,12 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const
if(a1 && a2 && a3) {
nrvo = geom::Normalize(a1->TransformedPos()-(a2->TransformedPos()+a3->TransformedPos())*.5);
} else {
- geom::Vec3 v0=GetCentralAtom()->TransformedPos();
- nrvo=geom::Cross(geom::Normalize(v0),
- geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1])));
+ AtomImplPtr a0 = GetCentralAtom();
+ if (a0) {
+ geom::Vec3 v0 = a0->TransformedPos();
+ nrvo = geom::Cross(geom::Normalize(v0),
+ geom::Normalize(geom::Vec3(-v0[2], v0[0], v0[1])));
+ }
LOG_VERBOSE("warning: could not find atoms for proper central normal calculation");
}
}
diff --git a/modules/mol/base/tests/test_residue.cc b/modules/mol/base/tests/test_residue.cc
index 89e8da33914110dc1d50150fb9dc19e2ab3dd50e..d828544b7eba1ef47d5c995b4341b438e68e8c9d 100644
--- a/modules/mol/base/tests/test_residue.cc
+++ b/modules/mol/base/tests/test_residue.cc
@@ -125,4 +125,92 @@ BOOST_AUTO_TEST_CASE(rename_res)
BOOST_CHECK_EQUAL(rA2B.GetName(), "B");
}
+BOOST_AUTO_TEST_CASE(test_centralatom)
+{
+ // COOK UP ENTITY FOR TEST
+ EntityHandle eh = CreateEntity();
+ XCSEditor e = eh.EditXCS();
+ ChainHandle ch = e.InsertChain("A");
+ // decent peptide with all entries
+ ResidueHandle rp1 = e.AppendResidue(ch, "A");
+ e.InsertAtom(rp1, "CA", geom::Vec3(2, 0, 0));
+ e.InsertAtom(rp1, "CB", geom::Vec3(1, 0, 0));
+ e.InsertAtom(rp1, "C", geom::Vec3(0, 0, 0));
+ e.InsertAtom(rp1, "O", geom::Vec3(0, 1, 0));
+ rp1.SetChemClass(ChemClass(ChemClass::PEPTIDE_LINKING));
+ // weird peptide with only CA and CB
+ ResidueHandle rp2 = e.AppendResidue(ch, "B");
+ e.InsertAtom(rp2, "CA", geom::Vec3(3, 0, 0));
+ e.InsertAtom(rp2, "CB", geom::Vec3(3, 1, 0));
+ rp2.SetChemClass(ChemClass(ChemClass::PEPTIDE_LINKING));
+ // CA-only peptide
+ ResidueHandle rp3 = e.AppendResidue(ch, "C");
+ e.InsertAtom(rp3, "CA", geom::Vec3(4, 0, 0));
+ rp3.SetChemClass(ChemClass(ChemClass::PEPTIDE_LINKING));
+ // peptide with custom atoms
+ ResidueHandle rp4 = e.AppendResidue(ch, "D");
+ AtomHandle rp4_ax = e.InsertAtom(rp4, "XX", geom::Vec3(5, 0, 0));
+ rp4.SetChemClass(ChemClass(ChemClass::PEPTIDE_LINKING));
+ // nucleotide with all needed entries
+ ResidueHandle rn1 = e.AppendResidue(ch, "E");
+ e.InsertAtom(rn1, "P", geom::Vec3(6, 0, 0));
+ e.InsertAtom(rn1, "OP1", geom::Vec3(6, 0.5, 0));
+ e.InsertAtom(rn1, "OP2", geom::Vec3(6, 1.5, 0));
+ rn1.SetChemClass(ChemClass(ChemClass::DNA_LINKING));
+ // nucleotide with only P
+ ResidueHandle rn2 = e.AppendResidue(ch, "F");
+ e.InsertAtom(rn2, "P", geom::Vec3(7, 0, 0));
+ rn2.SetChemClass(ChemClass(ChemClass::DNA_LINKING));
+ // nucleotide with custom atoms
+ ResidueHandle rn3 = e.AppendResidue(ch, "G");
+ AtomHandle rn3_ax = e.InsertAtom(rn3, "XX", geom::Vec3(8, 0, 0));
+ rn3.SetChemClass(ChemClass(ChemClass::DNA_LINKING));
+ // unknown chem class
+ ResidueHandle ru = e.AppendResidue(ch, "H");
+ e.InsertAtom(ru, "P", geom::Vec3(9, 0, 0));
+ e.InsertAtom(ru, "CA", geom::Vec3(9, 1, 0));
+ AtomHandle ru_ax = e.InsertAtom(ru, "XX", geom::Vec3(9, 2, 0));
+ ru.SetChemClass(ChemClass(ChemClass::UNKNOWN));
+
+ // CHECK CENTRAL ATOMS
+ BOOST_CHECK(rp1.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rp1.GetCentralAtom().GetQualifiedName(), "A.A1.CA");
+ BOOST_CHECK(rp2.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rp2.GetCentralAtom().GetQualifiedName(), "A.B2.CA");
+ BOOST_CHECK(rp3.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rp3.GetCentralAtom().GetQualifiedName(), "A.C3.CA");
+ BOOST_CHECK(!rp4.GetCentralAtom().IsValid());
+ BOOST_CHECK(rn1.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rn1.GetCentralAtom().GetQualifiedName(), "A.E5.P");
+ BOOST_CHECK(rn2.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rn2.GetCentralAtom().GetQualifiedName(), "A.F6.P");
+ BOOST_CHECK(!rn3.GetCentralAtom().IsValid());
+ BOOST_CHECK(!ru.GetCentralAtom().IsValid());
+
+ // CHECK NORMALS
+ BOOST_CHECK_EQUAL(rp1.GetCentralNormal(), geom::Vec3(0, 1, 0));
+ BOOST_CHECK_EQUAL(rp2.GetCentralNormal(), geom::Vec3(0, -1, 0));
+ BOOST_CHECK_EQUAL(rp3.GetCentralNormal(), geom::Vec3(0, 0, 1));
+ BOOST_CHECK_EQUAL(rp4.GetCentralNormal(), geom::Vec3(1, 0, 0));
+ BOOST_CHECK_EQUAL(rn1.GetCentralNormal(), geom::Vec3(0, -1, 0));
+ BOOST_CHECK_EQUAL(rn2.GetCentralNormal(), geom::Vec3(0, 0, 1));
+ BOOST_CHECK_EQUAL(rn3.GetCentralNormal(), geom::Vec3(1, 0, 0));
+ BOOST_CHECK_EQUAL(ru.GetCentralNormal(), geom::Vec3(1, 0, 0));
+
+ // CHECK SETTING CENTRAL ATOMS
+ rp4.SetCentralAtom(rp4_ax);
+ BOOST_CHECK(rp4.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rp4.GetCentralAtom().GetQualifiedName(), "A.D4.XX");
+ BOOST_CHECK_EQUAL(rp4.GetCentralNormal(), geom::Vec3(0, 0, 1));
+ rn3.SetCentralAtom(rn3_ax);
+ BOOST_CHECK(rn3.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(rn3.GetCentralAtom().GetQualifiedName(), "A.G7.XX");
+ BOOST_CHECK_EQUAL(rn3.GetCentralNormal(), geom::Vec3(0, 0, 1));
+ ru.SetCentralAtom(ru_ax);
+ BOOST_CHECK(ru.GetCentralAtom().IsValid());
+ BOOST_CHECK_EQUAL(ru.GetCentralAtom().GetQualifiedName(), "A.H8.XX");
+ // no normal for unknown residues
+ BOOST_CHECK_EQUAL(ru.GetCentralNormal(), geom::Vec3(1, 0, 0));
+}
+
BOOST_AUTO_TEST_SUITE_END();
diff --git a/modules/mol/mm/src/simulation.cc b/modules/mol/mm/src/simulation.cc
index 2a23ddb9a629a112b3301c35dd6756ad2ba08b92..dea985213bf993a1a49f0ec06c17f7b34d4e36fe 100644
--- a/modules/mol/mm/src/simulation.cc
+++ b/modules/mol/mm/src/simulation.cc
@@ -27,10 +27,10 @@ Simulation::Simulation(const ost::mol::EntityHandle& handle,
//note, that ent_ will be "completed" inside this function!
//(hydrogens and shit)
-
- ent_ = handle.Copy();
- TopologyPtr top = TopologyCreator::Create(ent_,settings);
- this->Init(top, settings);
+
+ ost::mol::EntityHandle ent = handle.Copy();
+ TopologyPtr top = TopologyCreator::Create(ent,settings);
+ this->Init(top, ent, settings);
}
Simulation::Simulation(const TopologyPtr top,
@@ -40,21 +40,16 @@ Simulation::Simulation(const TopologyPtr top,
if(static_cast<uint>(handle.GetAtomCount()) != top->GetNumParticles()){
throw ost::Error("Number of atoms in entity must be consistent with number of particles in topology!");
}
- ent_ = handle.Copy();
- this->Init(top, settings);
+ ost::mol::EntityHandle ent = handle.Copy();
+ this->Init(top, ent, settings);
}
void Simulation::Save(const String& filename){
+
std::ofstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSink ds(stream);
+
ds << *top_;
- geom::Vec3List positions = this->GetPositions(false,false);
- for(geom::Vec3List::iterator i = positions.begin();
- i != positions.end(); ++i){
- ds & (*i)[0];
- ds & (*i)[1];
- ds & (*i)[2];
- }
uint num_chains;
uint num_residues;
@@ -96,15 +91,11 @@ void Simulation::Save(const String& filename){
k != atom_list.end(); ++k){
atom_name = k->GetName();
atom_element = k->GetElement();
- geom::Vec3 pos = k->GetPos();
bfac = k->GetBFactor();
occ = k->GetOccupancy();
is_hetatm = k->IsHetAtom();
ds & atom_name;
ds & atom_element;
- ds & pos[0];
- ds & pos[1];
- ds & pos[2];
ds & bfac;
ds & occ;
ds & is_hetatm;
@@ -112,19 +103,18 @@ void Simulation::Save(const String& filename){
}
}
- ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
ost::mol::AtomHandleList bonded_atoms;
std::map<long,int> atom_indices;
int actual_index = 0;
- for(ost::mol::AtomHandleList::const_iterator i = atom_list.begin(), e = atom_list.end();
- i != e; ++i){
+ for(ost::mol::AtomHandleList::const_iterator i = atom_list_.begin(),
+ e = atom_list_.end(); i != e; ++i){
atom_indices[i->GetHashCode()] = actual_index;
++actual_index;
}
- for(ost::mol::AtomHandleList::iterator i = atom_list.begin();
- i != atom_list.end(); ++i){
+ for(ost::mol::AtomHandleList::iterator i = atom_list_.begin();
+ i != atom_list_.end(); ++i){
bonded_atoms = i->GetBondPartners();
num_bonded_atoms = bonded_atoms.size();
ds & num_bonded_atoms;
@@ -139,6 +129,7 @@ void Simulation::Save(const String& filename){
}
SimulationPtr Simulation::Load(const String& filename, SettingsPtr settings){
+
if (!boost::filesystem::exists(filename)) {
std::stringstream ss;
ss << "Could not open simulation File '"
@@ -146,71 +137,17 @@ SimulationPtr Simulation::Load(const String& filename, SettingsPtr settings){
throw ost::io::IOException(ss.str());
}
- SimulationPtr sim_ptr(new Simulation);
-
std::ifstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSource ds(stream);
- TopologyPtr top_p(new Topology);
- ds >> *top_p;
-
- sim_ptr->top_ = top_p;
-
- sim_ptr->system_ = SystemCreator::Create(sim_ptr->top_,settings,
- sim_ptr->system_force_mapper_);
-
- sim_ptr->integrator_ = settings->integrator;
-
- OpenMM::Platform::loadPluginsFromDirectory (settings->openmm_plugin_directory);
- OpenMM::Platform::loadPluginsFromDirectory (settings->custom_plugin_directory);
- OpenMM::Platform* platform;
-
- switch(settings->platform){
- case Reference:{
- platform = &OpenMM::Platform::getPlatformByName("Reference");
- break;
- }
- case OpenCL:{
- platform = &OpenMM::Platform::getPlatformByName("OpenCL");
- break;
- }
- case CUDA:{
- platform = &OpenMM::Platform::getPlatformByName("CUDA");
- break;
- }
- case CPU:{
- platform = &OpenMM::Platform::getPlatformByName("CPU");
- break;
- }
- default:{
- throw ost::Error("Invalid Platform when Loading simulation!");
- }
- }
-
- sim_ptr->context_ = ContextPtr(new OpenMM::Context(*(sim_ptr->system_),
- *(sim_ptr->integrator_),
- *platform));
- std::vector<OpenMM::Vec3> positions;
- OpenMM::Vec3 open_mm_vec;
- Real a,b,c;
- for(int i = 0; i < sim_ptr->system_->getNumParticles(); ++i){
- ds & a;
- ds & b;
- ds & c;
- open_mm_vec[0] = a;
- open_mm_vec[1] = b;
- open_mm_vec[2] = c;
- positions.push_back(open_mm_vec);
- }
- sim_ptr->context_->setPositions(positions);
+ SimulationPtr sim_ptr(new Simulation);
+ TopologyPtr top(new Topology);
+ ds >> *top;
uint num_chains;
uint num_residues;
uint num_atoms;
uint num_bonded_atoms;
- Real x_pos;
- Real y_pos;
- Real z_pos;
Real bfac;
Real occ;
bool is_hetatm;
@@ -239,29 +176,29 @@ SimulationPtr Simulation::Load(const String& filename, SettingsPtr settings){
for(uint k = 0; k < num_atoms; ++k){
ds & atom_name;
ds & atom_element;
- ds & x_pos;
- ds & y_pos;
- ds & z_pos;
ds & bfac;
ds & occ;
ds & is_hetatm;
- geom::Vec3 pos(x_pos,y_pos,z_pos);
- ed.InsertAtom(res,atom_name,pos,atom_element,occ,bfac,is_hetatm);
+ ed.InsertAtom(res, atom_name, geom::Vec3(0.0,0.0,0.0),
+ atom_element, occ, bfac, is_hetatm);
}
}
}
- ost::mol::AtomHandleList atom_list = ent.GetAtomList();
- for(uint i = 0; i < atom_list.size(); ++i){
+
+ sim_ptr->Init(top, ent, settings);
+
+ for(uint i = 0; i < sim_ptr->atom_list_.size(); ++i){
ds & num_bonded_atoms;
for(uint j = 0; j < num_bonded_atoms; ++j){
ds & atom_index;
- ed.Connect(atom_list[i],atom_list[atom_index]);
+ ed.Connect(sim_ptr->atom_list_[i], sim_ptr->atom_list_[atom_index]);
}
}
- sim_ptr->ent_ = ent;
-
+ // also loads the positions that have been set in the context
+ // they get mapped over to the attached entity
sim_ptr->context_->loadCheckpoint(stream);
+ sim_ptr->UpdatePositions();
return sim_ptr;
}
@@ -297,10 +234,12 @@ void Simulation::EnsurePluginsLoaded(const String& plugin_path) {
void Simulation::Init(const TopologyPtr top,
+ const ost::mol::EntityHandle& ent,
const SettingsPtr settings){
-
top_ = top;
+ ent_ = ent;
+ atom_list_ = ent_.GetAtomList();
if(!settings->integrator){
//user did not specify an integrator, so let's just use a standard integrator
@@ -358,12 +297,11 @@ void Simulation::Init(const TopologyPtr top,
context_ = ContextPtr(new OpenMM::Context(*system_,*integrator_,*platform,context_properties));
- ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
std::vector<OpenMM::Vec3> positions;
geom::Vec3 ost_vec;
OpenMM::Vec3 open_mm_vec;
- for(ost::mol::AtomHandleList::iterator i = atom_list.begin();
- i!=atom_list.end();++i){
+ for(ost::mol::AtomHandleList::iterator i = atom_list_.begin();
+ i!=atom_list_.end();++i){
ost_vec = i->GetPos();
open_mm_vec[0] = ost_vec[0]/10;
open_mm_vec[1] = ost_vec[1]/10;
@@ -446,14 +384,12 @@ void Simulation::UpdatePositions(bool enforce_periodic_box){
if(top_->GetNumParticles() != static_cast<uint>(ent_.GetAtomCount())){
throw ost::Error("Num particles in topology and num atoms in entity are not consistent!");
}
- geom::Vec3List positions = this->GetPositions(enforce_periodic_box, true);
+ geom::Vec3List positions;
+ StateExtractor::ExtractPositions(context_, positions, enforce_periodic_box,
+ true);
ost::mol::XCSEditor ed = ent_.EditXCS(ost::mol::BUFFERED_EDIT);
- ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
- ost::mol::AtomHandleList::iterator a = atom_list.begin();
- ost::mol::AtomHandleList::iterator ae = atom_list.end();
- geom::Vec3List::iterator v = positions.begin();
- for(; a != ae; ++a, ++v){
- ed.SetAtomPos(*a,*v);
+ for(uint i = 0; i < atom_list_.size(); ++i) {
+ ed.SetAtomPos(atom_list_[i], positions[i]);
}
}
diff --git a/modules/mol/mm/src/simulation.hh b/modules/mol/mm/src/simulation.hh
index 89c58b10eb96f8b2b0a5995986d2ea80a47a43b7..7cbf06d1aededa101f55658b14f50f42d73dd6ca 100644
--- a/modules/mol/mm/src/simulation.hh
+++ b/modules/mol/mm/src/simulation.hh
@@ -143,6 +143,7 @@ private:
Simulation() { } //hidden constructor...
void Init(const ost::mol::mm::TopologyPtr top,
+ const ost::mol::EntityHandle& ent,
const SettingsPtr settings);
int TimeToNextNotification();
@@ -160,6 +161,7 @@ private:
std::vector<int> time_to_notify_;
std::map<FuncType,uint> system_force_mapper_;
ost::mol::EntityHandle ent_;
+ ost::mol::AtomHandleList atom_list_;
};
}}} //ns