From 956808eb67951964fd7954525cafc6c082b984fd Mon Sep 17 00:00:00 2001
From: nij2003 <nij2003@mac140138.med.cornell.edu>
Date: Fri, 13 May 2011 16:06:03 -0400
Subject: [PATCH] Trajectory analysis tools added
---
modules/mol/alg/doc/molalg.rst | 114 ++++++++++++++++++
.../alg/pymod/export_trajectory_analysis.cc | 17 +--
modules/mol/alg/src/trajectory_analysis.cc | 63 +++++-----
modules/mol/alg/src/trajectory_analysis.hh | 14 +--
modules/mol/base/pymod/CMakeLists.txt | 1 +
modules/mol/base/pymod/wrap_mol.cc | 2 +
modules/mol/base/src/coord_frame.cc | 104 ++++++++++------
modules/mol/base/src/coord_frame.hh | 36 +++---
modules/mol/base/src/coord_group.cc | 5 +-
modules/mol/base/src/coord_group.hh | 10 +-
10 files changed, 256 insertions(+), 110 deletions(-)
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 578d9d1d1..f061e29af 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -82,3 +82,117 @@ The following function detects steric clashes in atomic structures. Two atoms ar
:returns: The filtered :class:`~ost.mol.EntityView`
+Trajectory Analyses
+--------------------------------------------------------------------------------
+
+This is a set of funcitons used for basic trajectory analysis such as extracting positions,
+distances, angles and RMSDs. The organization is such that most functions have their counterpart
+at the individual frame level (CoordFrame) so that they can alsobe called on one frame instead
+of the whole trajectory.
+
+All these functions have a "stride" argument that defaults to stride=1, which is used to skip frames in the anlysis.
+
+
+.. function:: AnalyzeAtomPos(traj, atom1, stride=1)
+
+ This function extracts the position of an atom from a trajectory. It returns
+ a vector containing the position of the atom for each analyzed frame.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeCenterOfMassPos(traj, sele, stride=1)
+
+ This function extracts the position of the center-of-mass of a selection
+ (:class:`~ost.mol.EntityView`) from a trajectory and returns it as a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param sele: The selection from which the center of mass is computed
+ :type sele: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeDistanceBetwAtoms(traj, atom1, atom2, stride=1)
+
+ This function extracts the distance between two atoms from a trajectory
+ and returns it as a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeAngle(traj, atom1, atom2, atom3, stride=1)
+
+ This function extracts the angle between three atoms from a trajectory
+ and returns it as a vector. The second atom is taken as being the central
+ atom, so that the angle is between the vectors (atom1.pos-atom2.pos) and
+ (atom3.pos-atom2.pos).
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param atom3: The third :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeDihedralAngle(traj, atom1, atom2, atom3, atom4, stride=1)
+
+ This function extracts the dihedral angle between four atoms from a trajectory
+ and returns it as a vector. The angle is between the planes containing the first
+ three and the last three atoms.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param atom3: The third :class:`~ost.mol.AtomHandle`.
+ :param atom4: The fourth :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+.. function:: AnalyzeDistanceBetwCenterOfMass(traj, sele1, sele2, stride=1)
+
+ This function extracts the distance between the center-of-mass of two selections
+ (:class:`~ost.mol.EntityView`) from a trajectory and returns it as a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param sele1: The selection from which the first center of mass is computed
+ :type sele1: :class:`~ost.mol.EntityView`.
+ :param sele2: The selection from which the second center of mass is computed
+ :type sele2: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeRMSD(traj, reference_view, sele_view)
+
+ This function extracts the rmsd between two :class:`~ost.mol.EntityView` and returns it as a vector
+ The views don't have to be from the same entity. The reference positions are taken directly
+ from the reference_view, evaluated only once. The positions from the sele_view are evaluated for
+ each frame.
+ 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.AnalyzeRMSD(t,sele,sele)
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param reference_view: The selection used as reference structure
+ :type reference_view: :class:`~ost.mol.EntityView`.
+ :param sele_view: The selection compared to the reference_view
+ :type sele_view: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+
+
+
+
+
+
+
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
index ff55f92ec..61e153bb3 100644
--- a/modules/mol/alg/pymod/export_trajectory_analysis.cc
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -25,6 +25,7 @@ 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)
@@ -32,14 +33,14 @@ 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());
+ def("AnalyzeAtomPos",&AnalyzeAtomPos, (arg("traj"), arg("atom"), arg("stride")=1));
+ def("AnalyzeCenterOfMassPos",&AnalyzeCenterOfMassPos, (arg("traj"), arg("selection"), arg("stride")=1));
+ def("AnalyzeDistanceBetwAtoms",&AnalyzeDistanceBetwAtoms, (arg("traj"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeAngle",&AnalyzeAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeDihedralAngle",&AnalyzeDihedralAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeDistanceBetwCenterOfMass",&AnalyzeDistanceBetwCenterOfMass, (arg("traj"), arg("selection"), arg("selection"), arg("stride")=1));
+ def("AnalyzeRMSD",&AnalyzeRMSD, (arg("traj"), arg("reference_view"), arg("selection"), arg("stride")=1));
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index 5f9e77623..c8660fc99 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -21,10 +21,8 @@
* 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>
@@ -34,72 +32,71 @@
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
+geom::Vec3List AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
+//This function extracts the position of an atom from a trajectory and returns it as a vector of geom::Vec3
{
- traj.CheckValidity();
+ CheckHandleValidity(traj);
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));
+ pos.push_back(frame->GetAtomPos(i1));
}
return pos;
}
-std::vector<Real> ExtractDistance(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+std::vector<Real> AnalyzeDistanceBetwAtoms(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();
+ CheckHandleValidity(traj);
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));
+ dist.push_back(frame->GetDistanceBetwAtoms(i1,i2));
}
return dist;
}
-std::vector<Real> ExtractAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+std::vector<Real> AnalyzeAngle(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();
+ CheckHandleValidity(traj);
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));
+ ang.push_back(frame->GetAngle(i1,i2,i3));
}
return ang;
}
-std::vector<Real> ExtractDihedral(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+std::vector<Real> AnalyzeDihedralAngle(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();
+ CheckHandleValidity(traj);
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));
+ 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
+geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride)
+//This function extracts the position of the CenterOfMass of a selection (entity view) from a trajectory
//and returns it as a vector.
{
- traj.CheckValidity();
+ CheckHandleValidity(traj);
geom::Vec3List pos;
pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices;
@@ -107,17 +104,17 @@ geom::Vec3List ExtractCMPosition(const CoordGroupHandle& traj, const EntityView&
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));
+ pos.push_back(frame->GetCenterOfMassPos(indices,masses));
}
return pos;
}
-std::vector<Real> ExtractCMDistance(const CoordGroupHandle& traj, const EntityView& Sele1,
+std::vector<Real> AnalyzeDistanceBetwCenterOfMass(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
+//This function extracts the distance between the CenterOfMass of two selection (entity views) from a trajectory
//and returns it as a vector.
{
- traj.CheckValidity();
+ CheckHandleValidity(traj);
std::vector<Real> dist;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices1,indices2;
@@ -126,21 +123,21 @@ std::vector<Real> ExtractCMDistance(const CoordGroupHandle& traj, const EntityVi
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));
+ dist.push_back(frame->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2));
}
return dist;
}
-std::vector<Real> ExtractRMSD(const CoordGroupHandle& traj, const EntityView& ReferenceView,
- const EntityView& SeleView, unsigned int stride)
+std::vector<Real> AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view,
+ const EntityView& sele_view, 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)
+// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,Sele,Sele)
{
- traj.CheckValidity();
- int count_ref=ReferenceView.GetAtomCount();
- int count_sele=SeleView.GetAtomCount();
+ CheckHandleValidity(traj);
+ 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");
}
@@ -148,11 +145,11 @@ std::vector<Real> ExtractRMSD(const CoordGroupHandle& traj, const EntityView& Re
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);
+ GetIndices(reference_view, sele_indices);
+ GetPositions(sele_view, 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));
+ rmsd.push_back(frame->GetRMSD(ref_pos,sele_indices));
}
return rmsd;
}
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
index d001d751f..d2947e63b 100644
--- a/modules/mol/alg/src/trajectory_analysis.hh
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -33,13 +33,13 @@
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);
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& sele,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& sele1, const EntityView& sele2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDihedralAngle(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 AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view, const EntityView& sele,unsigned int stride=1);
}}}//ns
diff --git a/modules/mol/base/pymod/CMakeLists.txt b/modules/mol/base/pymod/CMakeLists.txt
index 29c336dca..87369f260 100644
--- a/modules/mol/base/pymod/CMakeLists.txt
+++ b/modules/mol/base/pymod/CMakeLists.txt
@@ -5,6 +5,7 @@ export_bond.cc
export_chain.cc
export_chain_view.cc
export_coord_group.cc
+export_coord_frame.cc
export_editors.cc
export_entity.cc
export_entity_view.cc
diff --git a/modules/mol/base/pymod/wrap_mol.cc b/modules/mol/base/pymod/wrap_mol.cc
index 0ab6b4568..f29d1377f 100644
--- a/modules/mol/base/pymod/wrap_mol.cc
+++ b/modules/mol/base/pymod/wrap_mol.cc
@@ -39,6 +39,7 @@ void export_AtomView();
void export_ResidueView();
void export_Editors();
void export_CoordGroup();
+void export_CoordFrame();
void export_PropertyID();
void export_BoundingBox();
void export_QueryViewWrapper();
@@ -65,6 +66,7 @@ BOOST_PYTHON_MODULE(_ost_mol)
export_EntityView();
export_Editors();
export_CoordGroup();
+ export_CoordFrame();
export_PropertyID();
export_BoundingBox();
export_QueryViewWrapper();
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index 5a81953c7..80105c128 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -31,36 +31,53 @@
namespace ost { namespace mol {
- geom::Vec3 CoordFrame::GetAtomPosition(const AtomHandle& atom){
- return this->GetAtomPosition(atom.GetIndex());
+ CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos)
+ {
+ return CoordFrame(atom_pos);
}
- geom::Vec3 CoordFrame::GetAtomPosition(const int& index){
- return (*this)[index];
+
+ geom::Vec3 CoordFrame::GetAtomPos(const AtomHandle& atom){
+ //Returns the position of the atom in this frame
+ return this->GetAtomPos(atom.GetIndex());
+ }
+ geom::Vec3 CoordFrame::GetAtomPos(int i1){
+ //Returns the position in this frame of the atom with index i1
+ return (*this)[i1];
}
- Real CoordFrame::GetDistance(const AtomHandle& a1, const AtomHandle& a2){
- return this->GetDistance(a1.GetIndex(),a2.GetIndex());
+ Real CoordFrame::GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2){
+ //Return the distance in this frame between two atoms
+ return this->GetDistanceBetwAtoms(a1.GetIndex(),a2.GetIndex());
}
- Real CoordFrame::GetDistance(const int& i1, const int& i2){
+ Real CoordFrame::GetDistanceBetwAtoms(int i1, int i2){
+ //Return the distance in this frame between the two atoms with indices i1 and i2
return geom::Distance((*this)[i1],(*this)[i2]);
}
Real CoordFrame::GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3){
+ //Returns the Angle between the three atoms, a2 being considered as the central atom
+ //i.e the angle between the vector (a1.pos-a2.pos) and (a3.pos-a2.pos)
return this->GetAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex());
}
- Real CoordFrame::GetAngle(const int& i1, const int& i2, const int& i3){
+ Real CoordFrame::GetAngle(int i1, int i2, int i3){
+ //Returns the angl between the three atoms with indices i1,i2,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){
+ Real CoordFrame::GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2,
+ const AtomHandle& a3, const AtomHandle& a4){
+ //Returns the Dihedral angle between the four atoms a1,a2,a3,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){
+ Real CoordFrame::GetDihedralAngle(int i1, int i2, int i3, int i4){
+ //Returns the Dihedral angle between the four atoms with indices i1,i2,i3,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();
+ void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices){
+ // This function returns a vector containing the atom indices of the atoms in an EntityView
+ // It is used to accelerate the extraction of information from a trajectory
+ AtomViewList atoms=sele.GetAtomList();
indices.reserve(atoms.size());
for (AtomViewList::const_iterator i=atoms.begin(),
e=atoms.end(); i!=e; ++i) {
@@ -68,9 +85,11 @@ namespace ost { namespace mol {
}
}
- void GetMasses(const EntityView& Sele, std::vector<Real>& masses){
+ void GetMasses(const EntityView& sele, std::vector<Real>& masses){
+ // This function returns a vector containing the atom masses of the atoms in an EntityView
+ // It is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
Real mass_tot=0.0;
- AtomViewList atoms=Sele.GetAtomList();
+ AtomViewList atoms=sele.GetAtomList();
masses.reserve(atoms.size());
for (AtomViewList::const_iterator i=atoms.begin(),
e=atoms.end(); i!=e; ++i) {
@@ -84,26 +103,30 @@ namespace ost { namespace mol {
}
- void GetIndicesAndMasses(const EntityView& Sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
- GetIndices(Sele, indices);
- GetMasses(Sele, masses);
+ 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) {
+ void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos){
+ //Returns the positions of all the atoms in the EntityView
+ 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){
+ geom::Vec3 CoordFrame::GetCenterOfMassPos(const EntityView& sele){
+ //Returns the position of the centor of mass of the atoms in the EntityView
std::vector<unsigned long> indices;
std::vector<Real> masses;
- GetIndicesAndMasses(Sele,indices,masses);
- return this->GetCMPosition(indices,masses);
+ GetIndicesAndMasses(sele,indices,masses);
+ return this->GetCenterOfMassPos(indices,masses);
}
- geom::Vec3 CoordFrame::GetCMPosition(std::vector<unsigned long>& indices,std::vector<Real>& masses){
+ geom::Vec3 CoordFrame::GetCenterOfMassPos(std::vector<unsigned long>& indices,std::vector<Real>& masses){
+ //Returns the position of the centor of mass of the atoms from which the indices and masses are passed
+ //as vectors in argument
geom::Vec3 v;
for (unsigned int i=0,e=indices.size();i!=e; i++) {
v+=masses[i]*(*this)[indices[i]];
@@ -111,22 +134,27 @@ namespace ost { namespace mol {
return v;
}
- Real CoordFrame::GetCMDistance(const EntityView& Sele1,const EntityView& Sele2){
+ Real CoordFrame::GetDistanceBetwCenterOfMass(const EntityView& sele1,const EntityView& sele2){
+ //Returns the distance between the centers of mass of the two EntityViews
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);
+ GetIndicesAndMasses(sele1,indices1,masses1);
+ GetIndicesAndMasses(sele2,indices2,masses2);
+ return this->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2);
}
- Real CoordFrame::GetCMDistance(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,
+ Real CoordFrame::GetDistanceBetwCenterOfMass(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);
+ //Returns the distance between the centers of mass of the two groups of atoms from which the
+ //indices and masses are given as vectors as argument
+ geom::Vec3 v1=this->GetCenterOfMassPos(indices1, masses1);
+ geom::Vec3 v2=this->GetCenterOfMassPos(indices2, masses2);
return geom::Distance(v1,v2);
}
Real CoordFrame::GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele){
+ //Returns the RMSD between the positions of the atoms whose indices are given in indices_sele and the positions
+ //given in ref_pos
Real rmsd=0.0,val;
for (unsigned int i1=0; i1!=indices_sele.size(); ++i1) {
geom::Vec3 av_sele=(*this)[indices_sele[i1]];
@@ -134,19 +162,21 @@ namespace ost { namespace mol {
val=geom::Length2(av_ref-av_sele);
rmsd+=val;
}
- return rmsd/indices_sele.size();
+ return pow(rmsd/indices_sele.size(),0.5);
}
- Real CoordFrame::GetRMSD(const EntityView& Reference_View,const EntityView& Sele_View){
- int count_ref=Reference_View.GetAtomCount();
- int count_sele=Sele_View.GetAtomCount();
+ Real CoordFrame::GetRMSD(const EntityView& reference_view,const EntityView& sele_view){
+ //Return the rmsd between two EntityViews. The reference positions are taken directly from the reference_view
+ //whereas they are taken from this CoordFrame for the 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);
+ GetIndices(sele_view,indices_sele);
+ GetPositions(reference_view,ref_pos);
return this->GetRMSD(ref_pos,indices_sele);
}
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index a668c098d..c856e583c 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -40,30 +40,34 @@ public:
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);
+ geom::Vec3 GetAtomPos(const AtomHandle& atom);
+ geom::Vec3 GetAtomPos(int atom_index);
+ Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2);
+ Real GetDistanceBetwAtoms(int atom1_index, 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 GetAngle(int atom1_index, int atom2_index, 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);
+ Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index);
+ geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele);
+ geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices, std::vector<Real>& masses);
+ Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2);
+ Real GetDistanceBetwCenterOfMass(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);
+ 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;
+// factory method
+// create a frame froma Vec3List containing the positions of the atoms
+ DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos);
}}
diff --git a/modules/mol/base/src/coord_group.cc b/modules/mol/base/src/coord_group.cc
index eaf7bff4c..f92c58e22 100644
--- a/modules/mol/base/src/coord_group.cc
+++ b/modules/mol/base/src/coord_group.cc
@@ -84,8 +84,9 @@ CoordGroupHandle::operator bool() const
return this->IsValid();
}
-void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
-{
+//void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
+ void CoordGroupHandle::AddFrame(const geom::Vec3List& clist)
+ {
this->CheckValidity();
if (source_->IsMutable()) {
source_->AddFrame(clist);
diff --git a/modules/mol/base/src/coord_group.hh b/modules/mol/base/src/coord_group.hh
index b64d3ba2f..2b6d6fe7f 100644
--- a/modules/mol/base/src/coord_group.hh
+++ b/modules/mol/base/src/coord_group.hh
@@ -25,7 +25,6 @@
#include <vector>
#include <boost/shared_array.hpp>
-//#include <ost/mol/alg/trajectory_analysis.hh>
#include "atom_handle.hh"
#include "coord_source.hh"
@@ -65,7 +64,8 @@ public:
void Capture(uint frame);
/// \brief add frame
- void AddFrame(const std::vector<geom::Vec3>& clist);
+ //void AddFrame(const std::vector<geom::Vec3>& clist);
+ void AddFrame(const geom::Vec3List& clist);
void AddFrames(const CoordGroupHandle& cg);
/// \brief set an indidivial atom position in the given frame
@@ -89,13 +89,9 @@ public:
CoordGroupHandle(CoordSourcePtr source);
- /*
- //friend geom::Vec3List mol::alg::ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
- */
-//private:
- void CheckValidity() const;
private:
+ void CheckValidity() const;
CoordSourcePtr source_;
};
--
GitLab