diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc index b2c9dcd562559563b939512de48bcdfdf4db4e72..5671dbca091c49746976714679c04dd48a68fd36 100644 --- a/modules/mol/alg/pymod/export_trajectory_analysis.cc +++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc @@ -49,4 +49,5 @@ void export_TrajectoryAnalysis() 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("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 1688eed33212fefbafa35758a7848101e61d515e..7f81b91d13a68edaff14e5ea2fd3a20806443192 100644 --- a/modules/mol/alg/src/trajectory_analysis.cc +++ b/modules/mol/alg/src/trajectory_analysis.cc @@ -25,7 +25,9 @@ #include <ost/geom/vec3.hh> #include <ost/base.hh> #include <ost/geom/geom.hh> +#include <ost/gfx/entity.hh> #include <ost/mol/entity_view.hh> +#include <ost/mol/view_op.hh> #include <ost/mol/coord_group.hh> #include "trajectory_analysis.hh" @@ -293,6 +295,41 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c return helicity; } - + //This function constructs mean structures from a trajectory + EntityHandle CreateMeanStructure(const CoordGroupHandle& traj, const EntityView& selection, + int from, int to, unsigned int stride) + { + CheckHandleValidity(traj); + if (to==-1)to=traj.GetFrameCount(); + unsigned int n_atoms=selection.GetAtomCount(); + if (to<from) { + throw std::runtime_error("to smaller than from"); + } + unsigned int n_frames=to-from; + if (n_atoms==0){ + throw std::runtime_error("EntityView is empty"); + } + if (n_frames<=stride) { + throw std::runtime_error("number of frames is too small"); + } + std::vector<unsigned long> indices; + std::vector<geom::Vec3> mean_positions; + EntityHandle eh; + eh=CreateEntityFromView(selection,1); + GetIndices(selection,indices); + mean_positions.assign(n_atoms,geom::Vec3(0.0,0.0,0.0)); + for (int i=from; i<to; i+=stride) { + CoordFramePtr frame=traj.GetFrame(i); + for (unsigned int j=0; j<n_atoms; ++j) { + mean_positions[j]+=(*frame)[indices[j]]; + } + } + mol::XCSEditor edi=eh.EditXCS(mol::BUFFERED_EDIT); + AtomHandleList atoms=eh.GetAtomList(); + for (unsigned int j=0; j<n_atoms; ++j) { + edi.SetAtomPos(atoms[j],mean_positions[j]/Real(n_frames)); + } + return eh; + } }}} //ns diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh index 5c2ec6de2417b461eaf137da05874e62155c241f..d2ba49b6eb2862bb1f8fb7ec4a42da2d564aad26 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); + 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 #endif