diff --git a/modules/geom/src/composite3_op.cc b/modules/geom/src/composite3_op.cc
index 87f799d516ab3be540ee287dbb2ceed98fee6d4a..09fb2c72e7430c2fc938534ff4167052485a4b24 100644
--- a/modules/geom/src/composite3_op.cc
+++ b/modules/geom/src/composite3_op.cc
@@ -178,76 +178,5 @@ bool IsInSphere(const Sphere& s, const Vec3& v){
return Length(s.GetOrigin()-v)<=s.GetRadius();
}
-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));
- }
- 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.);
- }
- unsigned long k=0;
- err=2.0*prec;
- 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));
- }
- radius/=Real(n_res);
- for (int j=0; j!=3; ++j){
- res_sum=0.0;
- 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)/delta;
- }
- norm=Dot(gradient,gradient);
- if (norm>1.) {
- 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.);
- }
- 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;
-}
-
} // ns
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index ea07f3cca30b46f4aeef0139d7f4a0a0aab9345b..6e657081fb7ed1496d85514b29f390aa52d6028e 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -19,8 +19,12 @@
#include <algorithm>
#include <Eigen/SVD>
-#include "vec3.hh"
+#include "vec3.hh"
+// TODO: these are for the (misplaced) Vec3List algorithm functions
+#include "vecmat3_op.hh"
+#include "composite3.hh"
+#include "composite3_op.hh"
namespace geom {
@@ -68,18 +72,86 @@ Vec3 Vec3List::GetCenter() const
return center/=this->size();
}
-Line3 Vec3List::GetODRLine()
+Line3 Vec3List::GetODRLine() const
{
Vec3 center=this->GetCenter();
Vec3 direction=this->GetPrincipalAxes().GetRow(2);
return Line3(center,center+direction);
}
-Plane Vec3List::GetODRPlane()
+Plane Vec3List::GetODRPlane() const
{
Vec3 origin=this->GetCenter();
Vec3 normal=this->GetPrincipalAxes().GetRow(0);
return Plane(origin,normal);
}
+Line3 Vec3List::FitCylinder(const Vec3& initial_direction, const Vec3& center) const
+{
+ 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));
+ }
+ 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.);
+ }
+ unsigned long k=0;
+ err=2.0*prec;
+ 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));
+ }
+ radius/=Real(n_res);
+ for (int j=0; j!=3; ++j){
+ res_sum=0.0;
+ 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)/delta;
+ }
+ norm=Dot(gradient,gradient);
+ if (norm>1.) {
+ 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.);
+ }
+ 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 ac3458837fc20245ac7d185a56a58cbf69b8eb08..7b4bb5b33fd33cccfe285f6d602bfa8d968bb687 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -189,15 +189,13 @@ inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
os << "[" << v.x << ", " << v.y << ", " << v.z << "]";
return os;
}
-}
-
-#include <ost/geom/vec2.hh>
-#include <ost/geom/vec4.hh>
-#include <ost/geom/mat3.hh>
-#include <ost/geom/composite3.hh>
+} // ns geom
namespace geom {
+ // TODO: move to separate file
+ class Mat3;
+
class DLLEXPORT_OST_GEOM Vec3List : public std::vector<Vec3> {
public:
typedef std::vector<Vec3> base_type;
@@ -213,30 +211,40 @@ public:
base_type::operator=(rhs);
return *this;
}
+
+ // TODO: move some or all of these to stand-alone functions
Mat3 GetInertia() const;
-
Vec3 GetCenter() const;
-
Mat3 GetPrincipalAxes() const;
- Line3 GetODRLine();
- Plane GetODRPlane();
- Line3 FitCylinder(const Vec3 initial_direction, const Vec3 center);
-
+ Line3 GetODRLine() const;
+ Plane GetODRPlane() const;
+
+ //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 FitCylinder(const Vec3& initial_direction, const Vec3& center) const;
};
+} // ns geom
-inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
+#include <ost/geom/vec2.hh>
+#include <ost/geom/vec4.hh>
+#include <ost/geom/mat3.hh>
+#include <ost/geom/composite3.hh>
-inline Vec3::Vec3(const Vec4& v): x(v.x), y(v.y), z(v.z)
-{
- if (std::fabs(v.w)<1e-10) {
- throw DivideByZeroException();
- }
- x/=v.w;
- y/=v.w;
- z/=v.w;
-}
+namespace geom {
+ inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
+ inline Vec3::Vec3(const Vec4& v): x(v.x), y(v.y), z(v.z)
+ {
+ if (std::fabs(v.w)<1e-10) {
+ throw DivideByZeroException();
+ }
+ x/=v.w;
+ y/=v.w;
+ z/=v.w;
+ }
} // namespace geom
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index 15a290279063e9c2dd7dc9948f55c0d0d70437e6..5acc5bf4b79dcb0eca650ebe306f6d944cc78b68 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/modules/io/src/mol/dcd_io.cc
@@ -67,7 +67,7 @@ bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
}
bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
- bool& skip_flag, bool& gap_flag)
+ bool& ucell_flag, bool& gap_flag)
{
if (!istream) {
return false;
@@ -75,7 +75,7 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
char dummy[4];
gap_flag=true;
swap_flag=false;
- skip_flag=false;
+ ucell_flag=false;
if(gap_flag) istream.read(dummy,sizeof(dummy));
istream.read(header.hdrr,sizeof(char)*4);
istream.read(reinterpret_cast<char*>(header.icntrl),sizeof(int)*20);
@@ -91,8 +91,8 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
}
if(header.icntrl[19]!=0) { // CHARMM format
- skip_flag=(header.icntrl[10]!=0);
- if(skip_flag) {
+ ucell_flag=(header.icntrl[10]!=0);
+ if(ucell_flag) {
LOG_VERBOSE("LoadCHARMMTraj: using CHARMM format with per-frame header");
} else {
LOG_VERBOSE("LoadCHARMMTraj: using CHARMM format");
@@ -128,10 +128,10 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
}
-size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
+size_t calc_frame_size(bool ucell_flag, bool gap_flag, size_t num_atoms)
{
size_t frame_size=0;
- if (skip_flag) {
+ if (ucell_flag) {
frame_size+=14*sizeof(int);
}
if (gap_flag) {
@@ -143,33 +143,12 @@ size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
bool read_frame(std::istream& istream, const DCDHeader& header,
- size_t frame_size, bool skip_flag, bool gap_flag,
+ size_t frame_size, bool ucell_flag, bool gap_flag,
bool swap_flag, std::vector<float>& xlist,
std::vector<geom::Vec3>& frame,
uint frame_num,geom::Vec3& cell_size,geom::Vec3& cell_angles)
{
char dummy[4];
- //if(skip_flag) istream.seekg(14*4,std::ios_base::cur);
- if(skip_flag){
- istream.read(dummy,sizeof(dummy));
- double x,y,z,a,b,c;
- istream.read(reinterpret_cast<char*>(&x),sizeof(double));
- istream.read(reinterpret_cast<char*>(&a),sizeof(double));
- istream.read(reinterpret_cast<char*>(&y),sizeof(double));
- istream.read(reinterpret_cast<char*>(&b),sizeof(double));
- istream.read(reinterpret_cast<char*>(&c),sizeof(double));
- istream.read(reinterpret_cast<char*>(&z),sizeof(double));
- istream.read(dummy,sizeof(dummy));
- cell_size[0]=x;
- cell_size[1]=y;
- cell_size[2]=z;
- cell_angles[0]=acos(a);
- cell_angles[1]=acos(b);
- cell_angles[2]=acos(c);
- if(a!=0.0||b!=0.0||c!=0.0){
- LOG_ERROR("LoadCHARMMTraj: periodic cell not parallelepipedic, cell angles might be wrong, handle carefully")
- }
- }
// read each frame
if(!istream) {
@@ -178,6 +157,21 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
<< frame_num << "Nothing left to be read");
return false;
}
+
+ if(ucell_flag){
+ istream.read(dummy,sizeof(dummy));
+ double tmp[6];
+ istream.read(reinterpret_cast<char*>(tmp),sizeof(double)*6);
+ // a,alpha,b,beta,gamma,c (don't ask)
+ cell_size[0]=static_cast<Real>(tmp[0]);
+ cell_size[1]=static_cast<Real>(tmp[2]);
+ cell_size[2]=static_cast<Real>(tmp[5]);
+ cell_angles[0]=static_cast<Real>(acos(tmp[1]));
+ cell_angles[1]=static_cast<Real>(acos(tmp[3]));
+ cell_angles[2]=static_cast<Real>(acos(tmp[4]));
+ istream.read(dummy,sizeof(dummy));
+ }
+
// x coord
if(gap_flag) istream.read(dummy,sizeof(dummy));
istream.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
@@ -224,12 +218,12 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
}
bool read_frame(std::istream& istream, const DCDHeader& header,
- size_t frame_size, bool skip_flag, bool gap_flag,
+ size_t frame_size, bool ucell_flag, bool gap_flag,
bool swap_flag, std::vector<float>& xlist,
std::vector<geom::Vec3>& frame,uint frame_num)
{
geom::Vec3 cell_size=geom::Vec3(),cell_angles=geom::Vec3();
- return read_frame(istream,header, frame_size,skip_flag, gap_flag,
+ return read_frame(istream,header, frame_size,ucell_flag, gap_flag,
swap_flag, xlist,frame,frame_num, cell_size, cell_angles);
}
@@ -247,8 +241,8 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
Profile profile_load("LoadCHARMMTraj");
DCDHeader header;
- bool swap_flag=false, skip_flag=false, gap_flag=false;
- read_dcd_header(istream, header, swap_flag, skip_flag, gap_flag);
+ bool swap_flag=false, ucell_flag=false, gap_flag=false;
+ read_dcd_header(istream, header, swap_flag, ucell_flag, gap_flag);
if(alist.size() != static_cast<size_t>(header.t_atom_count)) {
LOG_ERROR("LoadCHARMMTraj: atom count missmatch: " << alist.size()
<< " in coordinate file, " << header.t_atom_count
@@ -260,18 +254,18 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
std::vector<geom::Vec3> clist(header.t_atom_count);
std::vector<float> xlist(header.t_atom_count);
geom::Vec3 cell_size, cell_angles;
- size_t frame_size=calc_frame_size(skip_flag, gap_flag, xlist.size());
+ size_t frame_size=calc_frame_size(ucell_flag, gap_flag, xlist.size());
int i=0;
for(;i<header.num;i+=stride) {
- std::cout << i << " " << header.num << std::endl;
- if (!read_frame(istream, header, frame_size, skip_flag, gap_flag,
+ if (!read_frame(istream, header, frame_size, ucell_flag, gap_flag,
swap_flag, xlist, clist, i,cell_size,cell_angles)) {
break;
}
- if(skip_flag) {
+ if(ucell_flag) {
cg.AddFrame(clist,cell_size,cell_angles);
+ } else {
+ cg.AddFrame(clist);
}
- else cg.AddFrame(clist);
// skip frames (defined by stride)
if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur);
@@ -321,7 +315,7 @@ private:
void FetchFrame(uint frame);
String filename_;
DCDHeader header_;
- bool skip_flag_;
+ bool ucell_flag_;
bool swap_flag_;
bool gap_flag_;
std::ifstream stream_;
@@ -337,17 +331,17 @@ private:
void DCDCoordSource::FetchFrame(uint frame)
{
if (!loaded_) {
- read_dcd_header(stream_, header_, swap_flag_, skip_flag_, gap_flag_);
+ read_dcd_header(stream_, header_, swap_flag_, ucell_flag_, gap_flag_);
frame_start_=stream_.tellg();
loaded_=true;
frame_count_=header_.num/stride_;
}
- size_t frame_size=calc_frame_size(skip_flag_, gap_flag_,
+ size_t frame_size=calc_frame_size(ucell_flag_, gap_flag_,
header_.t_atom_count);
size_t pos=frame_start_+frame_size*frame*stride_;
stream_.seekg(pos,std::ios_base::beg);
std::vector<float> xlist(header_.t_atom_count);
- if (!read_frame(stream_, header_, frame_size, skip_flag_, gap_flag_,
+ if (!read_frame(stream_, header_, frame_size, ucell_flag_, gap_flag_,
swap_flag_, xlist, *frame_.get(), frame)) {
}
}
@@ -379,48 +373,65 @@ void write_dcd_hdr(std::ofstream& out,
const mol::CoordGroupHandle& coord_group,
unsigned int stepsize)
{
+ // size of first header block in bytes
int32_t magic_number=84;
- char crd[]={'C', 'O', 'R', 'D'};
out.write(reinterpret_cast<char*>(&magic_number), 4);
+
+ // magic string
+ char crd[]={'C', 'O', 'R', 'D'};
out.write(crd, 4);
+
+ // icntrl[0], NSET, number of frames
int32_t num_frames=coord_group.GetFrameCount()/stepsize;
out.write(reinterpret_cast<char*>(&num_frames), 4);
+
+ // icntrl[1], ISTART, starting timestep
int32_t zero=0;
- // write zero for istart
out.write(reinterpret_cast<char*>(&zero), 4);
+
+ // icntrl[2], NSAVC, timesteps between DCD saves
int32_t one=1;
- // number of timesteps between dcd saves.
out.write(reinterpret_cast<char*>(&one), 4);
- // write spacer of 5 blank integers
- for (int i=0; i<5; ++i) {
+
+ // icntrl[3] to icntrl[7], unused
+ for (int i=3; i<=7; ++i) {
out.write(reinterpret_cast<char*>(&zero), 4);
}
- // write number of fixed atoms
+ // icntrl[8], NAMNF, number of fixed atoms
out.write(reinterpret_cast<char*>(&zero), 4);
+
+ // icntrl[9], DELTA, timestep as float for CHARMM format
float delta=1.0;
out.write(reinterpret_cast<char*>(&delta), 4);
- // write spacer of 10 blank integers
- for (int i=0; i <10; ++i) {
+
+ // icntrl[10], CHARMM format: ucell per frame
+ out.write(reinterpret_cast<char*>(&one), 4);
+
+ // icntrl[11] to icntrl[18], unused
+ for (int i=11; i<=18; ++i) {
out.write(reinterpret_cast<char*>(&zero), 4);
}
+
+ // icntrl[19], charmm version
+ int32_t charmm_version=24;
+ out.write(reinterpret_cast<char*>(&charmm_version), 4);
+ // bracket first header block
out.write(reinterpret_cast<char*>(&magic_number), 4);
- // we don't write any titles for now. This means that the block has only to
- // accomodate one int, harr, harr, harr.
- int32_t title_block_size=4;
- out.write(reinterpret_cast<char*>(&title_block_size), 4);
- out.write(reinterpret_cast<char*>(&zero), 4);
+
+ // no titles in title block
int32_t four=4;
- // again block size for titles?
out.write(reinterpret_cast<char*>(&four), 4);
- // has to be 4
+ out.write(reinterpret_cast<char*>(&zero), 4);
+ out.write(reinterpret_cast<char*>(&four), 4);
+
+ // atom count block
out.write(reinterpret_cast<char*>(&four), 4);
int32_t atom_count=coord_group.GetAtomCount();
out.write(reinterpret_cast<char*>(&atom_count), 4);
out.write(reinterpret_cast<char*>(&four), 4);
}
-
-}
+} // anon ns
void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group,
const String& pdb_filename, const String& dcd_filename,
@@ -442,10 +453,26 @@ void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group,
int32_t out_n=atom_count*4;
for (int i=0; i<frame_count; i+=stepsize) {
mol::CoordFramePtr frame=coord_group.GetFrame(i);
+
+ // ucell
+ int32_t bsize=48; // ucell block size, 6 doubles
+ geom::Vec3 csize=frame->GetCellSize();
+ geom::Vec3 cangles=frame->GetCellAngles();
+ // a,alpha,b,beta,gamma,c (don't ask)
+ double ucell[]={csize[0],
+ cos(cangles[0]),
+ csize[1],
+ cos(cangles[1]),
+ cos(cangles[2]),
+ csize[2]};
+ out.write(reinterpret_cast<char*>(&bsize),4);
+ out.write(reinterpret_cast<char*>(ucell),bsize);
+ out.write(reinterpret_cast<char*>(&bsize),4);
+
+
int k=0;
for (mol::CoordFrame::iterator j=frame->begin(),
e=frame->end(); j!=e; ++j, ++k) {
- //geom::Vec3 v=*j;
x[k]=float((*j)[0]);
y[k]=float((*j)[1]);
z[k]=float((*j)[2]);
diff --git a/modules/io/tests/CMakeLists.txt b/modules/io/tests/CMakeLists.txt
index 98f075bd9785c1c051f642e3be91dd2e125cc62f..9ac4ba1342db1995749b5cce02a1c89fa278b536 100644
--- a/modules/io/tests/CMakeLists.txt
+++ b/modules/io/tests/CMakeLists.txt
@@ -4,6 +4,7 @@ set(OST_IO_UNIT_TESTS
test_clustal.cc
test_io_pdb.cc
test_io_crd.cc
+ test_io_dcd.cc
test_io_sdf.cc
test_pir.cc
test_iomanager.cc
diff --git a/modules/io/tests/test_io_dcd.cc b/modules/io/tests/test_io_dcd.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5bfd0161f354155527dff2d79e7b24a02fb9ca8e
--- /dev/null
+++ b/modules/io/tests/test_io_dcd.cc
@@ -0,0 +1,88 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+
+#include <ost/io/mol/dcd_io.hh>
+#include <ost/mol/entity_handle.hh>
+#include <ost/mol/residue_handle.hh>
+#include <ost/mol/chain_handle.hh>
+#include <ost/mol/atom_handle.hh>
+#include <ost/mol/xcs_editor.hh>
+#include <ost/mol/coord_group.hh>
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+using namespace ost;
+using namespace ost::io;
+
+#include <boost/random.hpp>
+namespace {
+ boost::mt19937 RandomGenerator(time(NULL));
+ boost::uniform_01<boost::mt19937> UniformRandom(RandomGenerator);
+}
+
+BOOST_AUTO_TEST_SUITE( io )
+
+BOOST_AUTO_TEST_CASE(test_io_dcd_charmm_frames)
+{
+ mol::EntityHandle eh=mol::CreateEntity();
+ mol::XCSEditor ed=eh.EditXCS();
+ mol::ChainHandle chain=ed.InsertChain("A");
+ mol::ResidueHandle res=ed.AppendResidue(chain,mol::ResidueKey("UNK"));
+
+ static unsigned int natoms=13;
+
+ // create atoms and frame coords at the same time
+ mol::AtomHandleList atoms(natoms);
+ geom::Vec3List atom_pos(natoms);
+ std::ostringstream aname;
+ geom::Vec3 apos1,apos2;
+ for(size_t i=0;i<natoms;++i) {
+ for(size_t k=0;k<3;++k) {
+ apos1[k]=UniformRandom();
+ apos2[k]=UniformRandom();
+ }
+ aname.str("");
+ aname << "X" << i;
+ mol::AtomHandle atom=ed.InsertAtom(res,aname.str(),apos1);
+ atom_pos[i]=apos2;
+ atoms[i]=atom;
+ }
+ geom::Vec3 cell_size(UniformRandom(),UniformRandom(),UniformRandom());
+ geom::Vec3 cell_angles(M_PI*UniformRandom(),M_PI*UniformRandom(),M_PI*UniformRandom());
+
+ mol::CoordGroupHandle cg=mol::CreateCoordGroup(atoms);
+ cg.AddFrame(atom_pos,cell_size,cell_angles);
+
+ SaveCHARMMTraj(cg,"test_io_dcd_out.pdb","test_io_dcd_out.dcd");
+ mol::CoordGroupHandle cg2=LoadCHARMMTraj(eh,"test_io_dcd_out.dcd");
+ BOOST_CHECK_EQUAL(cg2.GetAtomCount(),natoms);
+ BOOST_CHECK_EQUAL(cg2.GetFrameCount(),1);
+
+ mol::CoordFramePtr cf2 = cg2.GetFrame(0);
+ BOOST_CHECK(geom::Distance(cf2->GetCellSize(),cell_size)<1e-6);
+ BOOST_CHECK(geom::Distance(cf2->GetCellAngles(),cell_angles)<1e-6);
+
+ geom::Vec3List atom_pos2=cg2.GetFramePositions(0);
+
+ for(size_t i=0;i<natoms;++i) {
+ BOOST_CHECK(geom::Distance(atom_pos[i],atom_pos2[i])<1e-6);
+ }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index 1f417064dd8cb2df0ee01956b1f0d0352c6c281f..e3e269e7d7054817d18a3ae642c0e3065846186c 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/modules/mol/base/pymod/export_coord_frame.cc
@@ -16,6 +16,11 @@
// along with this library; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
+
+/*
+ Authors: Marco Biasini, Niklaus Johner, Ansgar Philippsen
+*/
+
#include <boost/python.hpp>
using namespace boost::python;
@@ -26,28 +31,44 @@ using namespace boost::python;
using namespace ost;
using namespace ost::mol;
-geom::Vec3 (CoordFrame::*GetCellSize)() = &CoordFrame::GetCellSize;
-geom::Vec3 (CoordFrame::*GetCellAngles)() = &CoordFrame::GetCellAngles;
-geom::Vec3 (CoordFrame::*get_atom_pos)(const AtomHandle& atom) = &CoordFrame::GetAtomPos;
-Real (CoordFrame::*get_dist_atom)(const AtomHandle& a1, const AtomHandle& a2) = &CoordFrame::GetDistanceBetwAtoms;
-Real (CoordFrame::*get_angle)(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3) = &CoordFrame::GetAngle;
-Real (CoordFrame::*get_dihedral)(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4) = &CoordFrame::GetDihedralAngle;
-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_min_dist_cm_v)(const mol::EntityView& sele1, const mol::EntityView& sele2) = &CoordFrame::GetMinDistBetwCenterOfMassAndView;
-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;
-Real (CoordFrame::*get_alpha)(const mol::EntityView& segment) = &CoordFrame::GetAlphaHelixContent;
-geom::Line3 (CoordFrame::*get_odr_line)(const mol::EntityView& view1) = &CoordFrame::GetODRLine;
-geom::Line3 (CoordFrame::*get_odr_line2)() = &geom::Vec3List::GetODRLine;
-geom::Plane (CoordFrame::*get_odr_plane)(const mol::EntityView& view1) = &CoordFrame::GetODRPlane;
-geom::Line3 (CoordFrame::*fit_cylinder)(const mol::EntityView& view1) = &CoordFrame::FitCylinder;
+geom::Vec3 (CoordFrame::*get_atom_pos)(const AtomHandle&) const = &CoordFrame::GetAtomPos;
+Real (CoordFrame::*get_dist_atom)(const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetDistanceBetwAtoms;
+Real (CoordFrame::*get_angle)(const AtomHandle&, const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetAngle;
+Real (CoordFrame::*get_dihedral)(const AtomHandle&, const AtomHandle&, const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetDihedralAngle;
+geom::Vec3 (CoordFrame::*get_cm)(const mol::EntityView&) const = &CoordFrame::GetCenterOfMassPos;
+Real (CoordFrame::*get_dist_cm)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetDistanceBetwCenterOfMass;
+Real (CoordFrame::*get_min_dist_cm_v)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetMinDistBetwCenterOfMassAndView;
+Real (CoordFrame::*get_rmsd)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetRMSD;
+Real (CoordFrame::*get_min_dist)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetMinDistance;
+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;
+geom::Line3 (CoordFrame::*fit_cylinder)(const mol::EntityView&) const = &CoordFrame::FitCylinder;
+// TODO: move to geom
+geom::Line3 (CoordFrame::*get_odr_line2)() const = &geom::Vec3List::GetODRLine;
+
+CoordFrame create_coord_frame1(const geom::Vec3List& atom_pos)
+{
+ return CreateCoordFrame(atom_pos);
+}
+
+CoordFrame create_coord_frame2(const geom::Vec3List& atom_pos,
+ const geom::Vec3& cell_size,
+ const geom::Vec3& cell_angles)
+{
+ return CreateCoordFrame(atom_pos,cell_size,cell_angles);
+}
void export_CoordFrame()
{
- class_<CoordFrame>("CoordFrame",no_init)
- .def("GetCellSize",GetCellSize)
- .def("GetCellAngles",GetCellAngles)
+ // TODO: add ctors or factory functions that take python sequences
+ class_<CoordFrame>("CoordFrame",init<size_t, optional<geom::Vec3> >())
+ .def("SetCellSize",&CoordFrame::SetCellSize)
+ .def("GetCellSize",&CoordFrame::GetCellSize)
+ .add_property("cell_size",&CoordFrame::GetCellSize,&CoordFrame::SetCellSize)
+ .def("SetCellAngles",&CoordFrame::SetCellAngles)
+ .def("GetCellAngles",&CoordFrame::GetCellAngles)
+ .add_property("cell_angles",&CoordFrame::GetCellAngles,&CoordFrame::SetCellAngles)
.def("GetAtomPos", get_atom_pos)
.def("GetDistanceBetwAtoms", get_dist_atom)
.def("GetAngle", get_angle)
@@ -63,5 +84,6 @@ void export_CoordFrame()
.def("GetAlphaHelixContent",get_alpha)
.def("FitCylinder",fit_cylinder)
;
- def("CreateCoordFrame",CreateCoordFrame);
+ def("CreateCoordFrame",create_coord_frame1);
+ def("CreateCoordFrame",create_coord_frame2);
}
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index 71b9ed9a64d4f21eeed7630e886f47c02ee96d8e..1c30611b3b82a75ffd6660c6980183a0a94e5663 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -18,8 +18,8 @@
//------------------------------------------------------------------------------
/*
- Author: Niklaus Johner
- */
+ Authors: Niklaus Johner, Ansgar Philippsen
+*/
#include <ost/invalid_handle.hh>
#include <ost/integrity_error.hh>
@@ -32,147 +32,76 @@
namespace ost { namespace mol {
- CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos)
- {
- return CoordFrame(atom_pos);
+ namespace {
+ /*
+ Returns a score between 0 and 1 measuring the distance between
+ a pair of angles and a pair of reference angles
+ */
+ Real Eval2AngleDist(Real phi, Real phi0, Real psi, Real psi0, Real delta)
+ {
+ Real d1,d2,d;
+ d1=abs(phi-phi0);
+ d2=abs(psi-psi0);
+ if (d1>M_PI) d1=abs(d1-2.*M_PI);
+ if (d2>M_PI) d1=abs(d2-2.*M_PI);
+ d=(pow(d1,2.)+pow(d2,2.))/pow(delta,2.);
+ return 1./(1+d);
+ }
}
-
- geom::Vec3 CoordFrame::GetAtomPos(const AtomHandle& atom){
- //Returns the position of the atom in this frame
+
+ geom::Vec3 CoordFrame::GetAtomPos(const AtomHandle& atom) const
+ {
return this->GetAtomPos(atom.GetIndex());
}
- geom::Vec3 CoordFrame::GetAtomPos(int i1){
- //Returns the position in this frame of the atom with index i1
+
+ geom::Vec3 CoordFrame::GetAtomPos(int i1) const
+ {
return (*this)[i1];
}
- Real CoordFrame::GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2){
- //Return the distance in this frame between two atoms
+ Real CoordFrame::GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2) const
+ {
return this->GetDistanceBetwAtoms(a1.GetIndex(),a2.GetIndex());
}
- Real CoordFrame::GetDistanceBetwAtoms(int i1, int i2){
- //Return the distance in this frame between the two atoms with indices i1 and i2
+
+ Real CoordFrame::GetDistanceBetwAtoms(int i1, int i2) const
+ {
return geom::Distance((*this)[i1],(*this)[i2]);
}
- Real CoordFrame::GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3){
- //Returns the Angle between the three atoms, a2 being considered as the central atom
- //i.e the angle between the vector (a1.pos-a2.pos) and (a3.pos-a2.pos)
+ Real CoordFrame::GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3) const
+ {
return this->GetAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex());
}
- Real CoordFrame::GetAngle(int i1, int i2, int i3){
- //Returns the angl between the three atoms with indices i1,i2,i3
+
+ Real CoordFrame::GetAngle(int i1, int i2, int i3) const
+ {
return geom::Angle((*this)[i1]-(*this)[i2],(*this)[i3]-(*this)[i2]);
}
Real CoordFrame::GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2,
- const AtomHandle& a3, const AtomHandle& a4){
- //Returns the Dihedral angle between the four atoms a1,a2,a3,a4
+ const AtomHandle& a3, const AtomHandle& a4) const
+ {
return this->GetDihedralAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex(),a4.GetIndex());
}
- Real CoordFrame::GetDihedralAngle(int i1, int i2, int i3, int i4){
- //Returns the Dihedral angle between the four atoms with indices i1,i2,i3,i4
- return geom::DihedralAngle((*this)[i1],(*this)[i2],(*this)[i3],(*this)[i4]);
- }
- void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices){
- // This function returns a vector containing the atom indices of the atoms in an EntityView
- // It is used to accelerate the extraction of information from a trajectory
- AtomViewList atoms=sele.GetAtomList();
- indices.reserve(atoms.size());
- for (AtomViewList::const_iterator i=atoms.begin(),
- e=atoms.end(); i!=e; ++i) {
- indices.push_back(i->GetIndex());
- }
- }
-
- void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca){
- mol::AtomView ca;
- ResidueViewList residues=segment.GetResidueList();
- indices_ca.reserve(residues.size());
- //for (ResidueViewList::const_iterator res=residues.begin()+1,
- // e=residues.end(); res!=e-1; ++res){
- // if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
- //}
- for (ResidueViewList::const_iterator res=residues.begin(),
- e=residues.end(); res!=e; ++res) {
- ca=res->FindAtom("CA");
- CheckHandleValidity(ca);
- indices_ca.push_back(ca.GetIndex());
- }
- }
-
- void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca,
- std::vector<unsigned long>& indices_c, std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n){
- mol::AtomView a,n,c,o;
- ResidueViewList residues=segment.GetResidueList();
- indices_ca.reserve(residues.size());
- indices_n.reserve(residues.size());
- indices_c.reserve(residues.size());
- indices_o.reserve(residues.size());
- for (ResidueViewList::const_iterator res=residues.begin()+1,
- e=residues.end(); res!=e-1; ++res){
- if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
- }
- for (ResidueViewList::const_iterator res=residues.begin(),
- e=residues.end(); res!=e; ++res) {
- a=res->FindAtom("CA");
- CheckHandleValidity(a);
- indices_ca.push_back(a.GetIndex());
- c=res->FindAtom("C");
- CheckHandleValidity(c);
- indices_c.push_back(c.GetIndex());
- o=res->FindAtom("O");
- CheckHandleValidity(o);
- indices_o.push_back(o.GetIndex());
- n=res->FindAtom("N");
- CheckHandleValidity(n);
- indices_n.push_back(n.GetIndex());
- }
+ Real CoordFrame::GetDihedralAngle(int i1, int i2, int i3, int i4) const
+ {
+ return geom::DihedralAngle((*this)[i1],(*this)[i2],(*this)[i3],(*this)[i4]);
}
- void GetMasses(const EntityView& sele, std::vector<Real>& masses){
- // This function returns a vector containing the atom masses of the atoms in an EntityView
- // It is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
- Real mass_tot=0.0;
- AtomViewList atoms=sele.GetAtomList();
- masses.reserve(atoms.size());
- for (AtomViewList::const_iterator i=atoms.begin(),
- e=atoms.end(); i!=e; ++i) {
- masses.push_back(i->GetMass());
- mass_tot=mass_tot+i->GetMass();
- }
- for (std::vector<Real>::iterator
- j=masses.begin(), e2=masses.end(); j!=e2; ++j) {
- (*j)/=mass_tot;
- }
- }
-
- void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
- GetIndices(sele, indices);
- GetMasses(sele, masses);
- }
-
- void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos){
- //Returns the positions of all the atoms in the EntityView
- ref_pos.reserve(sele.GetAtomCount());
- for (mol::AtomViewIter a=sele.AtomsBegin(),e=sele.AtomsEnd(); a!=e; ++a) {
- ref_pos.push_back((*a).GetPos());
- }
- }
-
- geom::Vec3 CoordFrame::GetCenterOfMassPos(const EntityView& sele){
- //Returns the position of the centor of mass of the atoms in the EntityView
+ geom::Vec3 CoordFrame::GetCenterOfMassPos(const EntityView& sele) const
+ {
std::vector<unsigned long> indices;
std::vector<Real> masses;
GetIndicesAndMasses(sele,indices,masses);
return this->GetCenterOfMassPos(indices,masses);
}
- geom::Vec3 CoordFrame::GetCenterOfMassPos(std::vector<unsigned long>& indices,std::vector<Real>& masses){
- //Returns the position of the centor of mass of the atoms from which the indices and masses are passed
- //as vectors in argument
+ geom::Vec3 CoordFrame::GetCenterOfMassPos(std::vector<unsigned long>& indices,
+ std::vector<Real>& masses) const
+ {
geom::Vec3 v;
for (unsigned int i=0,e=indices.size();i!=e; i++) {
v+=masses[i]*(*this)[indices[i]];
@@ -180,8 +109,8 @@ namespace ost { namespace mol {
return v;
}
- Real CoordFrame::GetDistanceBetwCenterOfMass(const EntityView& sele1,const EntityView& sele2){
- //Returns the distance between the centers of mass of the two EntityViews
+ Real CoordFrame::GetDistanceBetwCenterOfMass(const EntityView& sele1,const EntityView& sele2) const
+ {
std::vector<unsigned long> indices1,indices2;
std::vector<Real> masses1,masses2;
GetIndicesAndMasses(sele1,indices1,masses1);
@@ -189,18 +118,19 @@ namespace ost { namespace mol {
return this->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2);
}
- Real CoordFrame::GetDistanceBetwCenterOfMass(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,
- std::vector<unsigned long>& indices2,std::vector<Real>& masses2){
- //Returns the distance between the centers of mass of the two groups of atoms from which the
- //indices and masses are given as vectors as argument
+ Real CoordFrame::GetDistanceBetwCenterOfMass(std::vector<unsigned long>& indices1,
+ std::vector<Real>& masses1,
+ std::vector<unsigned long>& indices2,
+ std::vector<Real>& masses2) const
+ {
geom::Vec3 v1=this->GetCenterOfMassPos(indices1, masses1);
geom::Vec3 v2=this->GetCenterOfMassPos(indices2, masses2);
return geom::Distance(v1,v2);
}
- Real CoordFrame::GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele){
- //Returns the RMSD between the positions of the atoms whose indices are given in indices_sele and the positions
- //given in ref_pos
+ Real CoordFrame::GetRMSD(const std::vector<geom::Vec3>& ref_pos,
+ const std::vector<unsigned long>& indices_sele) const
+ {
Real rmsd=0.0,val;
for (unsigned int i1=0; i1!=indices_sele.size(); ++i1) {
geom::Vec3 av_sele=(*this)[indices_sele[i1]];
@@ -211,9 +141,8 @@ namespace ost { namespace mol {
return pow(rmsd/indices_sele.size(),0.5);
}
- Real CoordFrame::GetRMSD(const EntityView& reference_view,const EntityView& sele_view){
- //Return the rmsd between two EntityViews. The reference positions are taken directly from the reference_view
- //whereas they are taken from this CoordFrame for the sele_view
+ Real CoordFrame::GetRMSD(const EntityView& reference_view,const EntityView& sele_view) const
+ {
int count_ref=reference_view.GetAtomCount();
int count_sele=sele_view.GetAtomCount();
if (count_ref!=count_sele){
@@ -226,8 +155,9 @@ namespace ost { namespace mol {
return this->GetRMSD(ref_pos,indices_sele);
}
- 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
+ Real CoordFrame::GetMinDistance(std::vector<unsigned long>& index_list1,
+ std::vector<unsigned long>& index_list2) const
+ {
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]);
@@ -237,18 +167,19 @@ namespace ost { namespace mol {
}
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)
+
+ Real CoordFrame::GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2) const
+ {
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)
+ Real CoordFrame::GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm,
+ std::vector<Real>& masses_cm,
+ std::vector<unsigned long>& indices_atoms) const
+ {
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) {
@@ -256,8 +187,10 @@ namespace ost { namespace mol {
}
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)
+
+ Real CoordFrame::GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm,
+ const mol::EntityView& view_atoms) const
+ {
std::vector<unsigned long> indices_cm,indices_atoms;
std::vector<Real> masses_cm;
GetIndices(view_atoms,indices_atoms);
@@ -265,8 +198,8 @@ namespace ost { namespace mol {
return this->GetMinDistBetwCenterOfMassAndView(indices_cm,masses_cm,indices_atoms);
}
- geom::Line3 CoordFrame::GetODRLine(std::vector<unsigned long>& indices_ca){
- //Returns the best fit line to atoms with indices in indices_ca
+ geom::Line3 CoordFrame::GetODRLine(std::vector<unsigned long>& indices_ca) const
+ {
geom::Vec3List atoms_pos_list;
atoms_pos_list.reserve(indices_ca.size());
for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
@@ -275,15 +208,14 @@ namespace ost { namespace mol {
return atoms_pos_list.GetODRLine();
}
- geom::Line3 CoordFrame::GetODRLine(const mol::EntityView& view1){
- //Returns the best fit line to atoms in the EntityView view1
+ geom::Line3 CoordFrame::GetODRLine(const mol::EntityView& view1) const {
std::vector<unsigned long> indices;
GetIndices(view1,indices);
return this->GetODRLine(indices);
}
- geom::Plane CoordFrame::GetODRPlane(std::vector<unsigned long>& indices_ca){
- //Returns the normal to the best fit plane to atoms with indices in indices_ca
+ geom::Plane CoordFrame::GetODRPlane(std::vector<unsigned long>& indices_ca) const
+ {
geom::Vec3List atoms_pos_list;
atoms_pos_list.reserve(indices_ca.size());
for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
@@ -292,17 +224,14 @@ namespace ost { namespace mol {
return atoms_pos_list.GetODRPlane();
}
- geom::Plane CoordFrame::GetODRPlane(const mol::EntityView& view1){
- //Returns the best fit plane to atoms in the EntityView view1
+ geom::Plane CoordFrame::GetODRPlane(const mol::EntityView& view1) const
+ {
std::vector<unsigned long> indices;
GetIndices(view1,indices);
return this->GetODRPlane(indices);
}
- 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::Line3 CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca) const {
geom::Vec3List atoms_pos_list;
int n_atoms=indices_ca.size();
atoms_pos_list.reserve(n_atoms);
@@ -326,99 +255,174 @@ namespace ost { namespace mol {
return atoms_pos_list.FitCylinder(initial_axis,center);
}
- geom::Line3 CoordFrame::FitCylinder(const mol::EntityView& view1){
- // Return a lin which is the axis of the best fit cylinder to the CA atoms
- // of the EntityView.
- // It is assumed that we fit an alpha-helix and that the CA atoms are oredered sequentially
+ geom::Line3 CoordFrame::FitCylinder(const mol::EntityView& view1) const {
CheckHandleValidity(view1);
std::vector<unsigned long> indices_ca;
GetCaIndices(view1, indices_ca);
return CoordFrame::FitCylinder(indices_ca);
}
- Real CoordFrame::GetAlphaHelixContent(const mol::EntityView& segment){
- //Returns the percentage of residues in the EntityView (segment) that are in an alpha-helix
- //Each residue is assigned as being in an alpha-helix or not, based on and backbone torsion angles
- //and the hydrogen bonding (O-N) (this term is used only if the segment contains at least 8 residues
- //The entity view should be a single segment with no gaps and no missing N,CA,C,O backbone atoms
- //the first and last residues will not be asessed for helicity
+ Real CoordFrame::GetAlphaHelixContent(const mol::EntityView& segment) const
+ {
CheckHandleValidity(segment);
std::vector<unsigned long> indices_c, indices_n, indices_ca, indices_o;
GetCaCONIndices(segment, indices_ca, indices_c, indices_o, indices_n);
return CoordFrame::GetAlphaHelixContent(indices_ca,indices_c,indices_o,indices_n);
}
- Real Eval2AngleDist(Real phi, Real phi0, Real psi, Real psi0, Real delta){
- //Returns a score between 0 and 1 measuring the distance between
- //a pair of angles and a pair of reference angles
- Real d1,d2,d;
- d1=abs(phi-phi0);
- d2=abs(psi-psi0);
- if (d1>M_PI) d1=abs(d1-2.*M_PI);
- if (d2>M_PI) d1=abs(d2-2.*M_PI);
- d=(pow(d1,2.)+pow(d2,2.))/pow(delta,2.);
- return 1./(1+d);
+ Real CoordFrame::GetAlphaHelixContent(std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_c,
+ std::vector<unsigned long>& indices_o,
+ std::vector<unsigned long>& indices_n) const
+ {
+ //See CoordFrame::GetAlphaHelixContent(const mol::EntityView) above.
+ unsigned long n_atoms=indices_ca.size();
+ geom::Vec3 c_previous,n,ca,c,n_next;
+ Real phi,psi,phi_0=-1.2,psi_0=-0.785,delta=0.8,d_0=3.3;
+ unsigned long n_helical_res=0;
+ std::vector<Real> score,dist,score2,dist2;
+ score.reserve(n_atoms-2);
+ score2.reserve(n_atoms-2);
+ dist.reserve(n_atoms-2);
+ dist2.reserve(n_atoms-2);
+ if (n_atoms!=indices_n.size()||n_atoms!=indices_c.size()||n_atoms!=indices_o.size()){
+ throw std::runtime_error("not same numbers of CA, C, O and N atoms in the selection");
+ }
+ if (n_atoms<=5){
+ throw std::runtime_error("At least five residues are needed to calulate an alpha helix similarity");
+ }
+ c=(*this)[indices_c[0]];
+ n_next=(*this)[indices_n[1]];
+ for (unsigned long i=1; i!=n_atoms-1; ++i) {
+ c_previous=c;
+ n=n_next;
+ ca=(*this)[indices_ca[i]];
+ c=(*this)[indices_c[i]];
+ n_next=(*this)[indices_n[i+1]];
+ phi=geom::DihedralAngle(c_previous,n,ca,c);
+ psi=geom::DihedralAngle(n,ca,c,n_next);
+ score.push_back(Eval2AngleDist(phi,phi_0,psi,psi_0,delta));
+ }
+ score2[0]=pow(score[0]*score[1],3./2.);
+ score2[n_atoms-3]=pow(score[n_atoms-3]*score[n_atoms-4],3./2.);
+ for (unsigned long i=1; i!=n_atoms-3; ++i) {
+ score2[i]=score[i-1]*score[i]*score[i+1];
+ }
+ if (n_atoms>=8){
+ for (unsigned long i=1; i!=n_atoms-1; ++i){
+ if (i<n_atoms-4){
+ dist.push_back(geom::Distance((*this)[indices_o[i]],(*this)[indices_n[i+4]]));
+ if (i>=5) {
+ dist2.push_back(std::min(dist[i-1],dist[i-5]));
+ } else {
+ dist2.push_back(dist[i-1]);
+ }
+ } else {
+ dist2.push_back(geom::Distance((*this)[indices_n[i]],(*this)[indices_o[i-4]]));
+ }
+ }
+ for (unsigned long i=0; i!=n_atoms-2; ++i){
+ if (dist2[i]<=d_0 || score2[i]>=0.3) {
+ n_helical_res+=1;
+ }
+ }
+ } else {
+ for (unsigned long i=0; i!=n_atoms-2; ++i){
+ if (score2[i]>=0.3) {
+ n_helical_res+=1;
+ }
+ }
+ }
+ return Real(n_helical_res)/Real(n_atoms-2);
+ }
+
+ CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos, const geom::Vec3& cs, const geom::Vec3& ca)
+ {
+ return CoordFrame(atom_pos,cs,ca);
}
- Real CoordFrame::GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
- std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n){
- //See CoordFrame::GetAlphaHelixContent(const mol::EntityView) above.
- unsigned long n_atoms=indices_ca.size();
- geom::Vec3 c_previous,n,ca,c,n_next;
- Real phi,psi,phi_0=-1.2,psi_0=-0.785,delta=0.8,d_0=3.3;
- unsigned long n_helical_res=0;
- std::vector<Real> score,dist,score2,dist2;
- score.reserve(n_atoms-2);
- score2.reserve(n_atoms-2);
- dist.reserve(n_atoms-2);
- dist2.reserve(n_atoms-2);
- if (n_atoms!=indices_n.size()||n_atoms!=indices_c.size()||n_atoms!=indices_o.size()){
- throw std::runtime_error("not same numbers of CA, C, O and N atoms in the selection");
- }
- if (n_atoms<=5){
- throw std::runtime_error("At least five residues are needed to calulate an alpha helix similarity");
- }
- c=(*this)[indices_c[0]];
- n_next=(*this)[indices_n[1]];
- for (unsigned long i=1; i!=n_atoms-1; ++i){
- c_previous=c;
- n=n_next;
- ca=(*this)[indices_ca[i]];
- c=(*this)[indices_c[i]];
- n_next=(*this)[indices_n[i+1]];
- phi=geom::DihedralAngle(c_previous,n,ca,c);
- psi=geom::DihedralAngle(n,ca,c,n_next);
- score.push_back(Eval2AngleDist(phi,phi_0,psi,psi_0,delta));
- }
- score2[0]=pow(score[0]*score[1],3./2.);
- score2[n_atoms-3]=pow(score[n_atoms-3]*score[n_atoms-4],3./2.);
- for (unsigned long i=1; i!=n_atoms-3; ++i){
- score2[i]=score[i-1]*score[i]*score[i+1];
- }
- if (n_atoms>=8){
- for (unsigned long i=1; i!=n_atoms-1; ++i){
- if (i<n_atoms-4){
- dist.push_back(geom::Distance((*this)[indices_o[i]],(*this)[indices_n[i+4]]));
- if (i>=5) {
- dist2.push_back(std::min(dist[i-1],dist[i-5]));
- }
- else dist2.push_back(dist[i-1]);
+ void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices)
+ {
+ AtomViewList atoms=sele.GetAtomList();
+ indices.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ indices.push_back(i->GetIndex());
+ }
+ }
+
+ void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca){
+ mol::AtomView ca;
+ ResidueViewList residues=segment.GetResidueList();
+ indices_ca.reserve(residues.size());
+ //for (ResidueViewList::const_iterator res=residues.begin()+1,
+ // e=residues.end(); res!=e-1; ++res){
+ // if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
+ //}
+ for (ResidueViewList::const_iterator res=residues.begin(),
+ e=residues.end(); res!=e; ++res) {
+ ca=res->FindAtom("CA");
+ CheckHandleValidity(ca);
+ indices_ca.push_back(ca.GetIndex());
}
- else dist2.push_back(geom::Distance((*this)[indices_n[i]],(*this)[indices_o[i-4]]));
+ }
+
+ void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_c, std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n){
+ mol::AtomView a,n,c,o;
+ ResidueViewList residues=segment.GetResidueList();
+ indices_ca.reserve(residues.size());
+ indices_n.reserve(residues.size());
+ indices_c.reserve(residues.size());
+ indices_o.reserve(residues.size());
+ for (ResidueViewList::const_iterator res=residues.begin()+1,
+ e=residues.end(); res!=e-1; ++res){
+ if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
}
- for (unsigned long i=0; i!=n_atoms-2; ++i){
- if (dist2[i]<=d_0 || score2[i]>=0.3)n_helical_res+=1;
+ for (ResidueViewList::const_iterator res=residues.begin(),
+ e=residues.end(); res!=e; ++res) {
+ a=res->FindAtom("CA");
+ CheckHandleValidity(a);
+ indices_ca.push_back(a.GetIndex());
+ c=res->FindAtom("C");
+ CheckHandleValidity(c);
+ indices_c.push_back(c.GetIndex());
+ o=res->FindAtom("O");
+ CheckHandleValidity(o);
+ indices_o.push_back(o.GetIndex());
+ n=res->FindAtom("N");
+ CheckHandleValidity(n);
+ indices_n.push_back(n.GetIndex());
}
}
- else {
- for (unsigned long i=0; i!=n_atoms-2; ++i){
- if (score2[i]>=0.3)n_helical_res+=1;
+
+ void GetMasses(const EntityView& sele, std::vector<Real>& masses){
+ Real mass_tot=0.0;
+ AtomViewList atoms=sele.GetAtomList();
+ masses.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ masses.push_back(i->GetMass());
+ mass_tot=mass_tot+i->GetMass();
}
+ for (std::vector<Real>::iterator
+ j=masses.begin(), e2=masses.end(); j!=e2; ++j) {
+ (*j)/=mass_tot;
+ }
}
- return Real(n_helical_res)/Real(n_atoms-2);
+
+
+ void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
+ GetIndices(sele, indices);
+ GetMasses(sele, masses);
+ }
+
+ void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos){
+ ref_pos.reserve(sele.GetAtomCount());
+ for (mol::AtomViewIter a=sele.AtomsBegin(),e=sele.AtomsEnd(); a!=e; ++a) {
+ ref_pos.push_back((*a).GetPos());
+ }
}
-
-
}}//ns
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index 98eb0c25c5560dbce0ad43e8450e79b8994b46df..ad76a5dc869d281775ca61fe0ce191fc6a4d8757 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -20,7 +20,7 @@
#define OST_MOL_COORD_FRAME_HH
/*
- Author: Marco Biasini and Niklaus Johner
+ Authors: Marco Biasini, Niklaus Johner, Ansgar Philippsen
*/
#include <boost/shared_ptr.hpp>
#include <ost/geom/geom.hh>
@@ -30,10 +30,22 @@
namespace ost { namespace mol {
+/*!
+ A coordinate frame in a trajectory, containing atom positions and
+ per-frame unit cell size and angles.
+*/
+
+/*
+ TODO:
+ - move algorithmic code to separate functions in mol/alg
+ - clean up mix of views and atom indices methods/functions
+ - use existing UnitCell class
+*/
+
class DLLEXPORT_OST_MOL CoordFrame : public geom::Vec3List{
private:
- geom::Vec3 periodic_cell_sizes_;
- geom::Vec3 periodic_cell_angles_;
+ geom::Vec3 ucell_size_;
+ geom::Vec3 ucell_angles_;
public:
typedef geom::Vec3List base_type;
@@ -42,62 +54,194 @@ public:
CoordFrame(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
CoordFrame(const base_type& rhs) : base_type(rhs) { }
CoordFrame(const std::vector<geom::Vec3>& rhs) : base_type(rhs) { }
- CoordFrame(const std::vector<geom::Vec3>& rhs,const geom::Vec3 box_size, const geom::Vec3 box_angles) : base_type(rhs) {
- this->periodic_cell_sizes_=box_size;
- this->periodic_cell_angles_=box_angles;
+ CoordFrame(const std::vector<geom::Vec3>& rhs,
+ const geom::Vec3 box_size,
+ const geom::Vec3 box_angles) : base_type(rhs) {
+ ucell_size_=box_size;
+ ucell_angles_=box_angles;
}
-
- geom::Vec3 GetCellSize() {
- return this->periodic_cell_sizes_;
+
+ void SetCellSize(const geom::Vec3& s) {
+ ucell_size_=s;
}
- geom::Vec3 GetCellAngles() {
- return this->periodic_cell_angles_;
+
+ geom::Vec3 GetCellSize() const {
+ return this->ucell_size_;
}
- geom::Vec3 GetAtomPos(const AtomHandle& atom);
- geom::Vec3 GetAtomPos(int atom_index);
- Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2);
- Real GetDistanceBetwAtoms(int atom1_index, int atom2_index);
- Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3);
- Real GetAngle(int atom1_index, int atom2_index, int atom3_index);
- Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4);
- Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index);
- geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele);
- geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices, std::vector<Real>& masses);
- Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2);
+
+ void SetCellAngles(const geom::Vec3& a) {
+ ucell_angles_=a;
+ }
+
+ geom::Vec3 GetCellAngles() const {
+ return this->ucell_angles_;
+ }
+
+ geom::Vec3 GetAtomPos(const AtomHandle& atom) const;
+ geom::Vec3 GetAtomPos(int atom_index) const;
+ Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2) const;
+ Real GetDistanceBetwAtoms(int atom1_index, int atom2_index) const;
+
+ /*!
+ Returns the angle between the three atoms in rad
+
+ a2 is considered the central atom, i.e the result is the angle between
+ the vectors (a1.pos-a2.pos) and (a3.pos-a2.pos)
+ */
+ Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3) const;
+
+ /*!
+ Returns the angle between the three atom indices in rad
+
+ atom2 is considered the central atom, i.e the result is the angle between
+ the vectors (atom1.pos-atom2.pos) and (atom3.pos-atom2.pos)
+ */
+ Real GetAngle(int atom1_index, int atom2_index, int atom3_index) const;
+
+ //! Returns the Dihedral angle between the four atoms a1,a2,a3,a4
+ Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4) const;
+
+ //! Returns the Dihedral angle between the four atom indices a1,a2,a3,a4
+ Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index) const;
+
+ //! Returns the position of the centor of mass of the atoms in the EntityView
+ geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele) const;
+
+ //! Calculates the center of mass given a list of indices and provided masses
+ /*!
+ the masses vector contains only the masses for the atoms indexed by indices, not for
+ all atoms in the frame
+ */
+ geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices,
+ std::vector<Real>& masses) const;
+
+ //! Returns the distance between the centers of mass of the two EntityViews
+ Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2) const;
+
+ /*! Returns the distance betweem two custom center of masses
+
+ Calculates the centers of mass of the two groups of atoms from given by index lists and
+ custom masses. See comment regarding masses vector for GetCenterOfMassPos
+ */
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 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);
- geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca);
- geom::Plane GetODRPlane(std::vector<unsigned long>& indices_ca);
- geom::Line3 GetODRLine(const mol::EntityView& view1);
- geom::Plane GetODRPlane(const mol::EntityView& view1);
- geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
- geom::Line3 FitCylinder(const mol::EntityView& view1);
- Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
- std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
- Real GetAlphaHelixContent(const mol::EntityView& segment);
-};
+ std::vector<unsigned long>& indices2, std::vector<Real>& masses2) const;
+
+ /*!
+ Returns the RMSD between the positions of the atoms whose indices are given in indices_sele
+ and the positions given in ref_pos
+ */
+ Real GetRMSD(const std::vector<geom::Vec3>& ref_pos,
+ const std::vector<unsigned long>& indices_sele) const;
+ /*!
+ Returns the rmsd between two EntityViews.
+
+ The reference positions are taken directly from the reference_view,
+ whereas they are taken from this CoordFrame for the sele_view
+ */
+ Real GetRMSD(const mol::EntityView& reference_view,
+ const mol::EntityView& sele_view) const;
+
+ //! Returns the minimal distance between two groups of atoms
+ /*!
+ with indices given by index_list1 and index_list2
+ */
+ Real GetMinDistance(std::vector<unsigned long>& index_list1,
+ std::vector<unsigned long>& index_list2) const;
+
+ //! Returns the minimal distance between the atoms of two views (view1 and view2)
+ Real GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2) const;
+
+ /*!
+ Returns the minimal distance between the center of mass of two groups of atoms
+
+ atoms and their masses are defined by indices and masses for those atoms
+ */
+ Real GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm,
+ std::vector<Real>& masses_cm,
+ std::vector<unsigned long>& indices_atoms) const;
+ //! Returns the minimal distance between the center of mass of two views
+ Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm,
+ const mol::EntityView& view_atoms) const;
+
+ //! Returns the best fit line to atoms with indices in indices_ca
+ geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca) const;
+
+ //! Returns the best fit line to atoms in the EntityView view1
+ geom::Plane GetODRPlane(std::vector<unsigned long>& indices_ca) const;
+
+ //! Returns the normal to the best fit plane to atoms with indices in indices_ca
+ geom::Line3 GetODRLine(const mol::EntityView& view1) const;
+ geom::Plane GetODRPlane(const mol::EntityView& view1) const;
+
+ /*!
+ 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::Line3 FitCylinder(std::vector<unsigned long>& indices_ca) const;
+
+ //! see FitCylinder(std::vector<unsigned long>&)
+ geom::Line3 FitCylinder(const mol::EntityView& view1) const;
+
+ /*!
+ Returns the percentage of residues in the EntityView (segment) that are in an alpha-helix
+ Each residue is assigned as being in an alpha-helix or not, based on and backbone torsion angles
+ and the hydrogen bonding (O-N) (this term is used only if the segment contains at least 8 residues
+ The entity view should be a single segment with no gaps and no missing N,CA,C,O backbone atoms
+ the first and last residues will not be asessed for helicity
+ */
+ Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_c,
+ std::vector<unsigned long>& indices_o,
+ std::vector<unsigned long>& indices_n) const;
+
+ //! see above
+ Real GetAlphaHelixContent(const mol::EntityView& segment) const;
+};
+
+ typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
+ typedef std::vector<CoordFramePtr> CoordFrameList;
+
+ // factory method
+ // create a frame from a Vec3List containing the positions of the atoms, optionally with unit cell
+ // TODO: why is this necessary? the ctor works fine
+ DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos,
+ const geom::Vec3& cell_size=geom::Vec3(),
+ const geom::Vec3& cell_angles=geom::Vec3());
+
+ // these should really be either in the entity view interface or in mol/alg
+
+ /*!
+ This function returns a vector containing the atom indices of the atoms in an EntityView;
+ it is used to accelerate the extraction of information from a trajectory
+ */
void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
+
+ /*!
+ This function returns a vector containing the atom masses of the atoms in an EntityView;
+ it is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
+ */
void GetMasses(const EntityView& sele, std::vector<Real>& masses);
- void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
+
+ //! conveniece for GetIndices and GetMasses in one call
+ void GetIndicesAndMasses(const EntityView& sele,
+ std::vector<unsigned long>& indices,
+ std::vector<Real>& masses);
+
+ //! Writes the positions of all atoms in the EntityView into the provided vec3 list
void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
+
+ //! Writes the indices of all atoms in the EntityView into the provided list
void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca);
- void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
- std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
-typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
-typedef std::vector<CoordFramePtr> CoordFrameList;
+ //! Writes the backbone indices of all residues in the EntityView into the provided list
+ void GetCaCONIndices(const EntityView& segment,
+ std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_c,
+ std::vector<unsigned long>& indices_o,
+ std::vector<unsigned long>& indices_n);
-// factory method
-// create a frame froma Vec3List containing the positions of the atoms
- DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos);
}}