diff --git a/modules/geom/pymod/export_vec3.cc b/modules/geom/pymod/export_vec3.cc
index 6ac28a6147e726f25d229da3ad100472ad349600..9a9b8f47c0850448a80a8033eb67c5df007b86e2 100644
--- a/modules/geom/pymod/export_vec3.cc
+++ b/modules/geom/pymod/export_vec3.cc
@@ -92,5 +92,7 @@ void export_Vec3()
.add_property("center", &Vec3List::GetCenter)
.add_property("inertia", &Vec3List::GetInertia)
.add_property("principal_axes", &Vec3List::GetPrincipalAxes)
+ .def("GetODRLine", &Vec3List::GetODRLine)
+ .def("FitCylinder", &Vec3List::FitCylinder)
;
}
diff --git a/modules/geom/src/composite3_op.cc b/modules/geom/src/composite3_op.cc
index 09fb2c72e7430c2fc938534ff4167052485a4b24..672896f0d9aaec18fad77a24aa8065c7c43e34e4 100644
--- a/modules/geom/src/composite3_op.cc
+++ b/modules/geom/src/composite3_op.cc
@@ -178,5 +178,50 @@ bool IsInSphere(const Sphere& s, const Vec3& v){
return Length(s.GetOrigin()-v)<=s.GetRadius();
}
+Line3 Vec3List::FitCylinder(){
+ Line3 axis=this->GetODRLine(),axis_old;
+ Real radius,res_sum_old,res_sum,delta=0.01,prec=0.001,err;
+ unsigned long n_step=1000, n_res=this->size();
+ Vec3 v,gradient;
+ radius=0.0;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ radius+=geom::Distance(axis,(*i));
+ }
+ radius/=Real(n_res);
+ res_sum=0.0;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ res_sum+=pow(Distance(axis,(*i))-radius,2.);
+ }
+ unsigned long k=0;
+ err=2.0*prec;
+ while (err>prec and k<n_step) {
+ res_sum_old=res_sum;
+ axis_old=axis;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ radius+=Distance(axis,(*i));
+ }
+ radius/=Real(n_res);
+ for (int j=0; j!=3; ++j){
+ res_sum=0.0;
+ v=Vec3(0.0,0.0,0.0);
+ v[j]=delta;
+ axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ res_sum+=pow(Distance(axis,(*i))-radius,2.);
+ }
+ gradient[j]=res_sum-res_sum_old;
+ }
+ gradient=Normalize(gradient);
+ axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
+ res_sum=0.0;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ res_sum+=pow(Distance(axis,(*i))-radius,2.);
+ }
+ err=fabs((res_sum-res_sum_old)/float(n_res));
+ k++;
+ }
+ return axis;
+}
+
} // ns
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index 84467ced5a4ddf66dd5ae23e1de054c7cb8e175d..bcc62d42ad0bb0814581f96b9750730e4a84ce7c 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -68,5 +68,12 @@ Vec3 Vec3List::GetCenter() const
return center/=this->size();
}
+Line3 Vec3List::GetODRLine()
+{
+ Vec3 center=this->GetCenter();
+ Vec3 direction=this->GetPrincipalAxes().GetRow(2);
+ return Line3(center,center+direction);
+}
+
}
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index db2b8bb7097dc76d0b74c7906a0ea4a63bda5929..2a27d500a8a0c5ca9c66127de47f0cf37f48f932 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -34,7 +34,7 @@ namespace geom {
// fw decl
class Vec2;
class Vec4;
-
+class Line3;
/// \brief Three dimensional vector class, using Real precision.
class DLLEXPORT_OST_GEOM Vec3:
@@ -193,6 +193,7 @@ inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
#include <ost/geom/vec2.hh>
#include <ost/geom/vec4.hh>
#include <ost/geom/mat3.hh>
+#include <ost/geom/composite3.hh>
namespace geom {
@@ -216,6 +217,9 @@ public:
Vec3 GetCenter() const;
Mat3 GetPrincipalAxes() const;
+ Line3 GetODRLine();
+ Line3 FitCylinder();
+
};
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
index 31377016a874dc6dcc42cc702d11a3a4058f0934..d53559215907bd13049e0afe2ab1bed29a9470c3 100644
--- a/modules/mol/alg/pymod/export_trajectory_analysis.cc
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -46,4 +46,6 @@ void export_TrajectoryAnalysis()
def("AnalyzeMinDistance", &AnalyzeMinDistance, (arg("traj"), arg("view1"), arg("view2"), arg("stride")=1));
def("AnalyzeMinDistanceBetwCenterOfMassAndView", &AnalyzeMinDistanceBetwCenterOfMassAndView, (arg("traj"), arg("view_cm"), arg("view_atoms"), arg("stride")=1));
def("AnalyzeAromaticRingInteraction", &AnalyzeAromaticRingInteraction, (arg("traj"), arg("view_ring1"), arg("view_ring2"), arg("stride")=1));
+ def("AnalyzeAlphaHelixAxis", &AnalyzeAlphaHelixAxis, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
+ def("AnalyzeBestFitLine", &AnalyzeBestFitLine, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index 2c1225e753ba76adacfa3ccfd8dc6a2ccce18a95..89b46a525c23c1776e157672390eb977337536b3 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -230,4 +230,51 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
}
return dist;
}
+
+ void AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
+ geom::Vec3List& centers, unsigned int stride)
+ {
+ CheckHandleValidity(traj);
+ if (prot_seg.GetAtomCount()==0){
+ throw std::runtime_error("EntityView is empty");
+ }
+ std::vector<unsigned long> indices_ca;
+ geom::Line3 axis;
+ directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ GetCaIndices(prot_seg, indices_ca);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ axis=frame->FitCylinder(indices_ca);
+ directions.push_back(axis.GetDirection());
+ centers.push_back(axis.GetOrigin());
+ }
+ return;
+ }
+
+ //std::vector<geom::Line3> AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg,
+ // unsigned int stride)
+ void AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
+ geom::Vec3List& centers, unsigned int stride)
+ {
+ CheckHandleValidity(traj);
+ if (prot_seg.GetAtomCount()==0){
+ throw std::runtime_error("EntityView is empty");
+ }
+ std::vector<unsigned long> indices_ca;
+ geom::Line3 axis;
+ directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ GetCaIndices(prot_seg, indices_ca);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ axis=frame->GetODRLine(indices_ca);
+ directions.push_back(axis.GetDirection());
+ centers.push_back(axis.GetOrigin());
+ }
+ return;
+ }
+
+
+
}}} //ns
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
index faa6c8ea06d657afe169ac688c8cd8e6b8148b11..6837e9fb3171a353167bb20b02210bea43a57ec6 100644
--- a/modules/mol/alg/src/trajectory_analysis.hh
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -43,5 +43,7 @@ namespace ost { namespace mol { namespace alg {
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistance(const CoordGroupHandle& traj, const EntityView& view1, const EntityView& view2,unsigned int stride=1);
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistanceBetwCenterOfMassAndView(const CoordGroupHandle& traj, const EntityView& view_cm, const EntityView& view_atoms,unsigned int stride=1);
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1, const EntityView& view_ring2,unsigned int stride=1);
+ void DLLEXPORT_OST_MOL_ALG AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
+ void DLLEXPORT_OST_MOL_ALG AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
}}}//ns
#endif
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index fdb41fce60b8032ecb6f527256b4b9c90d7cc5e1..5ec6f157b0212409cdc8d78b435208e20b814a9c 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/modules/mol/base/pymod/export_coord_frame.cc
@@ -46,6 +46,7 @@ void export_CoordFrame()
.def("GetDistanceBetwCenterOfMass", get_dist_cm)
.def("GetRMSD",get_rmsd)
.def("GetMinDistance",get_min_dist)
+ .def("GetODRLine",&geom::Vec3List::GetODRLine)
;
def("CreateCoordFrame",CreateCoordFrame);
}
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index 80ef4f58f3ea3aed2b1ae691b6cc35c70107faf0..01db3f0a2990d7cbd866e6aeb3c0922d7b91303c 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -84,6 +84,22 @@ namespace ost { namespace mol {
indices.push_back(i->GetIndex());
}
}
+
+ void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca){
+ mol::AtomView ca;
+ ResidueViewList residues=segment.GetResidueList();
+ indices_ca.reserve(residues.size());
+ //for (ResidueViewList::const_iterator res=residues.begin()+1,
+ // e=residues.end(); res!=e-1; ++res){
+ // if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
+ //}
+ for (ResidueViewList::const_iterator res=residues.begin(),
+ e=residues.end(); res!=e; ++res) {
+ ca=res->FindAtom("CA");
+ CheckHandleValidity(ca);
+ indices_ca.push_back(ca.GetIndex());
+ }
+ }
void GetMasses(const EntityView& sele, std::vector<Real>& masses){
// This function returns a vector containing the atom masses of the atoms in an EntityView
@@ -218,5 +234,29 @@ namespace ost { namespace mol {
GetIndicesAndMasses(view_cm,indices_cm,masses_cm);
return this->GetMinDistBetwCenterOfMassAndView(indices_cm,masses_cm,indices_atoms);
}
+
+ geom::Line3 CoordFrame::GetODRLine(std::vector<unsigned long>& indices_ca){
+ //Returns the best fit line to atoms with indices in indices_ca
+ geom::Vec3List atoms_pos_list;
+ atoms_pos_list.reserve(indices_ca.size());
+ for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
+ atoms_pos_list.push_back((*this)[*i1]);
+ }
+ return atoms_pos_list.GetODRLine();
+ }
+
+ geom::Line3 CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca){
+ //Returns a line which is the axis of a fitted cylinder to the atoms with indices given in indices_ca
+ geom::Vec3List atoms_pos_list;
+ atoms_pos_list.reserve(indices_ca.size());
+ for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
+ atoms_pos_list.push_back((*this)[*i1]);
+ }
+ return atoms_pos_list.FitCylinder();
+ }
+
+
+
+
}}//ns
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index 8b34cc52acb31086b8ee02d074fe0aac526e22f8..fe94ecc9f67bd4995927239aff5fd87656bcf82c 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -60,13 +60,16 @@ public:
Real GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm, std::vector<Real>& masses_cm,
std::vector<unsigned long>& indices_atoms);
Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm, const mol::EntityView& view_atoms);
+ geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca);
+ geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
};
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 GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca);
+
typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
typedef std::vector<CoordFramePtr> CoordFrameList;