Merge branch 'develop' of bc2-dng01.bc2.unibas.ch:/opt/git/GIT_DATA/ost into develop

modules/geom/pymod/export_vecmat3_op.cc

@@ -38,6 +38,7 @@ Real (*Mat3Comp)(const Mat3& m, unsigned int i, unsigned int j)       = &Comp;
 Real (*Mat3Minor)(const Mat3& m, unsigned int i, unsigned int j)      = &Minor;
 Vec3 (*Vec3Min)(const Vec3&, const Vec3&) = &Min;
 Vec3 (*Vec3Max)(const Vec3&, const Vec3&) = &Max;
+Real (*Vec3SignedAngle)(const Vec3& v1, const Vec3& v2, const Vec3& v3) = &SignedAngle;
 Real (*Vec3Distance2WithPBC)(const Vec3&, const Vec3&, const Vec3&)   = &Distance2WithPBC;
 Real (*Vec3DistanceWithPBC)(const Vec3&, const Vec3&, const Vec3&)    = &DistanceWithPBC;
 Vec3 (*wrap_vec3_1)(const Vec3&, const Vec3&, const Vec3&)               = &WrapVec3;
@@ -56,6 +57,7 @@ void export_VecMat3_op()
   def("Distance",Vec3Distance);  
   def("Equal",Mat3Equal, (arg("m1"), arg("m2"), arg("epsilon")=EPSILON));
   def("DihedralAngle", &DihedralAngle);
+  def("SignedAngle", Vec3SignedAngle, (arg("v1"), arg("v2"),arg("ref_normal")));
   def("Dot",Mat3Dot);
   def("Det",Mat3Det);
   def("Cross",Vec3Cross);

modules/mol/alg/pymod/export_structure_analysis.cc

@@ -34,7 +34,7 @@ void export_StructureAnalysis()
   def("CalculateAverageAgreementWithDensityMap",&CalculateAverageAgreementWithDensityMap,(arg("pos_list"),arg("density_map")));
   def("CalculateAgreementWithDensityMap",&CalculateAgreementWithDensityMap,(arg("pos_list"),arg("density_map")));
 #endif
-  def("WrapEntityInPeriodicCell",&WrapEntityInPeriodicCell,(arg("Entity"),arg("cell_center"),arg("ucell_size"),arg("ucell_angles")=geom::Vec3(),arg("group_res")=true));
+  def("WrapEntityInPeriodicCell",&WrapEntityInPeriodicCell,(arg("Entity"),arg("cell_center"),arg("ucell_size"),arg("ucell_angles")=geom::Vec3(),arg("group_res")=true,arg("follow_bonds")=false));
   
 
   def("PariwiseDistanceMatrix",pair_dist2,(arg("EntityView1"),arg("EntityView2")));

modules/mol/alg/pymod/structure_analysis.py

@@ -135,7 +135,7 @@ def CalculateDistanceDifferenceMatrix(sele1,sele2):
   This function calculates the pairwise distance differences between two EntityViews.
   The two EntityViews should have the same number of atoms
   It returns an NxN DistanceDifferenceMatrix M (where N is the number of atoms in sele1)
