From 45c7f9842352c8f7c6cfe524f3ab6c228de3f14f Mon Sep 17 00:00:00 2001 From: Niklaus Johner <nij2003@med.cornell.edu> Date: Wed, 6 Jul 2011 15:31:24 -0400 Subject: [PATCH] Added functions to calculate the minimal distance between two sets of atoms or between the center of mass of one set and the atoms in the other set. This function was added in geom::vecmat3_op and then on the CoordFrame and trajectory_analysis level. --- modules/geom/doc/composite.rst | 15 +++- modules/geom/src/fnc.hh | 3 +- modules/geom/src/vecmat3_op.cc | 15 ++++ modules/geom/src/vecmat3_op.hh | 3 + modules/geom/tests/test_op3.cc | 14 +++ modules/mol/alg/doc/molalg.rst | 42 ++++++++- .../alg/pymod/export_trajectory_analysis.cc | 3 + modules/mol/alg/src/trajectory_analysis.cc | 90 +++++++++++++++++-- modules/mol/alg/src/trajectory_analysis.hh | 5 +- modules/mol/base/pymod/export_coord_frame.cc | 3 +- modules/mol/base/src/coord_frame.cc | 41 ++++++++- modules/mol/base/src/coord_frame.hh | 7 +- 12 files changed, 226 insertions(+), 15 deletions(-) diff --git a/modules/geom/doc/composite.rst b/modules/geom/doc/composite.rst index 0271a206b..863763ccb 100644 --- a/modules/geom/doc/composite.rst +++ b/modules/geom/doc/composite.rst @@ -526,4 +526,17 @@ Operations on Geometrical Objects Check whether the `sphere` contains `point`. - :rtype: bool \ No newline at end of file + :rtype: bool + + +.. function:: MinDistance(list1, list2) + + Calculate the minimal distance between two sets of points + + :param list1: the first set of points + :type list1: :class:`~ost.geom.Vec3List` + + :param list2: the second set of points + :type list2: :class:`~ost.geom.Vec3List` + + \ No newline at end of file diff --git a/modules/geom/src/fnc.hh b/modules/geom/src/fnc.hh index 742fb30eb..0062db1f7 100644 --- a/modules/geom/src/fnc.hh +++ b/modules/geom/src/fnc.hh @@ -40,7 +40,8 @@ Real Distance(const Vec3& p1, const Vec3& p2); Real Distance(const Vec3& p, const Line& l); //! returns the distance between a point and a line Real Distance(const Line& l, const Vec3& p); - +//! returns the minimal distance between the points in two Vec3List +Real MinDistance(const Vec3List& l1, const Vec3List& l2); } // ns geom #endif diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc index c4b0a6c14..7fb8e2fef 100644 --- a/modules/geom/src/vecmat3_op.cc +++ b/modules/geom/src/vecmat3_op.cc @@ -192,4 +192,19 @@ Real DihedralAngle(const Vec3& p1, const Vec3& p2, const Vec3& p3, Dot(r12cross, r23cross)); } +Real MinDistance(const Vec3List& l1, const Vec3List& l2) +{ + // returns the minimal distance between two sets of points (Vec3List) + if (l1.size()==0 || l2.size()==0){throw std::runtime_error("cannot calculate minimal distance: empty Vec3List");} + Real min=Distance(*l1.begin(),*l2.begin()); + Real d; + for (Vec3List::const_iterator p1=l1.begin(),e1=l1.end(); p1!=e1; p1++) { + for (Vec3List::const_iterator p2=l2.begin(),e2=l2.end(); p2!=e2; p2++) { + d=Distance(*p1,*p2); + if (d<min) min=d; + } + } + return min; +} + } // ns diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh index 39d0f7117..373b9561b 100644 --- a/modules/geom/src/vecmat3_op.hh +++ b/modules/geom/src/vecmat3_op.hh @@ -194,6 +194,9 @@ inline Real Distance(const Vec3& p1, const Vec3& p2) return Length(p1-p2); } +//! returns the minimal distance between the points in two Vec3List +Real MinDistance(const Vec3List& l1, const Vec3List& l2); + } // ns #endif diff --git a/modules/geom/tests/test_op3.cc b/modules/geom/tests/test_op3.cc index 9672ff46a..a2e30db8b 100644 --- a/modules/geom/tests/test_op3.cc +++ b/modules/geom/tests/test_op3.cc @@ -48,4 +48,18 @@ BOOST_AUTO_TEST_CASE(vecmat_mult3) BOOST_CHECK(Equal(v3*m2,v2*m3)); } +BOOST_AUTO_TEST_CASE(vec3list_op) +{ + Vec3List l1,l2; + BOOST_CHECK_THROW(MinDistance(l1,l2),std::runtime_error); + l1.push_back(Vec3(0.0,0.0,0.0)); + l1.push_back(Vec3(1.0,0.0,0.0)); + BOOST_CHECK_THROW(MinDistance(l1,l2),std::runtime_error); + BOOST_CHECK_THROW(MinDistance(l2,l1),std::runtime_error); + l2.push_back(Vec3(4.0,4.0,0.0)); + l2.push_back(Vec3(2.0,4.0,5.0)); + l2.push_back(Vec3(0.0,3.0,6.0)); + BOOST_CHECK(MinDistance(l1,l2)==5.0); +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst index 4730eec07..7f186c7b1 100644 --- a/modules/mol/alg/doc/molalg.rst +++ b/modules/mol/alg/doc/molalg.rst @@ -182,7 +182,7 @@ used to skip frames in the analysis. consecutive frames analyzed. -.. function:: AnalyzeRMSD(traj, reference_view, sele_view) +.. function:: AnalyzeRMSD(traj, reference_view, sele_view, stride=1) 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 @@ -210,9 +210,49 @@ used to skip frames in the analysis. +.. function:: AnalyzeMinDistance(traj, view1, view2, stride=1) + This function extracts the minimal distance between two sets of atoms + (view1 and view2) for each frame in a trajectory and returns it as a vector. + :param traj: The trajectory to be analyzed. + :type traj: :class:`~ost.mol.CoordGroupHandle` + :param view1: The first group of atoms + :type view1: :class:`~ost.mol.EntityView`. + :param view2: The second group of atoms + :type view2: :class:`~ost.mol.EntityView`. + :param stride: Size of the increment of the frame's index between two + consecutive frames analyzed. + +.. function:: AnalyzeMinDistanceBetwCenterOfMassAndView(traj, view_cm, view_atoms, stride=1) + + This function extracts the minimal distance between a set of atoms + (view_atoms) and the center of mass of a second set of atoms (view_cm) + for each frame in a trajectory and returns it as a vector. + + :param traj: The trajectory to be analyzed. + :type traj: :class:`~ost.mol.CoordGroupHandle` + :param view_cm: The group of atoms from which the center of mass is taken + :type view_cm: :class:`~ost.mol.EntityView`. + :param view_atoms: The second group of atoms + :type view_atoms: :class:`~ost.mol.EntityView`. + :param stride: Size of the increment of the frame's index between two + consecutive frames analyzed. +.. function:: AnalyzeAromaticRingInteraction(traj, view_ring1, view_ring2, stride=1) + This function is a crude analysis of aromatic ring interactions. For each frame in a trajectory, it calculates + the minimal distance between the atoms in one view and the center of mass of the other + and vice versa, and returns the minimum between these two minimal distances. + Concretely, if the two views are the heavy atoms of two rings, then it returns the minimal + center of mass - heavy atom distance betweent he two rings + :param traj: The trajectory to be analyzed. + :type traj: :class:`~ost.mol.CoordGroupHandle` + :param view_ring1: First group of atoms + :type view_ring1: :class:`~ost.mol.EntityView`. + :param view_ring2: Second group of atoms + :type view_ring2: :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 61e153bb3..31377016a 100644 --- a/modules/mol/alg/pymod/export_trajectory_analysis.cc +++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc @@ -43,4 +43,7 @@ void export_TrajectoryAnalysis() 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)); + 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)); } diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc index c8660fc99..012b77b32 100644 --- a/modules/mol/alg/src/trajectory_analysis.cc +++ b/modules/mol/alg/src/trajectory_analysis.cc @@ -29,7 +29,6 @@ #include <ost/mol/coord_group.hh> #include "trajectory_analysis.hh" - namespace ost { namespace mol { namespace alg { geom::Vec3List AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride) @@ -92,7 +91,7 @@ std::vector<Real> AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomH return ang; } -geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride) +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. { @@ -101,7 +100,7 @@ geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const Entity pos.reserve(ceil(traj.GetFrameCount()/float(stride))); std::vector<unsigned long> indices; std::vector<Real> masses; - GetIndicesAndMasses(Sele, indices, masses); + GetIndicesAndMasses(sele, indices, masses); for (size_t i=0; i<traj.GetFrameCount(); i+=stride) { CoordFramePtr frame=traj.GetFrame(i); pos.push_back(frame->GetCenterOfMassPos(indices,masses)); @@ -109,8 +108,8 @@ geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const Entity return pos; } -std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& Sele1, - const EntityView& Sele2, unsigned int stride) +std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& sele1, + const EntityView& sele2, unsigned int stride) //This function extracts the distance between the CenterOfMass of two selection (entity views) from a trajectory //and returns it as a vector. { @@ -119,8 +118,8 @@ std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, dist.reserve(ceil(traj.GetFrameCount()/float(stride))); std::vector<unsigned long> indices1,indices2; std::vector<Real> masses1,masses2; - GetIndicesAndMasses(Sele1, indices1, masses1); - GetIndicesAndMasses(Sele2, indices2, masses2); + GetIndicesAndMasses(sele1, indices1, masses1); + GetIndicesAndMasses(sele2, indices2, masses2); for (size_t i=0; i<traj.GetFrameCount(); i+=stride) { CoordFramePtr frame=traj.GetFrame(i); dist.push_back(frame->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2)); @@ -153,5 +152,82 @@ std::vector<Real> AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& re } return rmsd; } + +std::vector<Real> AnalyzeMinDistance(const CoordGroupHandle& traj, const EntityView& view1, + const EntityView& view2,unsigned int stride) +// This function extracts the minimal distance between two sets of atoms (view1 and view2) for +// each frame in a trajectory and returns it as a vector. + { + CheckHandleValidity(traj); + if (view1.GetAtomCount()==0){ + throw std::runtime_error("first EntityView is empty"); + } + if (view2.GetAtomCount()==0){ + throw std::runtime_error("second EntityView is empty"); + } + std::vector<Real> dist; + dist.reserve(ceil(traj.GetFrameCount()/float(stride))); + std::vector<unsigned long> indices1,indices2; + GetIndices(view1, indices1); + GetIndices(view2, indices2); + for (size_t i=0; i<traj.GetFrameCount(); i+=stride) { + CoordFramePtr frame=traj.GetFrame(i); + dist.push_back(frame->GetMinDistance(indices1,indices2)); + } + return dist; +} +std::vector<Real> AnalyzeMinDistanceBetwCenterOfMassAndView(const CoordGroupHandle& traj, const EntityView& view_cm, + const EntityView& view_atoms,unsigned int stride) + // This function extracts the minimal distance between a set of atoms (view_atoms) and the center of mass + // of a second set of atoms (view_cm) for each frame in a trajectory and returns it as a vector. + { + CheckHandleValidity(traj); + if (view_cm.GetAtomCount()==0){ + throw std::runtime_error("first EntityView is empty"); + } + if (view_atoms.GetAtomCount()==0){ + throw std::runtime_error("second EntityView is empty"); + } + std::vector<Real> dist, masses_cm; + dist.reserve(ceil(traj.GetFrameCount()/float(stride))); + std::vector<unsigned long> indices_cm,indices_atoms; + GetIndicesAndMasses(view_cm, indices_cm,masses_cm); + GetIndices(view_atoms, indices_atoms); + for (size_t i=0; i<traj.GetFrameCount(); i+=stride) { + CoordFramePtr frame=traj.GetFrame(i); + dist.push_back(frame->GetMinDistBetwCenterOfMassAndView(indices_cm, masses_cm, indices_atoms)); + } + return dist; + } + +std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1, + const EntityView& view_ring2,unsigned int stride) + // This function is a crude analysis of aromatic ring interactions. For each frame in a trajectory, it calculates + // the minimal distance between the atoms in one view and the center of mass of the other + // and vice versa, and returns the minimum between these two minimal distances. + // Concretely, if the two views are the heavy atoms of two rings, then it returns the minimal + // center of mass - heavy atom distance betweent he two rings + { + CheckHandleValidity(traj); + if (view_ring1.GetAtomCount()==0){ + throw std::runtime_error("first EntityView is empty"); + } + if (view_ring2.GetAtomCount()==0){ + throw std::runtime_error("second EntityView is empty"); + } + std::vector<Real> dist, masses_ring1,masses_ring2; + dist.reserve(ceil(traj.GetFrameCount()/float(stride))); + std::vector<unsigned long> indices_ring1,indices_ring2; + Real d1,d2; + GetIndicesAndMasses(view_ring1, indices_ring1,masses_ring1); + GetIndicesAndMasses(view_ring2, indices_ring2,masses_ring2); + for (size_t i=0; i<traj.GetFrameCount(); i+=stride) { + CoordFramePtr frame=traj.GetFrame(i); + d1=frame->GetMinDistBetwCenterOfMassAndView(indices_ring1, masses_ring1, indices_ring2); + d2=frame->GetMinDistBetwCenterOfMassAndView(indices_ring2, masses_ring2, indices_ring1); + dist.push_back(std::min(d1,d2)); + } + return dist; + } }}} //ns diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh index d2947e63b..faa6c8ea0 100644 --- a/modules/mol/alg/src/trajectory_analysis.hh +++ b/modules/mol/alg/src/trajectory_analysis.hh @@ -40,7 +40,8 @@ namespace ost { namespace mol { namespace alg { 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); - - + 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); }}}//ns #endif diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc index 16665e0bb..fdb41fce6 100644 --- a/modules/mol/base/pymod/export_coord_frame.cc +++ b/modules/mol/base/pymod/export_coord_frame.cc @@ -33,7 +33,7 @@ Real (CoordFrame::*get_dihedral)(const AtomHandle& a1, const AtomHandle& a2, con geom::Vec3 (CoordFrame::*get_cm)(const mol::EntityView& sele) = &CoordFrame::GetCenterOfMassPos; Real (CoordFrame::*get_dist_cm)(const mol::EntityView& sele1, const mol::EntityView& sele2) = &CoordFrame::GetDistanceBetwCenterOfMass; Real (CoordFrame::*get_rmsd)(const mol::EntityView& Reference_View, const mol::EntityView& sele_View) = &CoordFrame::GetRMSD; - +Real (CoordFrame::*get_min_dist)(const mol::EntityView& view1, const mol::EntityView& view2) = &CoordFrame::GetMinDistance; void export_CoordFrame() { @@ -45,6 +45,7 @@ void export_CoordFrame() .def("GetCenterOfMassPos", get_cm) .def("GetDistanceBetwCenterOfMass", get_dist_cm) .def("GetRMSD",get_rmsd) + .def("GetMinDistance",get_min_dist) ; def("CreateCoordFrame",CreateCoordFrame); } diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc index 80105c128..80ef4f58f 100644 --- a/modules/mol/base/src/coord_frame.cc +++ b/modules/mol/base/src/coord_frame.cc @@ -180,4 +180,43 @@ namespace ost { namespace mol { return this->GetRMSD(ref_pos,indices_sele); } -}}//ns \ No newline at end of file + Real CoordFrame::GetMinDistance(std::vector<unsigned long>& index_list1, std::vector<unsigned long>& index_list2){ + //Returns the minimal distance between two groups of atoms with indices given by index_list1 and index_list2 + geom::Vec3List pos_list1,pos_list2; + for (std::vector<unsigned long>::const_iterator i1=index_list1.begin(),e=index_list1.end(); i1!=e; ++i1) { + pos_list1.push_back((*this)[*i1]); + } + for (std::vector<unsigned long>::const_iterator i1=index_list2.begin(),e=index_list2.end(); i1!=e; ++i1) { + pos_list2.push_back((*this)[*i1]); + } + return geom::MinDistance(pos_list1,pos_list2); + } + Real CoordFrame::GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2){ + //Returns the minimal distance between the atoms of two views (view1 and view2) + std::vector<unsigned long> index_list1,index_list2; + GetIndices(view1,index_list1); + GetIndices(view2,index_list2); + return this->GetMinDistance(index_list1,index_list2); + } + + Real CoordFrame::GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm, std::vector<Real>& masses_cm, + std::vector<unsigned long>& indices_atoms){ + //Returns the minimal distance between the center of mass of a group of atoms (with indices given in indices_cm and masses + //in masses_cm) and the atoms of another group (with indices given in indices_atoms) + geom::Vec3List cm_pos,atoms_pos_list; + cm_pos.push_back(this->GetCenterOfMassPos(indices_cm, masses_cm)); + for (std::vector<unsigned long>::const_iterator i1=indices_atoms.begin(),e=indices_atoms.end(); i1!=e; ++i1) { + atoms_pos_list.push_back((*this)[*i1]); + } + return geom::MinDistance(cm_pos,atoms_pos_list); + } + Real CoordFrame::GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm, const mol::EntityView& view_atoms){ + //Returns the minimal distance between the center of mass of a view (view_cm) and the atoms of another view (view_atoms) + std::vector<unsigned long> indices_cm,indices_atoms; + std::vector<Real> masses_cm; + GetIndices(view_atoms,indices_atoms); + GetIndicesAndMasses(view_cm,indices_cm,masses_cm); + return this->GetMinDistBetwCenterOfMassAndView(indices_cm,masses_cm,indices_atoms); + } +}}//ns + diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh index c856e583c..8b34cc52a 100644 --- a/modules/mol/base/src/coord_frame.hh +++ b/modules/mol/base/src/coord_frame.hh @@ -54,7 +54,12 @@ public: 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); + Real GetRMSD(const mol::EntityView& reference_view, const mol::EntityView& sele_view); + Real GetMinDistance(std::vector<unsigned long>& index_list1, std::vector<unsigned long>& index_list2); + Real GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2); + 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); }; void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices); -- GitLab