diff --git a/modules/geom/pymod/export_vecmat3_op.cc b/modules/geom/pymod/export_vecmat3_op.cc
index b2f7c17bf263cc2c2b7acecee3b6960947562883..0cdc3e37358693f63215c0bb8fe05adf17b6617e 100644
--- a/modules/geom/pymod/export_vecmat3_op.cc
+++ b/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);
diff --git a/modules/mol/alg/pymod/export_structure_analysis.cc b/modules/mol/alg/pymod/export_structure_analysis.cc
index 2016cb3053ca67bdbb1c57dc2ea90a9bd751612a..2714f5847ea4e40bd1ec5c63faec3ea317b57336 100644
--- a/modules/mol/alg/pymod/export_structure_analysis.cc
+++ b/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")));
diff --git a/modules/mol/alg/pymod/structure_analysis.py b/modules/mol/alg/pymod/structure_analysis.py
index 3e013e4791a975f7bbb38f7f15021da9992292ac..2a9dda6fcee74c1a12ff5670abb838c864351762 100644
--- a/modules/mol/alg/pymod/structure_analysis.py
+++ b/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
diff --git a/modules/mol/alg/pymod/trajectory_analysis.py b/modules/mol/alg/pymod/trajectory_analysis.py
index 064135c32684bf0fcb6cc068e38673ae432fcb02..23c53b5a70898d4826c393017137b368800984fb 100644
--- a/modules/mol/alg/pymod/trajectory_analysis.py
+++ b/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
diff --git a/modules/mol/alg/src/structure_analysis.cc b/modules/mol/alg/src/structure_analysis.cc
index 960acf47165868dc580161e548ed5c1ad60ef694..b561b0a7fb7d7607217190317c18d643aedd62ea 100644
--- a/modules/mol/alg/src/structure_analysis.cc
+++ b/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];
diff --git a/modules/mol/alg/src/structure_analysis.hh b/modules/mol/alg/src/structure_analysis.hh
index 9cdfa01f58016e2c4baec7fd1af62c9ebc17c5ab..2f450f6924fb9d0f50e558e2da1d54d1dce5606b 100644
--- a/modules/mol/alg/src/structure_analysis.hh
+++ b/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
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index 70ccadd92054bab3e54931c3a362c63e57290619..566f0a43c0c3562830401f339ea62865631d757f 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/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);