diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 7f186c7b102100c58b7599fe017223a52b3ded02..020ef98b6395f2ec74629e6fd4f3c8c9e4a88dc0 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -34,17 +34,39 @@
as a new coord group.
:param frames: The source coord group.
- :param from: index of the first frame
- :param to: index of the last frame plus one. If set to -1, the value is set to
- the number of frames in the coord group
+ :type frames: :class:`~ost.mol.CoordGroupHandle`
:param sel: An entity view containing the selection of atoms to be used for
superposition. If set to an invalid view, all atoms in the coord group are
used.
+ :type sel: :class:`ost.mol.EntityView`
+ :param from: index of the first frame
+ :param to: index of the last frame plus one. If set to -1, the value is set to
+ the number of frames in the coord group
:param ref: The index of the reference frame to use for superposition. If set
to -1, the each frame is superposed to the previous frame.
:returns: A newly created coord group containing the superposed frames.
+
+.. function:: SuperposeFrames(frames, sel, ref_view, from=0, to=-1)
+
+ Same as SuperposeFrames above, but the superposition is done on a reference
+ view and not on another frame of the trajectory.
+
+ :param frames: The source coord group.
+ :type frames: :class:`~ost.mol.CoordGroupHandle`
+ :param sel: An entity view containing the selection of atoms of the frames to be used for
+ superposition.
+ :type sel: :class:`ost.mol.EntityView`
+ :param ref_view: The reference view on which the frames will be superposed. The number
+ of atoms in this reference view should be equal to the number of atoms in sel.
+ :type ref_view: :class:`ost.mol.EntityView`
+ :param from: index of the first frame
+ :param to: index of the last frame plus one. If set to -1, the value is set to
+ the number of frames in the coord group
+
+ :returns: A newly created coord group containing the superposed frames.
+
.. autofunction:: ParseAtomNames
.. autofunction:: MatchResidueByNum
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 48e19d84fc7b35e599454f809c07b27abcf1c461..4f974d940ee2f090eb4d0a5c3d90f46a06616747 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -38,8 +38,12 @@ Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real, co
Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
mol::EntityView (*fc_a)(const mol::EntityView&, Real,bool)=&mol::alg::FilterClashes;
mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterClashes;
+mol::CoordGroupHandle (*superpose_frames1)(mol::CoordGroupHandle&, mol::EntityView&, int, int, int)=&mol::alg::SuperposeFrames;
+mol::CoordGroupHandle (*superpose_frames2)(mol::CoordGroupHandle&, mol::EntityView&, mol::EntityView&, int, int)=&mol::alg::SuperposeFrames;
}
+
+
BOOST_PYTHON_MODULE(_ost_mol_alg)
{
export_svdSuperPose();
@@ -52,7 +56,10 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("LocalDistTest", ldt_b, (arg("ref_index")=0, arg("mdl_index")=1));
def("FilterClashes", fc_a, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
def("FilterClashes", fc_b, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
- def("SuperposeFrames", &ost::mol::alg::SuperposeFrames,
+ def("SuperposeFrames", superpose_frames1,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
arg("end")=-1, arg("ref")=-1));
+ def("SuperposeFrames", superpose_frames2,
+ (arg("source"), arg("sel")=ost::mol::EntityView(), arg("ref_view")=ost::mol::EntityView(),arg("begin")=0,
+ arg("end")=-1));
}
diff --git a/modules/mol/alg/src/superpose_frames.cc b/modules/mol/alg/src/superpose_frames.cc
index aafec9d78f986b69b08f0c2e12c160e48a680bec..465f776e77558f14f2027d981d9c0c339a5c1e1e 100644
--- a/modules/mol/alg/src/superpose_frames.cc
+++ b/modules/mol/alg/src/superpose_frames.cc
@@ -75,9 +75,59 @@ inline EMatX3 col_sub(const EMatX3& m, const ECVec3& s)
return r;
}
-CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
+
+ void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& ref_centered,ECVec3& ref_center,CoordFramePtr frame,std::vector<unsigned long>& indices)
+{
+ // This function superposes and then adds a CoordFrame (frame) to a CoordGroup (superposed).
+ // ref_mat, ref_centered and ref_center contain respectively the positions, centered positions and
+ // vector to center of the reference points for superposition.
+ // indices is a vector of the indices of the atoms to be superposed on the reference positions
+ EMat3X frame_centered=EMat3X::Zero(3, indices.size());
+ EMat3X frame_mat=EMat3X::Zero(3, indices.size());
+ ERVec3 frame_center;
+ for (size_t j=0; j<indices.size(); ++j) {
+ frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
+ }
+ std::vector<geom::Vec3> frame_data=*frame;
+ frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
+ frame_centered=row_sub(frame_mat, frame_center);
+ //single value decomposition
+ Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
+ EMat3 matrixVT=svd.matrixV().transpose();
+ //determine rotation
+ Real detv=matrixVT.determinant();
+ Real dett=svd.matrixU().determinant();
+ Real det=detv*dett;
+ EMat3 e_rot;
+ if (det<0) {
+ EMat3 tmat=EMat3::Identity();
+ tmat(2,2)=-1;
+ e_rot=(svd.matrixU()*tmat)*matrixVT;
+ }else{
+ e_rot=svd.matrixU()*matrixVT;
+ }
+ // prepare rmsd calculation
+ geom::Vec3 shift=rvec_to_gvec(ref_center);
+ geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
+ geom::Mat3 rot=emat_to_gmat(e_rot);
+ geom::Mat4 mat4_com, mat4_rot, mat4_shift;
+ mat4_rot.PasteRotation(rot);
+ mat4_shift.PasteTranslation(shift);
+ mat4_com.PasteTranslation(com_vec);
+ geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
+ for (std::vector<geom::Vec3>::iterator c=frame_data.begin(),
+ e2=frame_data.end(); c!=e2; ++c) {
+ *c=geom::Vec3(tf*geom::Vec4(*c));
+ }
+ superposed.AddFrame(frame_data);
+}
+
+
+CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
int begin, int end, int ref)
{
+ //This function superposes the frames of a CoordGroup (cg) with indices between begin and end
+ //onto the frame with index ref. The superposition is done on a selection of atoms given by the EntityView sel.
int real_end=end==-1 ? cg.GetFrameCount() : end;
CoordFramePtr ref_frame;
std::vector<unsigned long> indices;
@@ -97,10 +147,7 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
- EMat3X frame_centered=EMat3X::Zero(3, indices.size());
- EMat3X frame_mat=EMat3X::Zero(3, indices.size());
ECVec3 ref_center;
- ERVec3 frame_center;
if (ref!=-1) {
ref_frame=cg.GetFrame(ref);
for (size_t i=0; i<indices.size(); ++i) {
@@ -125,46 +172,59 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
ref_centered=col_sub(ref_mat, ref_center);
}
CoordFramePtr frame=cg.GetFrame(i);
- for (size_t j=0; j<indices.size(); ++j) {
- frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
- }
- std::vector<geom::Vec3> frame_data=*frame;
- frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
- frame_centered=row_sub(frame_mat, frame_center);
- //single value decomposition
- Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
- EMat3 matrixVT=svd.matrixV().transpose();
- //determine rotation
- Real detv=matrixVT.determinant();
- Real dett=svd.matrixU().determinant();
- Real det=detv*dett;
- EMat3 e_rot;
- if (det<0) {
- EMat3 tmat=EMat3::Identity();
- tmat(2,2)=-1;
- e_rot=(svd.matrixU()*tmat)*matrixVT;
- }else{
- e_rot=svd.matrixU()*matrixVT;
- }
- // prepare rmsd calculation
- geom::Vec3 shift=rvec_to_gvec(ref_center);
- geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
- geom::Mat3 rot=emat_to_gmat(e_rot);
- geom::Mat4 mat4_com, mat4_rot, mat4_shift;
- mat4_rot.PasteRotation(rot);
- mat4_shift.PasteTranslation(shift);
- mat4_com.PasteTranslation(com_vec);
- geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
- for (std::vector<geom::Vec3>::iterator c=frame_data.begin(),
- e2=frame_data.end(); c!=e2; ++c) {
- *c=geom::Vec3(tf*geom::Vec4(*c));
- }
- superposed.AddFrame(frame_data);
+ AddSuperposedFrame(superposed,ref_mat,ref_centered,ref_center,frame,indices);
if (ref==-1) {
ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
}
}
return superposed;
}
-
+
+
+CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
+ EntityView& ref_view, int begin, int end)
+{
+ //This function superposes the frames of a CoordGroup (cg) with indices between begin and end,
+ //using a selection of atoms (sel), onto an EntityView (ref_view).
+ if (!ref_view.IsValid()){
+ throw std::runtime_error("Invalid reference view");
+ }
+ int real_end=end==-1 ? cg.GetFrameCount() : end;
+ CoordFramePtr ref_frame;
+ std::vector<unsigned long> indices;
+ if (!sel.IsValid()) {
+ indices.reserve(cg.GetAtomCount());
+ for (size_t i=0;i<cg.GetAtomCount(); ++i) {
+ indices.push_back(i);
+ }
+ } else {
+ AtomViewList atoms=sel.GetAtomList();
+ indices.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ indices.push_back(i->GetIndex());
+ }
+ }
+ if (int(indices.size())!=ref_view.GetAtomCount()){
+ throw std::runtime_error("atom counts of the two views are not equal");
+ }
+ CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+ EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
+ EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
+ ECVec3 ref_center;
+ AtomViewList atoms=ref_view.GetAtomList();
+ int i=0;
+ for (AtomViewList::const_iterator a=atoms.begin(),
+ e=atoms.end(); a!=e; a++, i++) {
+ ref_mat.row(i)=gvec_to_cvec((*a).GetPos());
+ }
+ ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+ ref_centered=col_sub(ref_mat, ref_center);
+ for (int i=begin; i<real_end; ++i) {
+ CoordFramePtr frame=cg.GetFrame(i);
+ AddSuperposedFrame(superposed,ref_mat,ref_centered,ref_center,frame,indices);
+ }
+ return superposed;
+}
+
}}}
diff --git a/modules/mol/alg/src/superpose_frames.hh b/modules/mol/alg/src/superpose_frames.hh
index c3475862b9e3d799fb30259e85bf044ef94b5276..8a75ef9aa9d8469d1b55a4518adca96c0f5d80b8 100644
--- a/modules/mol/alg/src/superpose_frames.hh
+++ b/modules/mol/alg/src/superpose_frames.hh
@@ -24,11 +24,15 @@
#include <ost/mol/alg/module_config.hh>
namespace ost { namespace mol { namespace alg {
-/// \brief returns a superposed version of coord group
-CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle cg,
- EntityView sel,
+/// \brief returns a superposed version of coord group, superposed on a reference frame
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle& cg,
+ EntityView& sel,
int begin=0, int end=-1,
int ref=-1);
+/// \brief returns a superposed version of coord group, superposed on a reference view
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle& cg,
+ EntityView& sel, EntityView& ref_view,
+ int begin=0, int end=-1);
}}}
#endif