From 1350fae62ac698ff038d1feeb40fe2ded4d9521b Mon Sep 17 00:00:00 2001
From: Ansgar Philippsen <ansgar.philippsen@gmail.com>
Date: Mon, 21 May 2012 10:32:50 -0400
Subject: [PATCH] fixes to dcd io and coord frame
fixed buggy dcd unit cell import
fixed missing dcd unit cell export
added unit cell for dcd io
some code cleanup for dcd export
some code cleanup for CoordFrame and Vec3List
---
modules/geom/src/composite3_op.cc | 71 ---
modules/geom/src/vec3.cc | 78 +++-
modules/geom/src/vec3.hh | 52 ++-
modules/io/src/mol/dcd_io.cc | 145 +++---
modules/io/tests/CMakeLists.txt | 1 +
modules/io/tests/test_io_dcd.cc | 88 ++++
modules/mol/base/pymod/export_coord_frame.cc | 62 ++-
modules/mol/base/src/coord_frame.cc | 444 ++++++++++---------
modules/mol/base/src/coord_frame.hh | 240 ++++++++--
9 files changed, 738 insertions(+), 443 deletions(-)
create mode 100644 modules/io/tests/test_io_dcd.cc
diff --git a/modules/geom/src/composite3_op.cc b/modules/geom/src/composite3_op.cc
index 87f799d51..09fb2c72e 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 ea07f3cca..6e657081f 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 ac3458837..7b4bb5b33 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 15a290279..5acc5bf4b 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 98f075bd9..9ac4ba134 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 000000000..5bfd0161f
--- /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 1f417064d..e3e269e7d 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 71b9ed9a6..1c30611b3 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 98eb0c25c..ad76a5dc8 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);
}}
--
GitLab