diff --git a/modules/geom/src/composite3_op.cc b/modules/geom/src/composite3_op.cc
index 672896f0d9aaec18fad77a24aa8065c7c43e34e4..87f799d516ab3be540ee287dbb2ceed98fee6d4a 100644
--- a/modules/geom/src/composite3_op.cc
+++ b/modules/geom/src/composite3_op.cc
@@ -178,12 +178,18 @@ bool IsInSphere(const Sphere& s, const Vec3& v){
return Length(s.GetOrigin()-v)<=s.GetRadius();
}
-Line3 Vec3List::FitCylinder(){
- Line3 axis=this->GetODRLine(),axis_old;
- Real radius,res_sum_old,res_sum,delta=0.01,prec=0.001,err;
+Line3 Vec3List::FitCylinder(const Vec3 initial_direction, const Vec3 center){
+ //This function fits a cylinder to the positions in Vec3List
+ //It takes as argument an initial guess for the direction and the geometric
+ //center of the atoms. The center is not changed during optimisation as the
+ //best fitting cylinder can be shown to have its axis pass through the geometric center
+ Line3 axis=Line3(center,center+initial_direction), axis_old;
+ Real radius,res_sum_old,res_sum,delta_0=0.01,prec=0.0000001,err,norm,delta;
unsigned long n_step=1000, n_res=this->size();
Vec3 v,gradient;
+
radius=0.0;
+ delta=delta_0;
for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
radius+=geom::Distance(axis,(*i));
}
@@ -197,6 +203,10 @@ Line3 Vec3List::FitCylinder(){
while (err>prec and k<n_step) {
res_sum_old=res_sum;
axis_old=axis;
+ radius=0.0;
+ if (k>50) {
+ delta=delta_0*pow((50./k),2.0);
+ }
for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
radius+=Distance(axis,(*i));
}
@@ -206,13 +216,26 @@ Line3 Vec3List::FitCylinder(){
v=Vec3(0.0,0.0,0.0);
v[j]=delta;
axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
+ radius=0.0;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ radius+=Distance(axis,(*i));
+ }
+ radius/=Real(n_res);
for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
res_sum+=pow(Distance(axis,(*i))-radius,2.);
}
- gradient[j]=res_sum-res_sum_old;
+ gradient[j]=(res_sum-res_sum_old)/delta;
+ }
+ norm=Dot(gradient,gradient);
+ if (norm>1.) {
+ gradient=Normalize(gradient);
}
- gradient=Normalize(gradient);
axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
+ radius=0.0;
+ for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+ radius+=Distance(axis,(*i));
+ }
+ radius/=Real(n_res);
res_sum=0.0;
for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
res_sum+=pow(Distance(axis,(*i))-radius,2.);
@@ -220,6 +243,9 @@ Line3 Vec3List::FitCylinder(){
err=fabs((res_sum-res_sum_old)/float(n_res));
k++;
}
+ if (err>prec) {
+ std::cout<<"axis fitting did not converge"<<std::endl;
+ }
return axis;
}
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index 2a27d500a8a0c5ca9c66127de47f0cf37f48f932..a2d1d3cea92d4fb79dfb71650eaa17fa1cadce1a 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -218,7 +218,7 @@ public:
Mat3 GetPrincipalAxes() const;
Line3 GetODRLine();
- Line3 FitCylinder();
+ Line3 FitCylinder(const Vec3 initial_direction, const Vec3 center);
};
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index 2d0e6f32f9d5873dd16ed823c6ce4a8d01ac024f..c7222e599fc043cc57ddc94dca9495bd664c175a 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -233,6 +233,9 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
void AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
geom::Vec3List& centers, unsigned int stride)
+ //This function calculates the best fitting cylinder to the C-alpha atoms of an EntityView and returns
+ //the geometric center as well as the axis of that cylinder. We take care to have the axis point towards
+ //the last residue of the selection, usually the direction of the alpha-helix
{
CheckHandleValidity(traj);
if (prot_seg.GetAtomCount()==0){
@@ -240,13 +243,17 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
}
std::vector<unsigned long> indices_ca;
geom::Line3 axis;
+ Real sign;
directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
GetCaIndices(prot_seg, indices_ca);
+ unsigned int n_atoms=indices_ca.size();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
axis=frame->FitCylinder(indices_ca);
- directions.push_back(axis.GetDirection());
+ sign=geom::Dot(axis.GetDirection(),(*frame)[indices_ca[n_atoms-1]]-axis.GetOrigin());
+ sign=sign/fabs(sign);
+ directions.push_back(sign*axis.GetDirection());
centers.push_back(axis.GetOrigin());
}
return;
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index efeb37e20fc086ebb425751b258e74334f86d3d1..bb70ac62ad0acbf8c877433990875ddb8987e01b 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -277,12 +277,29 @@ namespace ost { namespace mol {
geom::Line3 CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca){
//Returns a line which is the axis of a fitted cylinder to the atoms with indices given in indices_ca
+ //It is assumed that we fit an alpha-helix and that the CA atoms are oredered sequentially
+ //This is used for the initial guess of the helix axis
geom::Vec3List atoms_pos_list;
- atoms_pos_list.reserve(indices_ca.size());
+ int n_atoms=indices_ca.size();
+ atoms_pos_list.reserve(n_atoms);
for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
atoms_pos_list.push_back((*this)[*i1]);
}
- return atoms_pos_list.FitCylinder();
+ //Initial guess
+ geom::Vec3 initial_axis=geom::Vec3(0.0,0.0,0.0),center=geom::Vec3(0.0,0.0,0.0);
+ if (n_atoms<5) {
+ initial_axis=atoms_pos_list[n_atoms-1]-atoms_pos_list[0];
+ }
+ else {
+ for (geom::Vec3List::const_iterator i=atoms_pos_list.begin(),e=atoms_pos_list.end()-4; i!=e; ++i) {
+ initial_axis+=(*(i+4))-(*i);
+ }
+ }
+ for (geom::Vec3List::const_iterator i=atoms_pos_list.begin(),e=atoms_pos_list.end(); i!=e; ++i) {
+ center+=(*i);
+ }
+ center/=atoms_pos_list.size();
+ return atoms_pos_list.FitCylinder(initial_axis,center);
}
Real CoordFrame::GetAlphaHelixContent(const mol::EntityView& segment){