diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index ec9a6cb75d2b269bc229ef0a2311bdc9fe15b6da..db2b8bb7097dc76d0b74c7906a0ea4a63bda5929 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -205,7 +205,7 @@ public:
Vec3List(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
Vec3List(const Vec3List& rhs) : base_type(rhs) { }
-
+ Vec3List(const base_type& rhs) : base_type(rhs) { }
Vec3List& operator=(const Vec3List& rhs)
{
base_type::operator=(rhs);
diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index 2fbf108cbe9133fcab98cd55b30fb99c74dabd2a..0ffb3f65183cb217c5fcfe6b1093a884ea42bca9 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -2,6 +2,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
wrap_mol_alg.cc
export_svd_superpose.cc
export_clash.cc
+ export_trajectory_analysis.cc
)
set(OST_MOL_ALG_PYMOD_MODULES
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ff55f92ec1369a4e3890e476bb67dcfb36e3b38e
--- /dev/null
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -0,0 +1,45 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 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>
+using namespace boost::python;
+
+#include <ost/mol/alg/trajectory_analysis.hh>
+
+using namespace ost;
+using namespace ost::mol::alg;
+
+//"thin wrappers" for default parameters
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractapos, ExtractAtomPosition, 2, 3)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmpos, ExtractCMPosition, 2, 3)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractdist, ExtractDistance, 3, 4)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractang, ExtractAngle, 4, 5)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractdih, ExtractDihedral, 5, 6)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmdist, ExtractCMDistance, 3, 4)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractrmsd, ExtractRMSD, 3, 4)
+
+void export_TrajectoryAnalysis()
+{
+ def("ExtractAtomPosition",ExtractAtomPosition,extractapos());
+ def("ExtractCMPosition",&ExtractCMPosition,extractcmpos());
+ def("ExtractDistance",&ExtractDistance,extractdist());
+ def("ExtractAngle",&ExtractAngle,extractang());
+ def("ExtractDihedral",&ExtractDihedral,extractdih());
+ def("ExtractCMDistance",&ExtractCMDistance,extractcmdist());
+ def("ExtractRMSD",&ExtractRMSD,extractrmsd());
+}
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 82803c983ac9eed3ecc0364326310476be57ece9..48e19d84fc7b35e599454f809c07b27abcf1c461 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -26,7 +26,8 @@ using namespace boost::python;
using namespace ost;
void export_svdSuperPose();
-
+void export_TrajectoryAnalysis();
+void export_Clash();
#if OST_IMG_ENABLED
void export_entity_to_density();
#endif
@@ -42,6 +43,7 @@ mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterC
BOOST_PYTHON_MODULE(_ost_mol_alg)
{
export_svdSuperPose();
+ export_TrajectoryAnalysis();
#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 ae27e33e6f0160d5556d0028cc3d75fe32e16c62..1bd52057e93d59a5de0ec09edab397b5bc812261 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -7,6 +7,7 @@ set(OST_MOL_ALG_HEADERS
filter_clashes.hh
construct_cbeta.hh
clash_score.hh
+ trajectory_analysis.hh
)
set(OST_MOL_ALG_SOURCES
@@ -17,6 +18,7 @@ set(OST_MOL_ALG_SOURCES
superpose_frames.cc
filter_clashes.cc
construct_cbeta.cc
+ trajectory_analysis.cc
)
set(MOL_ALG_DEPS ost_mol ost_seq)
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5f9e7762356a1120819f423e6f59bfe2126dcff1
--- /dev/null
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -0,0 +1,160 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 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
+//------------------------------------------------------------------------------
+
+/*
+ * Author Niklaus Johner
+ */
+#include <stdexcept>
+//#include <boost/bind.hpp>
+#include <ost/base.hh>
+#include <ost/geom/vec3.hh>
+//#include <ost/mol/alg/svd_superpose.hh>
+#include <ost/base.hh>
+#include <ost/geom/geom.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/coord_group.hh>
+#include "trajectory_analysis.hh"
+
+
+namespace ost { namespace mol { namespace alg {
+
+geom::Vec3List ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
+//This function extracts the position of an atom returns it as a vector of geom::Vec3
+//Doesn't work in python, because it cannot create the vector of geom::Vec3
+{
+ traj.CheckValidity();
+ geom::Vec3List pos;
+ pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ pos.push_back(frame->GetAtomPosition(i1));
+ }
+ return pos;
+}
+
+std::vector<Real> ExtractDistance(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ unsigned int stride)
+//This function extracts the distance between two atoms from a trajectory and returns it as a vector
+{
+ traj.CheckValidity();
+ std::vector<Real> dist;
+ dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex();
+ int i2=a2.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ dist.push_back((*frame).GetDistance(i1,i2));
+ }
+ return dist;
+}
+
+std::vector<Real> ExtractAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ const AtomHandle& a3, unsigned int stride)
+//This function extracts the angle between three atoms from a trajectory and returns it as a vector
+{
+ traj.CheckValidity();
+ std::vector<Real> ang;
+ ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ ang.push_back((*frame).GetAngle(i1,i2,i3));
+ }
+ return ang;
+}
+
+std::vector<Real> ExtractDihedral(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ const AtomHandle& a3, const AtomHandle& a4, unsigned int stride)
+//This function extracts the dihedral angle between four atoms from a trajectory and returns it as a vector
+{
+ traj.CheckValidity();
+ std::vector<Real> ang;
+ ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex(),i4=a4.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ ang.push_back((*frame).GetDihedralAngle(i1,i2,i3,i4));
+ }
+ return ang;
+}
+
+geom::Vec3List ExtractCMPosition(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride)
+//This function extracts the position of the CM of two selection (entity views) from a trajectory
+//and returns it as a vector.
+ {
+ traj.CheckValidity();
+ geom::Vec3List pos;
+ pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> indices;
+ std::vector<Real> masses;
+ GetIndicesAndMasses(Sele, indices, masses);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ pos.push_back(frame->GetCMPosition(indices,masses));
+ }
+ return pos;
+}
+
+std::vector<Real> ExtractCMDistance(const CoordGroupHandle& traj, const EntityView& Sele1,
+ const EntityView& Sele2, unsigned int stride)
+//This function extracts the distance between the CM of two selection (entity views) from a trajectory
+//and returns it as a vector.
+ {
+ traj.CheckValidity();
+ std::vector<Real> dist;
+ dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> indices1,indices2;
+ std::vector<Real> masses1,masses2;
+ GetIndicesAndMasses(Sele1, indices1, masses1);
+ GetIndicesAndMasses(Sele2, indices2, masses2);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ dist.push_back((*frame).GetCMDistance(indices1,masses1,indices2,masses2));
+ }
+ return dist;
+}
+
+std::vector<Real> ExtractRMSD(const CoordGroupHandle& traj, const EntityView& ReferenceView,
+ const EntityView& SeleView, unsigned int stride)
+// This function extracts the rmsd between two entity views and returns it as a vector
+// The views don't have to be from the same entity
+// If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
+// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.ExtractRMSD(t,Sele,Sele)
+ {
+ traj.CheckValidity();
+ int count_ref=ReferenceView.GetAtomCount();
+ int count_sele=SeleView.GetAtomCount();
+ if (count_ref!=count_sele){
+ throw std::runtime_error("atom counts of the two views are not equal");
+ }
+ std::vector<Real> rmsd;
+ rmsd.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> sele_indices;
+ std::vector<geom::Vec3> ref_pos;
+ GetIndices(ReferenceView, sele_indices);
+ GetPositions(SeleView, ref_pos);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ rmsd.push_back((*frame).GetRMSD(ref_pos,sele_indices));
+ }
+ return rmsd;
+}
+
+}}} //ns
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d001d751f0da5deaeae209792973ae5adbff9bb8
--- /dev/null
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -0,0 +1,46 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 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
+//------------------------------------------------------------------------------
+
+/*
+ * Niklaus Johner
+ */
+#ifndef OST_TRAJECTORY_ANALYSIS_HH
+#define OST_TRAJECTORY_ANALYSIS_HH
+
+#include <ost/mol/alg/module_config.hh>
+
+#include <ost/base.hh>
+#include <ost/geom/geom.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/coord_group.hh>
+
+
+namespace ost { namespace mol { namespace alg {
+
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG ExtractCMPosition(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractDistance(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractCMDistance(const CoordGroupHandle& traj, const EntityView& Sele1, const EntityView& Sele2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractDihedral(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractRMSD(const CoordGroupHandle& traj, const EntityView& Reference_View, const EntityView& Sele,unsigned int stride=1);
+
+
+}}}//ns
+#endif
diff --git a/modules/mol/base/pymod/export_coord_group.cc b/modules/mol/base/pymod/export_coord_group.cc
index ebd0f1f33daa0fc703d92f4e5351626741b52f61..c1fe4d11cf28ba2df88cf34a32af582b1da4fedc 100644
--- a/modules/mol/base/pymod/export_coord_group.cc
+++ b/modules/mol/base/pymod/export_coord_group.cc
@@ -39,7 +39,15 @@ namespace {
void (CoordGroupHandle::*capture1)()=&CoordGroupHandle::Capture;
void (CoordGroupHandle::*capture2)(uint)=&CoordGroupHandle::Capture;
}
-
+/*
+//"thin wrappers" for default parameters
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractapos, ExtractAtomPosition, 1, 2)
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractdist, ExtractDistance, 2, 3)
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractang, ExtractAngle, 3, 4)
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractdih, ExtractDihedral, 4, 5)
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractcmdist, ExtractCMDistance, 2, 3)
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(extractrmsd, ExtractRMSD, 2, 3)
+*/
void export_CoordGroup()
{
class_<CoordGroupHandle>("CoordGroupHandle",no_init)
@@ -61,6 +69,14 @@ void export_CoordGroup()
.def("__getitem__",cg_getitem)
.def("__setitem__",cg_setitem)
.def("Filter", &CoordGroupHandle::Filter)
+ /*
+ .def("ExtractAtomPosition",&CoordGroupHandle::ExtractAtomPosition,extractapos())
+ .def("ExtractDistance",&CoordGroupHandle::ExtractDistance,extractdist())
+ .def("ExtractAngle",&CoordGroupHandle::ExtractAngle,extractang())
+ .def("ExtractDihedral",&CoordGroupHandle::ExtractDihedral,extractdih())
+ .def("ExtractCMDistance",&CoordGroupHandle::ExtractCMDistance,extractcmdist())
+ .def("ExtractRMSD",&CoordGroupHandle::ExtractRMSD,extractrmsd())
+ */
;
def("CreateCoordGroup",CreateCoordGroup);
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..5a81953c71ff6eb4ed11983cc9671b8d7a4ff45b 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -0,0 +1,153 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 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
+//------------------------------------------------------------------------------
+
+/*
+ Author: Niklaus Johner
+ */
+
+#include <ost/invalid_handle.hh>
+#include <ost/integrity_error.hh>
+#include <ost/log.hh>
+#include <ost/mol/in_mem_coord_source.hh>
+#include <ost/mol/view_op.hh>
+#include <ost/mol/mol.hh>
+#include "coord_frame.hh"
+
+namespace ost { namespace mol {
+
+ geom::Vec3 CoordFrame::GetAtomPosition(const AtomHandle& atom){
+ return this->GetAtomPosition(atom.GetIndex());
+ }
+ geom::Vec3 CoordFrame::GetAtomPosition(const int& index){
+ return (*this)[index];
+ }
+
+ Real CoordFrame::GetDistance(const AtomHandle& a1, const AtomHandle& a2){
+ return this->GetDistance(a1.GetIndex(),a2.GetIndex());
+ }
+ Real CoordFrame::GetDistance(const int& i1, const int& i2){
+ return geom::Distance((*this)[i1],(*this)[i2]);
+ }
+
+ Real CoordFrame::GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3){
+ return this->GetAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex());
+ }
+ Real CoordFrame::GetAngle(const int& i1, const int& i2, const int& i3){
+ return geom::Angle((*this)[i1]-(*this)[i2],(*this)[i3]-(*this)[i2]);
+ }
+
+ Real CoordFrame::GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4){
+ return this->GetDihedralAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex(),a4.GetIndex());
+ }
+ Real CoordFrame::GetDihedralAngle(const int& i1, const int& i2, const int& i3, const int& i4){
+ return geom::DihedralAngle((*this)[i1],(*this)[i2],(*this)[i3],(*this)[i4]);
+ }
+
+ void GetIndices(const EntityView& Sele, std::vector<unsigned long>& indices){
+ AtomViewList atoms=Sele.GetAtomList();
+ indices.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ indices.push_back(i->GetIndex());
+ }
+ }
+
+ void GetMasses(const EntityView& Sele, std::vector<Real>& masses){
+ Real mass_tot=0.0;
+ AtomViewList atoms=Sele.GetAtomList();
+ masses.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ masses.push_back(i->GetMass());
+ mass_tot=mass_tot+i->GetMass();
+ }
+ for (std::vector<Real>::iterator
+ j=masses.begin(), e2=masses.end(); j!=e2; ++j) {
+ (*j)/=mass_tot;
+ }
+ }
+
+
+ void GetIndicesAndMasses(const EntityView& Sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
+ GetIndices(Sele, indices);
+ GetMasses(Sele, masses);
+ }
+
+ void GetPositions(const EntityView& Sele, std::vector<geom::Vec3>& ref_pos){
+ ref_pos.reserve(Sele.GetAtomCount());
+ for (mol::AtomViewIter a=Sele.AtomsBegin(),e=Sele.AtomsEnd(); a!=e; ++a) {
+ ref_pos.push_back((*a).GetPos());
+ }
+ }
+
+ geom::Vec3 CoordFrame::GetCMPosition(const EntityView& Sele){
+ std::vector<unsigned long> indices;
+ std::vector<Real> masses;
+ GetIndicesAndMasses(Sele,indices,masses);
+ return this->GetCMPosition(indices,masses);
+ }
+
+ geom::Vec3 CoordFrame::GetCMPosition(std::vector<unsigned long>& indices,std::vector<Real>& masses){
+ geom::Vec3 v;
+ for (unsigned int i=0,e=indices.size();i!=e; i++) {
+ v+=masses[i]*(*this)[indices[i]];
+ }
+ return v;
+ }
+
+ Real CoordFrame::GetCMDistance(const EntityView& Sele1,const EntityView& Sele2){
+ std::vector<unsigned long> indices1,indices2;
+ std::vector<Real> masses1,masses2;
+ GetIndicesAndMasses(Sele1,indices1,masses1);
+ GetIndicesAndMasses(Sele2,indices2,masses2);
+ return this->GetCMDistance(indices1,masses1,indices2,masses2);
+ }
+
+ Real CoordFrame::GetCMDistance(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,
+ std::vector<unsigned long>& indices2,std::vector<Real>& masses2){
+ geom::Vec3 v1=this->GetCMPosition(indices1, masses1);
+ geom::Vec3 v2=this->GetCMPosition(indices2, masses2);
+ return geom::Distance(v1,v2);
+ }
+
+ Real CoordFrame::GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele){
+ Real rmsd=0.0,val;
+ for (unsigned int i1=0; i1!=indices_sele.size(); ++i1) {
+ geom::Vec3 av_sele=(*this)[indices_sele[i1]];
+ geom::Vec3 av_ref=ref_pos[i1];
+ val=geom::Length2(av_ref-av_sele);
+ rmsd+=val;
+ }
+ return rmsd/indices_sele.size();
+ }
+
+ Real CoordFrame::GetRMSD(const EntityView& Reference_View,const EntityView& Sele_View){
+ int count_ref=Reference_View.GetAtomCount();
+ int count_sele=Sele_View.GetAtomCount();
+ if (count_ref!=count_sele){
+ throw std::runtime_error("atom counts of the two views are not equal");
+ }
+ std::vector<unsigned long> indices_sele;
+ std::vector<geom::Vec3> ref_pos;
+ GetIndices(Sele_View,indices_sele);
+ GetPositions(Reference_View,ref_pos);
+ return this->GetRMSD(ref_pos,indices_sele);
+ }
+
+}}//ns
\ No newline at end of file
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index c55b855a45ee8508aacece0145a04678dcd68ef0..a668c098d5de4b48be2c00479f468de713fe0e0e 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -20,19 +20,51 @@
#define OST_MOL_COORD_FRAME_HH
/*
- Author: Marco Biasini
+ Author: Marco Biasini and Niklaus Johner
*/
#include <boost/shared_ptr.hpp>
#include <ost/geom/geom.hh>
#include <ost/mol/module_config.hh>
+#include <ost/mol/entity_view.hh>
+#include "atom_handle.hh"
namespace ost { namespace mol {
-
-typedef std::vector<geom::Vec3> CoordFrame;
+class DLLEXPORT_OST_MOL CoordFrame : public geom::Vec3List{
+public:
+ typedef geom::Vec3List base_type;
+ CoordFrame() : base_type() {}
+
+ CoordFrame(size_t size, const geom::Vec3& value=geom::Vec3()) : base_type(size, value) {}
+ CoordFrame(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
+ CoordFrame(const base_type& rhs) : base_type(rhs) { }
+ CoordFrame(const std::vector<geom::Vec3>& rhs) : base_type(rhs) { }
+
+ geom::Vec3 GetAtomPosition(const AtomHandle& atom);
+ geom::Vec3 GetAtomPosition(const int& atom_index);
+ Real GetDistance(const AtomHandle& a1, const AtomHandle& a2);
+ Real GetDistance(const int& atom1_index, const int& atom2_index);
+ Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3);
+ Real GetAngle(const int& atom1_index, const int& atom2_index, const int& atom3_index);
+ Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4);
+ Real GetDihedralAngle(const int& a1_index, const int& a2_index, const int& a3_index, const int& a4_index);
+ geom::Vec3 GetCMPosition(const mol::EntityView& Sele);
+ geom::Vec3 GetCMPosition(std::vector<unsigned long>& indices,std::vector<Real>& masses);
+ Real GetCMDistance(const mol::EntityView& Sele1, const mol::EntityView& Sele2);
+ Real GetCMDistance(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,std::vector<unsigned long>& indices2,std::vector<Real>& masses2);
+ Real GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele);
+ Real GetRMSD(const mol::EntityView& Reference_View, const mol::EntityView& Sele_View);
+};
+
+ void GetIndices(const EntityView& Sele, std::vector<unsigned long>& indices);
+ void GetMasses(const EntityView& Sele,std::vector<Real>& masses);
+ void GetIndicesAndMasses(const EntityView& Sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
+ void GetPositions(const EntityView& Sele, std::vector<geom::Vec3>& ref_pos);
+
typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
typedef std::vector<CoordFramePtr> CoordFrameList;
+
}}
#endif
diff --git a/modules/mol/base/src/coord_group.cc b/modules/mol/base/src/coord_group.cc
index a68985d6ceab558e821b6d6f05dfb8eae19cc00b..eaf7bff4cbafca4c6ed3091ed728dd101f11846d 100644
--- a/modules/mol/base/src/coord_group.cc
+++ b/modules/mol/base/src/coord_group.cc
@@ -24,6 +24,7 @@
#include <ost/mol/mol.hh>
#include "coord_group.hh"
+
namespace ost { namespace mol {
CoordGroupHandle CreateCoordGroup(const AtomHandleList& atomlist)
@@ -203,5 +204,6 @@ CoordGroupHandle CoordGroupHandle::Filter(const EntityView& selected) const
}
return filtered_cg;
}
-
+
+
}} // ns
diff --git a/modules/mol/base/src/coord_group.hh b/modules/mol/base/src/coord_group.hh
index a3e582b15c12bfcc787d21a30b743cc2c04c1488..b64d3ba2fce31b2d6964fe41209ff4ca29a7c9e4 100644
--- a/modules/mol/base/src/coord_group.hh
+++ b/modules/mol/base/src/coord_group.hh
@@ -25,7 +25,7 @@
#include <vector>
#include <boost/shared_array.hpp>
-
+//#include <ost/mol/alg/trajectory_analysis.hh>
#include "atom_handle.hh"
#include "coord_source.hh"
@@ -88,9 +88,14 @@ public:
CoordGroupHandle Filter(const EntityView& selected) const;
CoordGroupHandle(CoordSourcePtr source);
-private:
- void CheckValidity() const;
+ /*
+ //friend geom::Vec3List mol::alg::ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
+ */
+
+//private:
+ void CheckValidity() const;
+private:
CoordSourcePtr source_;
};