-  where M[i,j]=(sele2.atoms[i].pos-sele2.atoms[j].pos)-(sele1.atoms[i].pos-sele1.atoms[j].pos)
+  where M[i,j]=|sele2.atoms[i].pos-sele2.atoms[j].pos|-|sele1.atoms[i].pos-sele1.atoms[j].pos|
   """
   try:import numpy as npy
   except ImportError:
@@ -150,10 +150,10 @@ def CalculateDistanceDifferenceMatrix(sele1,sele2):
   n_atoms=sele1.GetAtomCount()
   M=npy.zeros([n_atoms,n_atoms])
   for i,a1 in enumerate(sele1.atoms):
-    for j,a2 in enumerate(sele2.atoms):
+    for j,a2 in enumerate(sele1.atoms):
       if i>=j:continue
-      d1=geom.Distance(a1.pos,a2.pos)
-      d2=geom.Distance(sele2.atoms[i].pos,sele2.atoms[j].pos)
+      d1=ost.geom.Distance(a1.pos,a2.pos)
+      d2=ost.geom.Distance(sele2.atoms[i].pos,sele2.atoms[j].pos)
       M[i,j]=d2-d1
       M[j,i]=d2-d1
   return M

modules/mol/alg/pymod/trajectory_analysis.py

@@ -206,10 +206,35 @@ def AverageDistanceMatrixFromTraj(t,sele,first=0,last=-1):
   M=npy.zeros([n_atoms,n_atoms])
   for i,a1 in enumerate(sele.atoms):
     for j,a2 in enumerate(sele.atoms):
+      if j>i:continue
       d=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
       M[i,j]=npy.mean(d)
       M[j,i]=npy.mean(d)
   return M
 
-
+def AnalyzeDistanceFluctuationMatrix(t,sele,first=0,last=-1):
+  try:
+    import numpy as npy
+  except ImportError:
+    LogError("Function needs numpy, but I could not import it.")
+    raise
+  n_atoms=sele.GetAtomCount()
+  M=npy.zeros([n_atoms,n_atoms])
+  for i,a1 in enumerate(sele.atoms):
+    for j,a2 in enumerate(sele.atoms):
+      if i>j:continue
+      d=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
+      M[j,i]=npy.std(d)
+      M[i,j]=npy.std(d)
+  return M
+  
+def IterativeSuperposition(t,sele,threshold=1.0,initial_sele=None,iterations=5,ref_frame=0):
+  if initial_sele:current_sele=initial_sele
+  else: current_sele=sele
+  for i in range(iterations):
+    t=ost.mol.alg.SuperposeFrames(t,current_sele,ref=ref_frame)
+    al=[a for a in sele.atoms if ost.mol.alg.AnalyzeRMSF(t,ost.mol.CreateViewFromAtoms([a]))<threshold]
+    if len(al)==0:return
+    current_sele=ost.mol.CreateViewFromAtoms(al)
+  return current_sele
     
\ No newline at end of file

modules/mol/alg/src/structure_analysis.cc

@@ -90,7 +90,7 @@ std::vector< std::vector<Real> > PariwiseDistanceMatrix(const EntityView& view){
   
 #endif
 void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom::Vec3 cell_center, const geom::Vec3 ucell_size, \
-                                                    const geom::Vec3 ucell_angles, bool group_residues){
+                                                    const geom::Vec3 ucell_angles, bool group_residues,bool follow_bonds){
   mol::XCSEditor edi=eh.EditXCS(mol::BUFFERED_EDIT);
   bool orthogonal;
   if (ucell_angles==geom::Vec3()){ orthogonal=true; }
@@ -103,14 +103,37 @@ void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom:
       ResidueHandle r=residues[i];
       AtomHandleList atoms=r.GetAtomList();
       geom::Vec3 ref_pos=atoms[0].GetPos();
-      if (orthogonal) {
-        for (AtomHandleList::iterator a=atoms.begin(), e=atoms.end(); a!=e; ++a) {
-          edi.SetAtomPos((*a),geom::WrapVec3((*a).GetPos(),ref_pos,ucell_size));
-        }}
-      else {
-        for (AtomHandleList::iterator a=atoms.begin(), e=atoms.end(); a!=e; ++a) {
-          edi.SetAtomPos((*a),geom::WrapVec3((*a).GetPos(),ref_pos,ucell_size,ucell_angles));
-        }}
+      if (follow_bonds) {
+        int natoms=r.GetAtomCount(),n=1,cycles=1;
+        AtomHandle a=atoms[0];
+        AtomHandleList wrapped_atoms;
+        wrapped_atoms.reserve(natoms);
+        wrapped_atoms.push_back(a);
+        while (n<natoms and cycles<natoms) {
+          ++cycles;
+          for (AtomHandleList::iterator a=wrapped_atoms.begin(), e=wrapped_atoms.end(); a!=e; ++a) {
+            AtomHandleList al=(*a).GetBondPartners();
+            ref_pos=(*a).GetPos();
+            for (AtomHandleList::iterator a2=al.begin(), e2=al.end(); a2!=e2; ++a2) {
+              if (std::find(wrapped_atoms.begin(),wrapped_atoms.end(), *a2) !=wrapped_atoms.end()){
+                continue;
+              } else {
+                ++n;
+                wrapped_atoms.push_back(*a2);
+                if (orthogonal) {
+                    edi.SetAtomPos((*a2),geom::WrapVec3((*a2).GetPos(),ref_pos,ucell_size));
+                } else {
+                  edi.SetAtomPos((*a2),geom::WrapVec3((*a2).GetPos(),ref_pos,ucell_size,ucell_angles));
+          }}}}}
+      } else {
+        if (orthogonal) {
+          for (AtomHandleList::iterator a=atoms.begin(), e=atoms.end(); a!=e; ++a) {
+            edi.SetAtomPos((*a),geom::WrapVec3((*a).GetPos(),ref_pos,ucell_size));
+          }}
+        else {
+          for (AtomHandleList::iterator a=atoms.begin(), e=atoms.end(); a!=e; ++a) {
+            edi.SetAtomPos((*a),geom::WrapVec3((*a).GetPos(),ref_pos,ucell_size,ucell_angles));
+        }}}
     }
     for (unsigned int i=0; i<n_residues; ++i) {
       ResidueHandle r=residues[i];

modules/mol/alg/src/structure_analysis.hh

@@ -37,7 +37,7 @@ namespace ost { namespace mol { namespace alg {
   std::vector<Real> DLLEXPORT_OST_MOL_ALG CalculateAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map);
   Real DLLEXPORT_OST_MOL_ALG CalculateAverageAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map);
 #endif
-  void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom::Vec3 cell_center, const geom::Vec3 ucell_size, const geom::Vec3 ucell_angles=geom::Vec3(), bool group_res=true);
+  void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom::Vec3 cell_center, const geom::Vec3 ucell_size, const geom::Vec3 ucell_angles=geom::Vec3(), bool group_res=true, bool follow_bond=true);
   std::vector< std::vector<Real> > DLLEXPORT_OST_MOL_ALG PariwiseDistanceMatrix(const EntityView& view1, const EntityView& view2);
   std::vector< std::vector<Real> > DLLEXPORT_OST_MOL_ALG PariwiseDistanceMatrix(const EntityView& view);
 }}}//ns

modules/mol/base/pymod/export_coord_frame.cc

@@ -43,7 +43,11 @@ Real (CoordFrame::*get_min_dist)(const mol::EntityView&, const mol::EntityView&)
 Real (CoordFrame::*get_alpha)(const mol::EntityView&) const = &CoordFrame::GetAlphaHelixContent;
 geom::Line3 (CoordFrame::*get_odr_line)(const mol::EntityView&) const = &CoordFrame::GetODRLine;
 geom::Plane (CoordFrame::*get_odr_plane)(const mol::EntityView&) const = &CoordFrame::GetODRPlane;
-std::pair<geom::Line3,Real> (CoordFrame::*fit_cylinder)(const mol::EntityView&) const = &CoordFrame::FitCylinder;
+boost::python::tuple wrap_FitCylinder(const CoordFrame& f,const mol::EntityView& ev) {
+  std::pair<geom::Line3,Real> pair=f.FitCylinder(ev);
+  return boost::python::make_tuple<geom::Line3,Real>(pair.first,pair.second);
+}  
+
 // TODO: move to geom
 geom::Line3 (CoordFrame::*get_odr_line2)() const = &geom::Vec3List::GetODRLine;
 
@@ -83,7 +87,7 @@ void export_CoordFrame()
     .def("GetODRLine",get_odr_line)
     .def("GetODRLine",get_odr_line2)
     .def("GetAlphaHelixContent",get_alpha)
-    .def("FitCylinder",fit_cylinder)
+    .def("FitCylinder",wrap_FitCylinder)
   ;
   def("CreateCoordFrame",create_coord_frame1);
   def("CreateCoordFrame",create_coord_frame2);