diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index bcc62d42ad0bb0814581f96b9750730e4a84ce7c..ea07f3cca30b46f4aeef0139d7f4a0a0aab9345b 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -75,5 +75,11 @@ Line3 Vec3List::GetODRLine()
return Line3(center,center+direction);
}
-
+Plane Vec3List::GetODRPlane()
+{
+ Vec3 origin=this->GetCenter();
+ Vec3 normal=this->GetPrincipalAxes().GetRow(0);
+ return Plane(origin,normal);
+}
+
}
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index a2d1d3cea92d4fb79dfb71650eaa17fa1cadce1a..ac3458837fc20245ac7d185a56a58cbf69b8eb08 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -35,6 +35,7 @@ namespace geom {
class Vec2;
class Vec4;
class Line3;
+class Plane;
/// \brief Three dimensional vector class, using Real precision.
class DLLEXPORT_OST_GEOM Vec3:
@@ -218,6 +219,7 @@ public:
Mat3 GetPrincipalAxes() const;
Line3 GetODRLine();
+ Plane GetODRPlane();
Line3 FitCylinder(const Vec3 initial_direction, const Vec3 center);
};
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
index 5671dbca091c49746976714679c04dd48a68fd36..7e42d3eee3afb3df3c9dac7680c9216f93665635 100644
--- a/modules/mol/alg/pymod/export_trajectory_analysis.cc
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -48,6 +48,7 @@ void export_TrajectoryAnalysis()
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));
+ def("AnalyzeBestFitPlane", &AnalyzeBestFitPlane, (arg("traj"), arg("protein_segment"), arg("normals"), arg("origins"), arg("stride")=1));
def("AnalyzeHelicity", &AnalyzeHelicity, (arg("traj"), arg("protein_segment"), arg("stride")=1));
def("CreateMeanStructure", &CreateMeanStructure, (arg("traj"), arg("selection"), arg("from")=0, arg("to")=-1, arg("stride")=1 ));
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index c7222e599fc043cc57ddc94dca9495bd664c175a..c5cc63c682c2015d559cc149f320a09130cb4d90 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -272,7 +272,7 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
geom::Line3 axis;
directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
- GetCaIndices(prot_seg, indices_ca);
+ GetIndices(prot_seg, indices_ca);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
axis=frame->GetODRLine(indices_ca);
@@ -281,7 +281,28 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
}
return;
}
-
+
+ void AnalyzeBestFitPlane(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& normals,
+ geom::Vec3List& origins, unsigned int stride)
+ {
+ CheckHandleValidity(traj);
+ if (prot_seg.GetAtomCount()==0){
+ throw std::runtime_error("EntityView is empty");
+ }
+ std::vector<unsigned long> indices_ca;
+ geom::Plane best_plane;
+ normals.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ origins.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ GetIndices(prot_seg, indices_ca);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ best_plane=frame->GetODRPlane(indices_ca);
+ normals.push_back(best_plane.GetNormal());
+ origins.push_back(best_plane.GetOrigin());
+ }
+ return;
+ }
+
std::vector<Real> AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg,
unsigned int stride)
{
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
index d2ba49b6eb2862bb1f8fb7ec4a42da2d564aad26..501638925983170caf5e1b9a991f85900aec60d3 100644
--- a/modules/mol/alg/src/trajectory_analysis.hh
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -45,6 +45,7 @@ namespace ost { namespace mol { namespace alg {
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);
+ void DLLEXPORT_OST_MOL_ALG AnalyzeBestFitPlane(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& normals, geom::Vec3List& origins, unsigned int stride=1);
EntityHandle DLLEXPORT_OST_MOL_ALG CreateMeanStructure(const CoordGroupHandle& traj, const EntityView& selection, int from=0, int to=-1, unsigned int stride=1);
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg, unsigned int stride=1);
}}}//ns
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index bb70ac62ad0acbf8c877433990875ddb8987e01b..c6d13aa78a7ff1ba9b8a6d8920d00e6902cdb30b 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -274,6 +274,16 @@ namespace ost { namespace mol {
}
return atoms_pos_list.GetODRLine();
}
+
+ geom::Plane CoordFrame::GetODRPlane(std::vector<unsigned long>& indices_ca){
+ //Returns the normal to the best fit plane 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.GetODRPlane();
+ }
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
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index db9645dc688ca66db3111d56f158bd9ae0d77766..dbd51c29b6030dfdf687b590824e48678a628fc0 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -61,6 +61,7 @@ public:
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::Plane GetODRPlane(std::vector<unsigned long>& indices_ca);
geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);