diff --git a/modules/geom/pymod/export_vecmat3_op.cc b/modules/geom/pymod/export_vecmat3_op.cc
index 048cdd170b4ca4e45c724db50f21d7a93fdcabd6..b858001faac7bdaf6485c4433e6acb6b0da69cbd 100644
--- a/modules/geom/pymod/export_vecmat3_op.cc
+++ b/modules/geom/pymod/export_vecmat3_op.cc
@@ -38,7 +38,8 @@ Real (*Mat3Comp)(const Mat3& m, unsigned int i, unsigned int j) = &Comp;
Real (*Mat3Minor)(const Mat3& m, unsigned int i, unsigned int j) = &Minor;
Vec3 (*Vec3Min)(const Vec3&, const Vec3&) = &Min;
Vec3 (*Vec3Max)(const Vec3&, const Vec3&) = &Max;
-
+Real (*Vec3Distance2WithPBC)(const Vec3&, const Vec3&, const Vec3&) = &Distance2WithPBC;
+Real (*Vec3DistanceWithPBC)(const Vec3&, const Vec3&, const Vec3&) = &DistanceWithPBC;
void export_VecMat3_op()
{
@@ -64,4 +65,10 @@ void export_VecMat3_op()
def("OrthogonalVector",OrthogonalVector);
def("Min",Vec3Min);
def("Max",Vec3Max);
+ def("Distance2WithPBC",Vec3Distance2WithPBC);
+ def("DistanceWithPBC",Vec3DistanceWithPBC);
+ def("MinDistance",MinDistance);
+ def("MinDistanceWithPBC",MinDistanceWithPBC);
+ def("WrapVec3",WrapVec3);
+ def("WrapVec3List",WrapVec3List);
}
diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc
index 0ce26dde3f1684710265b3c0cc1c79f69ee8bf01..50b5b612e939d66e6d8119ddced416ff87e52f99 100644
--- a/modules/geom/src/vecmat3_op.cc
+++ b/modules/geom/src/vecmat3_op.cc
@@ -192,19 +192,61 @@ Real DihedralAngle(const Vec3& p1, const Vec3& p2, const Vec3& p3,
Dot(r12cross, r23cross));
}
+
Real MinDistance(const Vec3List& l1, const Vec3List& l2)
{
// returns the minimal distance between two sets of points (Vec3List)
if (l1.size()==0 || l2.size()==0){throw std::runtime_error("cannot calculate minimal distance: empty Vec3List");}
- Real min=Distance(*l1.begin(),*l2.begin());
+ Real min=Length2(*l1.begin()-*l2.begin());
+ Real d;
+ for (Vec3List::const_iterator p1=l1.begin(),e1=l1.end(); p1!=e1; p1++) {
+ for (Vec3List::const_iterator p2=l2.begin(),e2=l2.end(); p2!=e2; p2++) {
+ d=Length2(*p1-*p2);
+ if (d<min) min=d;
+ }
+ }
+ return std::sqrt(min);
+}
+
+Real MinDistanceWithPBC(const Vec3List& l1, const Vec3List& l2, Vec3& basis_vec)
+{
+ // returns the minimal distance between two sets of points (Vec3List)
+ // given the periodic boundary condition along x,y,z given in the basis_vec
+ if (l1.size()==0 || l2.size()==0){throw std::runtime_error("cannot calculate minimal distance: empty Vec3List");}
+ Real min=Length2(*l1.begin()-*l2.begin());
Real d;
+ Vec3 v;
+ for (int i=0; i<3; i++) {
+ basis_vec[i]=std::fabs(basis_vec[i]);
+ }
for (Vec3List::const_iterator p1=l1.begin(),e1=l1.end(); p1!=e1; p1++) {
for (Vec3List::const_iterator p2=l2.begin(),e2=l2.end(); p2!=e2; p2++) {
- d=Distance(*p1,*p2);
+ d=Distance2WithPBC(*p1, *p2, basis_vec);
if (d<min) min=d;
}
}
- return min;
+ return std::sqrt(min);
+}
+
+Vec3 WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& basis_vec){
+ Vec3 v;
+ Real r;
+ for (int i=0; i<3; i++) {
+ r=(v1[i]-box_center[i])/basis_vec[i];
+ r=(r > 0.0) ? std::floor(r + 0.5) : std::ceil(r - 0.5);
+ v[i]=v1[i]-basis_vec[i]*r;
+ }
+ return v;
}
+Vec3List WrapVec3List(const Vec3List& vl, const Vec3& box_center,const Vec3& basis_vec){
+ Vec3List vl_out;
+ vl_out.reserve(vl_out.size());
+ for (Vec3List::const_iterator v1=vl.begin(),e=vl.end();v1!=e ; v1++) {
+ vl_out.push_back(WrapVec3(*v1,box_center,basis_vec));
+ }
+ return vl_out;
+}
+
+
} // ns
diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh
index 373b9561b1338fe3907e5e37809b4d71b501911f..5489ea9e8f1c48f669053fafa45f14220b52384f 100644
--- a/modules/geom/src/vecmat3_op.hh
+++ b/modules/geom/src/vecmat3_op.hh
@@ -194,9 +194,34 @@ inline Real Distance(const Vec3& p1, const Vec3& p2)
return Length(p1-p2);
}
+
+//! return the squared distance between two points with periodic boundaries in x,y,z given in basis_vec
+inline Real Distance2WithPBC(const Vec3& v1, const Vec3& v2, const Vec3& basis_vec){
+ Vec3 v;
+ v=v1-v2;
+ for (int i=0; i<3; i++) {
+ if (std::fabs(v[i])>basis_vec[i]/2.){
+ v[i]=std::fabs(v[i])-basis_vec[i]*int(std::fabs(v[i])/basis_vec[i]+0.5);
+ }
+ }
+ return Length2(v);
+}
+//! return the distance between two points with periodic boundaries in x,y,z given in basis_vec
+inline Real DistanceWithPBC(const Vec3& v1, const Vec3& v2, const Vec3& basis_vec){
+ return sqrt(Distance2WithPBC(v1, v2, basis_vec));
+}
//! returns the minimal distance between the points in two Vec3List
Real MinDistance(const Vec3List& l1, const Vec3List& l2);
+//! returns the minimal distance between the points in two Vec3List
+// with periodic boundaries in x,y,z given in basis_vec
+Real MinDistanceWithPBC(const Vec3List& l1, const Vec3List& l2, Vec3& basis_vec);
+
+//!wraps a vector in a box with periodic boundaries
+Vec3 WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& basis_vec);
+//!wraps all the verctors in a Vec3List in a box with periodic boundaries
+Vec3List WrapVec3List(const Vec3List& vl,const Vec3& box_center,const Vec3& basis_vec);
+
} // ns
#endif
diff --git a/modules/io/pymod/__init__.py b/modules/io/pymod/__init__.py
index e22106b51fc4aba65bb38950d0709c631e020621..95ea418b344a21910a925f463151d34af4a8eb22 100644
--- a/modules/io/pymod/__init__.py
+++ b/modules/io/pymod/__init__.py
@@ -194,10 +194,12 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
def SavePDB(models, filename, dialect=None, pqr=False, profile='DEFAULT'):
"""
- Save entity or list of entities to disk. If a list of entities is supplied the
- PDB file will be saved as a multi PDB file. Each of the entities is wrapped
- into a MODEL/ENDMDL pair.
-
+ Save entity or list of entities to disk. If a list of entities is supplied
+ the PDB file will be saved as a multi PDB file. Each of the entities is
+ wrapped into a MODEL/ENDMDL pair.
+
+ If the atom number exceeds 99999, '*****' is used.
+
:param models: The entity or list of entities (handles or views) to be saved
:param filename: The filename
:type filename: string
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index 58c1127a1a68617471bb8af2b6e1544138944188..e5422bda63ed4db572a15749cb692706452ec924 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/modules/io/src/mol/dcd_io.cc
@@ -141,14 +141,22 @@ size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
return frame_size;
}
+
bool read_frame(std::istream& istream, const DCDHeader& header,
size_t frame_size, bool skip_flag, bool gap_flag,
bool swap_flag, std::vector<float>& xlist,
std::vector<geom::Vec3>& frame,
- uint frame_num)
+ 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.seekg(14*4,std::ios_base::cur);
+ if(skip_flag){
+ istream.read(dummy,sizeof(dummy));
+ istream.read(reinterpret_cast<char*>(&cell_size[0]),sizeof(float)*3);
+ istream.read(reinterpret_cast<char*>(&cell_angles[0]),sizeof(float)*3);
+ istream.read(dummy,sizeof(dummy));
+ }
+
// read each frame
if(!istream) {
/* premature EOF */
@@ -201,7 +209,17 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
return true;
}
-
+bool read_frame(std::istream& istream, const DCDHeader& header,
+ size_t frame_size, bool skip_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,
+ swap_flag, xlist,frame,frame_num, cell_size, cell_angles);
+}
+
+
mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted!
const String& trj_fn,
unsigned int stride)
@@ -227,20 +245,23 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
mol::CoordGroupHandle cg=CreateCoordGroup(alist);
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());
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,
- swap_flag, xlist, clist, i)) {
+ swap_flag, xlist, clist, i,cell_size,cell_angles)) {
break;
}
- cg.AddFrame(clist);
+ if(skip_flag) {
+ cg.AddFrame(clist,cell_size,cell_angles);
+ }
+ else cg.AddFrame(clist);
// skip frames (defined by stride)
if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur);
}
-
istream.get();
if(!istream.eof()) {
LOG_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: "
@@ -279,6 +300,7 @@ public:
}
virtual void AddFrame(const std::vector<geom::Vec3>& coords) {}
+ virtual void AddFrame(const std::vector<geom::Vec3>& coords,const geom::Vec3& box_size,const geom::Vec3& box_angles) {}
virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords) {}
private:
diff --git a/modules/io/src/mol/pdb_reader.cc b/modules/io/src/mol/pdb_reader.cc
index 8fa9b53b41572eb827e365b1090cfcb5096f7e02..c52e5cf3151e1fc1e2c59f39616e822a4b7dd90f 100644
--- a/modules/io/src/mol/pdb_reader.cc
+++ b/modules/io/src/mol/pdb_reader.cc
@@ -602,16 +602,6 @@ bool PDBReader::ParseAtomIdent(const StringRef& line, int line_num,
}
}
- std::pair<bool, int> a_num=line.substr(6, 5).ltrim().to_int();
- if (!a_num.first) {
- if (!(profile_.fault_tolerant)) {
- throw IOException(str(format("invalid atom number on line %d") %line_num));
- }
- if (!(charmm_style_)) {
- LOG_WARNING("invalid atom number on line " << line_num);
- }
- }
-
alt_loc=line[16];
res_name=line.substr(17, charmm_style_ ? 4 : 3).trim();
std::pair<bool, int> res_num=line.substr(22, 4).ltrim().to_int();;
diff --git a/modules/io/src/mol/pdb_writer.cc b/modules/io/src/mol/pdb_writer.cc
index e18d96d7fbc551494ca0277a9f05b5fda9b7c226..f3c863594c926c330426c5322f0638715cb676ce 100644
--- a/modules/io/src/mol/pdb_writer.cc
+++ b/modules/io/src/mol/pdb_writer.cc
@@ -80,12 +80,7 @@ void write_atom(std::ostream& ostr, FormattedLine& line,
line( 0, 6)=record_name;
// Avoid writing out atomnumbers larger than 5 digits
if (atomnum > 99999) {
- if (charmm_style) {
- line( 6, 5)=fmt::LPadded("*****");
- } else {
- throw IOException("Atom number is too long for PDB output."
- " At most 5 digits are allowed");
- }
+ line( 6, 5)=fmt::LPadded("*****");
} else {
line( 6, 5)=fmt::LPaddedInt(atomnum);
}
diff --git a/modules/io/tests/test_io_pdb.cc b/modules/io/tests/test_io_pdb.cc
index 48ba5667caf3a5dfe279f5d2852c52693bf3e02e..6c6092198ad106deb8d0cb7135d0f5acce9280f0 100644
--- a/modules/io/tests/test_io_pdb.cc
+++ b/modules/io/tests/test_io_pdb.cc
@@ -457,6 +457,64 @@ BOOST_AUTO_TEST_CASE(write_atom)
" 1.00128.00 C ");
}
+BOOST_AUTO_TEST_CASE(write_atom_100000)
+{
+ char c_names[] = "ABCDEFGHIJ";
+ std::stringstream out;
+ PDBWriter writer(out, IOProfile());
+
+ mol::EntityHandle ent=mol::CreateEntity();
+ mol::XCSEditor edi=ent.EditXCS();
+ mol::ChainHandle ch;
+ mol::ResidueHandle r;
+ mol::AtomHandle a;
+
+ for (unsigned long i = 0; i < 10; ++i) {
+ ch=edi.InsertChain(String(1, c_names[i]));
+ for (unsigned long j = 1; j < 1000; ++j) {
+ r=edi.AppendResidue(ch, "ARG");
+ a=edi.InsertAtom(r,"N", geom::Vec3(26.861, 50.841, 38.803), "N");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"CA", geom::Vec3(27.437, 49.969, 37.786), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"C", geom::Vec3(26.336, 48.959, 37.429), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"O", geom::Vec3(25.745, 48.313, 38.312), "O");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"CB", geom::Vec3(28.653, 49.266, 38.349), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"CG", geom::Vec3(29.870, 50.188, 38.416), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"CD", geom::Vec3(31.033, 49.532, 39.173), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"NE", geom::Vec3(32.318, 50.244, 39.125), "N");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"CZ", geom::Vec3(33.462, 49.750, 39.679), "C");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"NH1", geom::Vec3(33.522, 48.572, 40.308), "N");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ a=edi.InsertAtom(r,"NH2", geom::Vec3(34.610, 50.427, 39.597), "N");
+ a.SetOccupancy(1.0);
+ a.SetBFactor(128.0);
+ }
+ }
+
+ writer.Write(ent);
+ String s=out.str();
+ BOOST_CHECK_EQUAL(s.substr(8099844, 5), "99999");
+ BOOST_CHECK_EQUAL(s.substr(8099925, 5), "*****");
+}
+
BOOST_AUTO_TEST_CASE(write_hetatom)
{
std::stringstream out;
diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index b4f97c5cb8b1ce228e9d1bc5fabec9ab1253cd22..004c475efdf99c095aa693c133f7a2af51e84248 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -3,6 +3,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
export_svd_superpose.cc
export_clash.cc
export_trajectory_analysis.cc
+ export_structure_analysis.cc
)
set(OST_MOL_ALG_PYMOD_MODULES
diff --git a/modules/mol/alg/pymod/export_structure_analysis.cc b/modules/mol/alg/pymod/export_structure_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5a91310949eab28d1ba49dfdaa42763c8beff47e
--- /dev/null
+++ b/modules/mol/alg/pymod/export_structure_analysis.cc
@@ -0,0 +1,33 @@
+//------------------------------------------------------------------------------
+// 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 <boost/python.hpp>
+using namespace boost::python;
+
+#include <ost/mol/alg/structure_analysis.hh>
+
+using namespace ost;
+using namespace ost::mol::alg;
+
+void export_StructureAnalysis()
+{
+ def("GetPosListFromView",&GetPosListFromView, (arg("view")));
+ def("CalculateAverageAgreementWithDensityMap",&CalculateAverageAgreementWithDensityMap,(arg("pos_list"),arg("density_map")));
+ def("CalculateAgreementWithDensityMap",&CalculateAgreementWithDensityMap,(arg("pos_list"),arg("density_map")));
+ def("WrapEntityInPeriodicCell",&WrapEntityInPeriodicCell,(arg("Entity"),arg("cell_center"),arg("nasis_vec")));
+}
\ No newline at end of file
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index fb68ece360a33a3729c11902f1df26a746bd43d8..20a78f90e3ede2920d09b7d1cadc70cfc1fe65df 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -28,6 +28,7 @@ using namespace ost;
void export_svdSuperPose();
void export_TrajectoryAnalysis();
+void export_StructureAnalysis();
void export_Clash();
#if OST_IMG_ENABLED
void export_entity_to_density();
@@ -76,6 +77,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
{
export_svdSuperPose();
export_TrajectoryAnalysis();
+ export_StructureAnalysis();
#if OST_IMG_ENABLED
export_entity_to_density();
#endif
@@ -143,4 +145,4 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("FillStereoChemicalParams",&fill_stereochemical_params_wrapper);
def("IsStandardResidue",&mol::alg::IsStandardResidue);
-}
\ No newline at end of file
+}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 99c1cd6778bbf5781f06f5079a32ad9e06d8016d..050d4822712f8789bfac97938e18c327334f7d92 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -8,6 +8,7 @@ set(OST_MOL_ALG_HEADERS
construct_cbeta.hh
clash_score.hh
trajectory_analysis.hh
+ structure_analysis.hh
)
set(OST_MOL_ALG_SOURCES
@@ -19,6 +20,7 @@ set(OST_MOL_ALG_SOURCES
filter_clashes.cc
construct_cbeta.cc
trajectory_analysis.cc
+ structure_analysis.cc
)
set(MOL_ALG_DEPS ost_mol ost_seq)
diff --git a/modules/mol/alg/src/structure_analysis.cc b/modules/mol/alg/src/structure_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..73bfca844b831e7afac45ef09ec73d0f4345233b
--- /dev/null
+++ b/modules/mol/alg/src/structure_analysis.cc
@@ -0,0 +1,83 @@
+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+
+/*
+ * Author Niklaus Johner
+ */
+
+#include <ost/base.hh>
+#include <ost/mol/mol.hh>
+#include "structure_analysis.hh"
+
+namespace ost { namespace mol { namespace alg {
+
+geom::Vec3List GetPosListFromView(const EntityView& view){
+ CheckHandleValidity(view);
+ geom::Vec3List vl;
+ AtomViewList atoms=view.GetAtomList();
+ vl.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ vl.push_back(i->GetPos());
+ }
+ return vl;
+}
+
+std::vector<Real> CalculateAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map){
+ CheckHandleValidity(density_map);
+ std::vector<Real> v;
+ v.reserve(vl.size());
+ for (geom::Vec3List::const_iterator v1=vl.begin(),e=vl.end(); v1!=e; ++v1) {
+ img::Point p(density_map.CoordToIndex(*v1));
+ v.push_back(density_map.GetReal(p));
+ }
+ return v;
+}
+
+Real CalculateAverageAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map){
+ CheckHandleValidity(density_map);
+ std::vector<Real> v=CalculateAgreementWithDensityMap(vl, density_map);
+ Real sum=0.0;
+ for (std::vector<Real>::const_iterator i=v.begin(),e=v.end(); i!=e; ++i) {
+ sum=sum+*i;
+ }
+ return sum/float(vl.size());
+}
+
+void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom::Vec3 cell_center, const geom::Vec3 basis_vec){
+ mol::XCSEditor edi=eh.EditXCS(mol::BUFFERED_EDIT);
+ geom::Vec3 cm,wrapped_cm,shift;
+ edi=eh.EditXCS();
+ ResidueHandleList residues=eh.GetResidueList();
+ unsigned int n_residues=eh.GetResidueCount();
+ for (unsigned int i=0; i<n_residues; ++i) {
+ ResidueHandle r=residues[i];
+ cm=r.GetCenterOfMass();
+ wrapped_cm=geom::WrapVec3(cm,cell_center,basis_vec);
+ if (wrapped_cm==cm) continue;
+ AtomHandleList atoms=r.GetAtomList();
+ unsigned int n_atoms=r.GetAtomCount();
+ shift=wrapped_cm-cm;
+ for (unsigned int j=0; j<n_atoms; ++j) {
+ edi.SetAtomPos(atoms[j],atoms[j].GetPos()+shift);
+ }
+ }
+}
+
+}}} //ns
diff --git a/modules/mol/alg/src/structure_analysis.hh b/modules/mol/alg/src/structure_analysis.hh
new file mode 100644
index 0000000000000000000000000000000000000000..798dc09739f93f8c9cf5921119be590a103c4458
--- /dev/null
+++ b/modules/mol/alg/src/structure_analysis.hh
@@ -0,0 +1,39 @@
+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+
+/*
+ * Niklaus Johner
+ */
+#ifndef OST_STRUCTURE_ANALYSIS_HH
+#define OST_STRUCTURE_ANALYSIS_HH
+
+#include <ost/mol/alg/module_config.hh>
+
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/entity_handle.hh>
+#include <ost/img/map.hh>
+
+
+namespace ost { namespace mol { namespace alg {
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG GetPosListFromView(const EntityView& view);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG CalculateAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map);
+ Real DLLEXPORT_OST_MOL_ALG CalculateAverageAgreementWithDensityMap(const geom::Vec3List& vl, img::MapHandle& density_map);
+ void DLLEXPORT_OST_MOL_ALG WrapEntityInPeriodicCell(EntityHandle eh, const geom::Vec3 cell_center, const geom::Vec3 basis_vec);
+}}}//ns
+#endif
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index 6636421ad702ab8dc41f88f865ae8de3fd2ab451..3f760c65c78491c75001a6348d7c246c2a4c733d 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/modules/mol/base/pymod/export_coord_frame.cc
@@ -26,6 +26,8 @@ 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;
@@ -35,9 +37,12 @@ Real (CoordFrame::*get_dist_cm)(const mol::EntityView& sele1, const mol::EntityV
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;
+
void export_CoordFrame()
{
class_<CoordFrame>("CoordFrame",no_init)
+ .def("GetCellSize",GetCellSize)
+ .def("GetCellAngles",GetCellAngles)
.def("GetAtomPos", get_atom_pos)
.def("GetDistanceBetwAtoms", get_dist_atom)
.def("GetAngle", get_angle)
diff --git a/modules/mol/base/pymod/export_coord_group.cc b/modules/mol/base/pymod/export_coord_group.cc
index fa3a06132292cc2e65d727197e2888832c54da0b..6bdeb1d6958f797898aea3a2915a3063daf10dea 100644
--- a/modules/mol/base/pymod/export_coord_group.cc
+++ b/modules/mol/base/pymod/export_coord_group.cc
@@ -38,6 +38,8 @@ namespace {
void (CoordGroupHandle::*capture1)()=&CoordGroupHandle::Capture;
void (CoordGroupHandle::*capture2)(uint)=&CoordGroupHandle::Capture;
+ void (CoordGroupHandle::*add_frame1)(const geom::Vec3List&)=&CoordGroupHandle::AddFrame;
+ void (CoordGroupHandle::*add_frame2)(const geom::Vec3List&,const geom::Vec3&,const geom::Vec3&)=&CoordGroupHandle::AddFrame;
}
void export_CoordGroup()
@@ -46,8 +48,10 @@ void export_CoordGroup()
.def("IsValid",&CoordGroupHandle::IsValid)
.def("GetEntity",&CoordGroupHandle::GetEntity)
.def("GetAtomCount",&CoordGroupHandle::GetAtomCount)
+ .def("GetFrame",&CoordGroupHandle::GetFrame2)
.def("AddFrames", &CoordGroupHandle::AddFrames)
- .def("AddFrame", &CoordGroupHandle::AddFrame)
+ .def("AddFrame", add_frame1)
+ .def("AddFrame", add_frame2)
.def("GetFrameCount",&CoordGroupHandle::GetFrameCount)
.def("GetFramePositions",&CoordGroupHandle::GetFramePositions)
.def("SetFramePositions",&CoordGroupHandle::SetFramePositions)
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index dbd51c29b6030dfdf687b590824e48678a628fc0..3c5b2b28040628c7cf0c83fc0e25da4d3d6de2c7 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -31,15 +31,28 @@
namespace ost { namespace mol {
class DLLEXPORT_OST_MOL CoordFrame : public geom::Vec3List{
+private:
+ geom::Vec3 periodic_cell_sizes_;
+ geom::Vec3 periodic_cell_angles_;
public:
typedef geom::Vec3List base_type;
- CoordFrame() : base_type() {}
+ CoordFrame() : base_type() {}
CoordFrame(size_t size, const geom::Vec3& value=geom::Vec3()) : base_type(size, value) {}
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;
+ }
+ geom::Vec3 GetCellSize() {
+ return this->periodic_cell_sizes_;
+ }
+ geom::Vec3 GetCellAngles() {
+ return this->periodic_cell_angles_;
+ }
geom::Vec3 GetAtomPos(const AtomHandle& atom);
geom::Vec3 GetAtomPos(int atom_index);
Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2);
diff --git a/modules/mol/base/src/coord_group.cc b/modules/mol/base/src/coord_group.cc
index 127c82463bf20c1df65f2f0234d3474bdc3c78ed..db3fd341ce04a76c6e928ca1e0fd93e869b583c7 100644
--- a/modules/mol/base/src/coord_group.cc
+++ b/modules/mol/base/src/coord_group.cc
@@ -103,6 +103,16 @@ CoordGroupHandle::operator bool() const
}
}
+ void CoordGroupHandle::AddFrame(const geom::Vec3List& clist, const geom::Vec3& cell_size, const geom::Vec3& cell_angles)
+{
+ this->CheckValidity();
+ if (source_->IsMutable()) {
+ source_->AddFrame(clist,cell_size,cell_angles);
+ } else {
+ throw IntegrityError("Can't add frame to immutable CoordGroup");
+ }
+}
+
void CoordGroupHandle::AddFrames(const CoordGroupHandle& cg)
{
this->CheckValidity();
@@ -143,6 +153,12 @@ CoordFramePtr CoordGroupHandle::GetFrame(uint frame) const
return source_->GetFrame(frame);
}
+CoordFrame CoordGroupHandle::GetFrame2(uint frame)
+{
+ this->CheckValidity();
+ return *(source_->GetFrame(frame));
+}
+
AtomHandleList CoordGroupHandle::GetAtomList() const
{
this->CheckValidity();
diff --git a/modules/mol/base/src/coord_group.hh b/modules/mol/base/src/coord_group.hh
index 22a829dd662f50d8bee3cb6490806ca120039f01..652fe34636d3dc11c7fa5168c83c05d906b6524a 100644
--- a/modules/mol/base/src/coord_group.hh
+++ b/modules/mol/base/src/coord_group.hh
@@ -69,7 +69,8 @@ public:
/// \brief add frame
//void AddFrame(const std::vector<geom::Vec3>& clist);
void AddFrame(const geom::Vec3List& clist);
-
+ void AddFrame(const geom::Vec3List& clist,const geom::Vec3& cell_size,const geom::Vec3& cell_angles);
+
void AddFrames(const CoordGroupHandle& cg);
/// \brief set an indidivial atom position in the given frame
void SetAtomPos(uint frame, AtomHandle atom, const geom::Vec3& pos);
@@ -85,6 +86,7 @@ public:
AtomHandleList GetAtomList() const;
CoordFramePtr GetFrame(uint frame) const;
+ CoordFrame GetFrame2(uint frame);
/// \brief return a filtered coord group, containing only the atoms in the
/// view
diff --git a/modules/mol/base/src/coord_source.hh b/modules/mol/base/src/coord_source.hh
index 604817535d9bc2fc0bd12c20739135e287642862..fec0d9946cb2ecd275abb23d5ebb58b080422734 100644
--- a/modules/mol/base/src/coord_source.hh
+++ b/modules/mol/base/src/coord_source.hh
@@ -70,6 +70,7 @@ public:
void Capture(uint f);
virtual void AddFrame(const std::vector<geom::Vec3>& coords) = 0;
+ virtual void AddFrame(const std::vector<geom::Vec3>& coords,const geom::Vec3& cell_size,const geom::Vec3& cell_angles) = 0;
virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords) = 0;
protected:
void SetMutable(bool flag);
diff --git a/modules/mol/base/src/in_mem_coord_source.cc b/modules/mol/base/src/in_mem_coord_source.cc
index 1690c11048713795f5bf3d86bbef29b9fb82702f..ea01db894a69efa04189e126b3f2d6429df46f12 100644
--- a/modules/mol/base/src/in_mem_coord_source.cc
+++ b/modules/mol/base/src/in_mem_coord_source.cc
@@ -30,6 +30,12 @@ void InMemCoordSource::AddFrame(const std::vector<geom::Vec3>& coords)
frames_.push_back(fp);
}
+void InMemCoordSource::AddFrame(const std::vector<geom::Vec3>& coords,const geom::Vec3& cell_size,const geom::Vec3& cell_angles)
+{
+ CoordFramePtr fp(new CoordFrame(coords,cell_size,cell_angles));
+ frames_.push_back(fp);
+}
+
void InMemCoordSource::InsertFrame(int pos, const std::vector<geom::Vec3>& coords)
{
CoordFrameList::iterator it = frames_.begin();
diff --git a/modules/mol/base/src/in_mem_coord_source.hh b/modules/mol/base/src/in_mem_coord_source.hh
index dc032ebbea294fdc6b4bb7ba86d79847b5e818ca..0984791aae634bcf315efe08f55f10850ef77c5f 100644
--- a/modules/mol/base/src/in_mem_coord_source.hh
+++ b/modules/mol/base/src/in_mem_coord_source.hh
@@ -44,6 +44,7 @@ public:
void AddFrame(const CoordFramePtr& frame);
virtual void AddFrame(const std::vector<geom::Vec3>& coords);
+ virtual void AddFrame(const std::vector<geom::Vec3>& coords,const geom::Vec3& cell_size,const geom::Vec3& cell_angles);
virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords);
private: