Skip to content
Snippets Groups Projects
Commit 122c97ea authored by Niklaus Johner's avatar Niklaus Johner
Browse files

Added a function to calculate an average structure from a trajectory

parent 6fb3c0dc
No related branches found
No related tags found
No related merge requests found
......@@ -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 ));
}
......@@ -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
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment