diff --git a/.gitignore b/.gitignore
index 34ea03ce96564582673101299044af3e7792dac1..53babcd731562cfbea3ea0cfbfd8222dfabc2a89 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,4 @@ CPackConfig.cmake
CPackSourceConfig.cmake
deployment/macos/DNG.app
deployment/macos/standalone
+ChemLib
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5ed6abafccc0d47d949b11811f9a45fdebdddf84..08575372b9f41c3d2b6091bad61327232eb529c2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,9 +8,6 @@ include(OST)
project(OpenStructure CXX C)
-INCLUDE(CTest)
-ENABLE_TESTING()
-
option(USE_SHADER "whether to compile with shader support"
OFF)
option(SET_RPATH "embed rpath upon make install"
@@ -112,7 +109,6 @@ if (NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE Debug)
endif()
-enable_testing()
set(STAGE_DIR "${CMAKE_BINARY_DIR}/stage")
set(EXECUTABLE_OUTPUT_PATH ${STAGE_DIR}/bin)
set(HEADER_STAGE_PATH ${STAGE_DIR}/include)
diff --git a/modules/base/pymod/export_logger.cc b/modules/base/pymod/export_logger.cc
index 4c2b79b20f3b95dcdf8f0576a3bd166c49916a37..d0b0800224b65da68a2082729e181e271d891389 100644
--- a/modules/base/pymod/export_logger.cc
+++ b/modules/base/pymod/export_logger.cc
@@ -60,6 +60,12 @@ void pop_verb()
Logger::Instance().PopVerbosityLevel();
}
+int get_verb()
+{
+ return Logger::Instance().GetVerbosityLevel();
+}
+
+
void push_log_sink(LogSinkPtr sink)
{
Logger::Instance().PushSink(sink);
@@ -105,6 +111,7 @@ void export_Logger()
def("PushVerbosityLevel",push_verb);
def("PopVerbosityLevel",pop_verb);
+ def("GetVerbosityLevel",get_verb);
def("PushLogSink",push_log_sink);
def("PopLogSink",pop_log_sink);
def("LogError",log_error);
diff --git a/modules/base/src/log.hh b/modules/base/src/log.hh
index b535346f3caeae6f56c2809d5453a0ca78ef3ff6..4b129767d7462b6c9e0ba5f8b3d4b5ca9e620a1b 100644
--- a/modules/base/src/log.hh
+++ b/modules/base/src/log.hh
@@ -52,7 +52,7 @@ public:
static Logger& Instance();
LogSinkPtr GetCurrentSink() { return sink_stack_.top(); }
- int GetLogLevel() const {return level_;}
+ int GetVerbosityLevel() const {return level_;}
void ResetSinks() {
while (sink_stack_.size()>1) {
@@ -71,7 +71,7 @@ private:
};
-#define OST_DO_LOGGING_(m, l) if (::ost::Logger::Instance().GetLogLevel()>=l) {\
+#define OST_DO_LOGGING_(m, l) if (::ost::Logger::Instance().GetVerbosityLevel()>=l) {\
std::stringstream tmp_s__; \
tmp_s__ << m << std::endl; \
::ost::Logger::Instance().GetCurrentSink()->LogMessage(tmp_s__.str(), l); \
diff --git a/modules/conop/data/charmm.cif b/modules/conop/data/charmm.cif
index ced6449d30ad68b02a4026cf8e752a479637d445..5368434308251855fa93842ebf3db635bceb3867 100644
--- a/modules/conop/data/charmm.cif
+++ b/modules/conop/data/charmm.cif
@@ -1040,7 +1040,7 @@ ILE O O O 0 N Y N O ILE 4
ILE CB CB C 0 N N N CB ILE 5
ILE CG2 CG2 C 0 N N N CG2 ILE 6
ILE CG1 CG1 C 0 N N N CG1 ILE 7
-ILE CD CD C 0 N N N CD ILE 8
+ILE CD CD1 C 0 N N N CD ILE 8
ILE HA HA H 0 N N N HA ILE 9
ILE HB HB H 0 N N N HB ILE 10
ILE HG21 HG21 H 0 N N N HG21 ILE 11
diff --git a/modules/conop/pymod/export_compound.cc b/modules/conop/pymod/export_compound.cc
index b0c6e11d59655283b781ce5d34ff6d2cda3daba1..379f2059832561f9aef461589378e0571950b193 100644
--- a/modules/conop/pymod/export_compound.cc
+++ b/modules/conop/pymod/export_compound.cc
@@ -51,7 +51,7 @@ void export_Compound() {
register_ptr_to_python<CompoundPtr>();
class_<CompoundLib>("CompoundLib", no_init)
- .def("Load", &CompoundLib::Load).staticmethod("Load")
+ .def("Load", &CompoundLib::Load, arg("readonly")=true).staticmethod("Load")
.def("FindCompound", &CompoundLib::FindCompound)
.def("ClearCache", &CompoundLib::ClearCache)
;
diff --git a/modules/conop/src/compound_lib.cc b/modules/conop/src/compound_lib.cc
index 10e46b16ba0dabe90274699accd01d8c832c548b..b1878a973fd1c75913dbae05d6acca2cac9162d2 100644
--- a/modules/conop/src/compound_lib.cc
+++ b/modules/conop/src/compound_lib.cc
@@ -237,10 +237,11 @@ CompoundLibPtr CompoundLib::Create(const String& database)
}
-CompoundLibPtr CompoundLib::Load(const String& database)
+CompoundLibPtr CompoundLib::Load(const String& database, bool readonly)
{
+ int flags=readonly ? SQLITE_OPEN_READONLY : SQLITE_OPEN_READWRITE;
CompoundLibPtr lib(new CompoundLib);
- int retval=sqlite3_open(database.c_str(), &lib->conn_);
+ int retval=sqlite3_open_v2(database.c_str(), &lib->conn_, flags, NULL);
if (SQLITE_OK==retval) {
return lib;
}
diff --git a/modules/conop/src/compound_lib.hh b/modules/conop/src/compound_lib.hh
index 8948e3479585f61e0fd9fa3d8b4b5c747fa97ccf..3d642d512c6d50ae4c0c92ae6ee9edbbae7bdedc 100644
--- a/modules/conop/src/compound_lib.hh
+++ b/modules/conop/src/compound_lib.hh
@@ -37,7 +37,7 @@ typedef std::map<String, CompoundPtr> CompoundMap;
class DLLEXPORT_OST_CONOP CompoundLib {
public:
- static CompoundLibPtr Load(const String& database);
+ static CompoundLibPtr Load(const String& database, bool readonly=true);
static CompoundLibPtr Create(const String& database);
~CompoundLib();
diff --git a/modules/io/pymod/__init__.py b/modules/io/pymod/__init__.py
index 2cf6900a10edd878f78a2b748830d33f2f102afe..5534c540dd5f1e025ff22f791577e5152f81b021 100644
--- a/modules/io/pymod/__init__.py
+++ b/modules/io/pymod/__init__.py
@@ -187,11 +187,38 @@ def LoadImageList (files):
LoadMapList=LoadImageList
+def LoadCHARMMTraj(crd, dcd_file=None, lazy_load=False, stride=1,
+ dialect='CHARMM'):
+ """
+ Load CHARMM trajectory file.
+
+ :param crd: EntityHandle or filename of the CRD (PDB) file containing the
+ structure. The structure must have the same number of atoms as the
+ trajectory
+ :param dcd_file: The filename of the DCD file. If not set, and crd is a string,
+ the filename is set to the <crd>.dcd
+ :param layz_load: Whether the trajectory should be loaded on demand. Instead of
+ loading the complete trajectory into memory, the trajectory frames are
+ loaded from disk when requested.
+ :param stride: The spacing of the frames to load. When set to 2, for example,
+ every second frame is loaded from the trajectory
+ """
+ if not isinstance(crd, mol.EntityHandle):
+ if dcd_file==None:
+ dcd_file='%s.dcd' % os.path.splitext(crd)[0]
+ crd=LoadPDB(crd, dialect=dialect)
+
+ else:
+ if not dcd_file:
+ raise ValueError("No DCD filename given")
+ return LoadCHARMMTraj_(crd, dcd_file, stride, lazy_load)
+
## \example fft_li.py
#
-# This scripts loads one or more images and shows their Fourier Transforms on the screen. A viewer
-# is opened for each loaded image. The Fourier Transform honors the origin of the reference system,
-# which is assumed to be at the center of the image.
+# This scripts loads one or more images and shows their Fourier Transforms on
+# the screen. A viewer is opened for each loaded image. The Fourier Transform
+# honors the origin of the reference system, which is assumed to be at the
+# center of the image.
#
# Usage:
#
diff --git a/modules/io/pymod/wrap_io.cc b/modules/io/pymod/wrap_io.cc
index dc5ce987c38cb9dd4f2cd692cc306b8b2dd567a4..6eefe172870f52e934a27d6e0f204b6d868ec61b 100644
--- a/modules/io/pymod/wrap_io.cc
+++ b/modules/io/pymod/wrap_io.cc
@@ -69,11 +69,6 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_view_ov,
BOOST_PYTHON_FUNCTION_OVERLOADS(save_charmm_trj_ov,
SaveCHARMMTraj, 3, 4)
-mol::CoordGroupHandle load_dcd1(const String& c, const String& t) {return LoadCHARMMTraj(c,t);}
-mol::CoordGroupHandle load_dcd2(const String& c, const String& t, unsigned int s) {return LoadCHARMMTraj(c,t,s);}
-mol::CoordGroupHandle load_dcd3(const mol::EntityHandle& e, const String& t) {return LoadCHARMMTraj(e,t);}
-mol::CoordGroupHandle load_dcd4(const mol::EntityHandle& e, const String& t, unsigned int s) {return LoadCHARMMTraj(e,t,s);}
-
}
void export_pdb_io();
@@ -115,10 +110,8 @@ BOOST_PYTHON_MODULE(_io)
def("LoadSDF", &LoadSDF);
def("LoadCRD", &LoadCRD);
- def("LoadCHARMMTraj",load_dcd1);
- def("LoadCHARMMTraj",load_dcd2);
- def("LoadCHARMMTraj",load_dcd3);
- def("LoadCHARMMTraj",load_dcd4);
+ def("LoadCHARMMTraj_", &LoadCHARMMTraj, (arg("ent"), arg("trj_filename"),
+ arg("stride")=1, arg("lazy_load")=false));
def("SaveCHARMMTraj",SaveCHARMMTraj,save_charmm_trj_ov());
def("LoadMAE", &LoadMAE);
diff --git a/modules/io/src/img/map_io_mrc_handler.cc b/modules/io/src/img/map_io_mrc_handler.cc
index aa4a1cf40be7817932edd7662357bb52722f6e73..f5e5d1dc7d7333db5acdf60678678bae507ac400 100644
--- a/modules/io/src/img/map_io_mrc_handler.cc
+++ b/modules/io/src/img/map_io_mrc_handler.cc
@@ -775,7 +775,7 @@ void import_helper(img::MapHandle& image, std::istream& in,const MRC& formatmrc)
f >> header;
if(header.mode==5) header.mode=0;
- if(Logger::Instance().GetLogLevel()>=4) {
+ if(Logger::Instance().GetVerbosityLevel()>=4) {
header.Print();
}
if(header.mode==3 || header.mode==4) {
@@ -864,7 +864,7 @@ void export_helper(const img::MapHandle& image,
throw(IOException("MRC/CCP4 export: full complex export not supported."));
}
}
- if(Logger::Instance().GetLogLevel()>=4) {
+ if(Logger::Instance().GetVerbosityLevel()>=4) {
header.Print();
}
}
diff --git a/modules/io/src/img/map_io_nanoscope_handler.cc b/modules/io/src/img/map_io_nanoscope_handler.cc
index a18a25c56ee5d52050d3ae5637d0b978056836c6..1ade5a6f0e73f5b5103224d2ca63450bcab3797e 100644
--- a/modules/io/src/img/map_io_nanoscope_handler.cc
+++ b/modules/io/src/img/map_io_nanoscope_handler.cc
@@ -178,7 +178,7 @@ void header_filler(NHeader& h, std::istream& in, int inum)
void print_header (const NHeader& header, int inum)
{
- if(Logger::Instance().GetLogLevel()>2) {
+ if(Logger::Instance().GetVerbosityLevel()>2) {
LOG_INFO("io_nanoscope: header dump for image " << inum);
LOG_INFO(" px : " << header.px);
LOG_INFO(" py : " << header.py);
diff --git a/modules/io/src/img/map_io_tiff_handler.cc b/modules/io/src/img/map_io_tiff_handler.cc
index a4c3fd191312e5ec94caa4eb11a4dda00a97e2f4..b80256d0f7908ab9af8c23a04c5ce2253f222d91 100644
--- a/modules/io/src/img/map_io_tiff_handler.cc
+++ b/modules/io/src/img/map_io_tiff_handler.cc
@@ -450,7 +450,7 @@ TIFF* MapIOTiffHandler::open_subimage_file(const boost::filesystem::path& locati
if(directory_number>0) {
TIFFSetDirectory(tfile,directory_number);
}
- if(Logger::Instance().GetLogLevel()>=5) {
+ if(Logger::Instance().GetVerbosityLevel()>=5) {
TIFFPrintDirectory(tfile,stderr,directory_number);
}
return tfile;
@@ -481,7 +481,7 @@ TIFF* MapIOTiffHandler::open_subimage_stream(std::istream& location,const TIF& f
if(directory_number>0) {
TIFFSetDirectory(tfile,directory_number);
}
- if(Logger::Instance().GetLogLevel()>=5) {
+ if(Logger::Instance().GetVerbosityLevel()>=5) {
TIFFPrintDirectory(tfile,stderr,directory_number);
}
return tfile;
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index 177cd57af0e7a7144d40123c17c091d675f77c7f..86a59ed42d845a1a215892694bd3ba088d2cb8d9 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/modules/io/src/mol/dcd_io.cc
@@ -66,42 +66,30 @@ bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
return a1.GetIndex()<a2.GetIndex();
}
-mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
- const String& trj_fn,
- unsigned int stride)
+bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
+ bool& skip_flag, bool& gap_flag)
{
- Profile profile_load("LoadCHARMMTraj");
-
- mol::AtomHandleList alist(alist2);
- std::sort(alist.begin(),alist.end(),less_index);
-
- bool gap_flag = true;
-
- boost::filesystem::path trj_f(trj_fn);
- boost::filesystem::ifstream ff(trj_f, std::ios::binary);
-
- DCDHeader header;
- char dummy[4];
- bool swap_flag=false;
-
- LOG_INFO("importing trajectory data");
-
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(header.hdrr,sizeof(char)*4);
- ff.read(reinterpret_cast<char*>(header.icntrl),sizeof(int)*20);
+ if (!istream) {
+ return false;
+ }
+ char dummy[4];
+ gap_flag=true;
+ swap_flag=false;
+ skip_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);
if(header.icntrl[1]<0 || header.icntrl[1]>1e8) {
// nonsense atom count, try swapping
swap_int(header.icntrl,20);
if(header.icntrl[1]<0 || header.icntrl[1]>1e8) {
throw(IOException("LoadCHARMMTraj: nonsense atom count in header"));
} else {
- LOG_INFO("LoadCHARMMTraj: byte-swapping");
+ LOG_VERBOSE("LoadCHARMMTraj: byte-swapping");
swap_flag=true;
}
}
- bool skip_flag=false;
-
if(header.icntrl[19]!=0) { // CHARMM format
skip_flag=(header.icntrl[10]!=0);
if(skip_flag) {
@@ -113,24 +101,102 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
// XPLOR format
LOG_VERBOSE("LoadCHARMMTraj: using XPLOR format");
}
-
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(reinterpret_cast<char*>(&header.ntitle),sizeof(int));
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ if (!istream) {
+ return false;
+ }
+ istream.read(reinterpret_cast<char*>(&header.ntitle),sizeof(int));
+ if (!istream) {
+ return false;
+ }
if(swap_flag) swap_int(&header.ntitle,1);
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(header.title,sizeof(char)*header.ntitle);
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+
+ istream.read(header.title,sizeof(char)*header.ntitle);
header.title[header.ntitle]='\0';
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(reinterpret_cast<char*>(&header.t_atom_count),sizeof(int));
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ istream.read(reinterpret_cast<char*>(&header.t_atom_count),sizeof(int));
if(swap_flag) swap_int(&header.t_atom_count,1);
- if(gap_flag) ff.read(dummy,sizeof(dummy));
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
header.num=header.icntrl[0];
header.istep=header.icntrl[1];
header.freq=header.icntrl[2];
header.nstep=header.icntrl[3];
header.f_atom_count=header.icntrl[8];
header.atom_count=header.t_atom_count-header.f_atom_count;
+ return true;
+}
+
+
+size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
+{
+ size_t frame_size=0;
+ if (skip_flag) {
+ frame_size+=14*sizeof(int);
+ }
+ if (gap_flag) {
+ frame_size+=6*sizeof(int);
+ }
+ frame_size+=3*sizeof(float)*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)
+{
+ char dummy[4];
+ if(skip_flag) istream.seekg(14*4,std::ios_base::cur);
+ // read each frame
+ if(!istream) {
+ /* premature EOF */
+ LOG_ERROR("LoadCHARMMTraj: premature end of file, frames read");
+ return false;
+ }
+ // x coord
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ istream.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+ if(swap_flag) swap_float(&xlist[0],xlist.size());
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ for(uint j=0;j<frame.size();++j) {
+ frame[j].x=xlist[j];
+ }
+
+ // y coord
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ istream.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+ if(swap_flag) swap_float(&xlist[0],xlist.size());
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ for(uint j=0;j<frame.size();++j) {
+ frame[j].y=xlist[j];
+ }
+
+ // z coord
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ istream.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+ if(swap_flag) swap_float(&xlist[0],xlist.size());
+ if(gap_flag) istream.read(dummy,sizeof(dummy));
+ for(uint j=0;j<frame.size();++j) {
+ frame[j].z=xlist[j];
+ }
+ return true;
+}
+
+
+mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
+ const String& trj_fn,
+ unsigned int stride)
+{
+ Profile profile_load("LoadCHARMMTraj");
+ mol::AtomHandleList alist(alist2);
+ std::sort(alist.begin(),alist.end(),less_index);
+
+ std::ifstream istream(trj_fn.c_str(), std::ios::binary);
+ DCDHeader header;
+ bool swap_flag=false, skip_flag=false, gap_flag=false;
+ read_dcd_header(istream, header, swap_flag, skip_flag, gap_flag);
LOG_DEBUG("LoadCHARMMTraj: " << header.num << " trajectories with "
<< header.atom_count << " atoms (" << header.f_atom_count
<< " fixed) each");
@@ -145,63 +211,23 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
mol::CoordGroupHandle cg=CreateCoordGroup(alist);
std::vector<geom::Vec3> clist(header.t_atom_count);
std::vector<float> xlist(header.t_atom_count);
-
- size_t frame_size=0;
- if (skip_flag) {
- frame_size+=14*4;
- }
- if (gap_flag) {
- frame_size+=6*sizeof(dummy);
- }
- frame_size+=3*sizeof(float)*xlist.size();
-
+ size_t frame_size=calc_frame_size(skip_flag, gap_flag, xlist.size());
int i=0;
for(;i<header.num;i+=stride) {
- if(skip_flag) ff.seekg(14*4,std::ios_base::cur);
- // read each frame
- if(!ff) {
- /* premature EOF */
- LOG_ERROR("LoadCHARMMTraj: premature end of file, " << i
- << " frames read");
+ if (!read_frame(istream, header, frame_size, skip_flag, gap_flag,
+ swap_flag, xlist, clist)) {
break;
}
- // x coord
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
- if(swap_flag) swap_float(&xlist[0],xlist.size());
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- for(uint j=0;j<clist.size();++j) {
- clist[j].x=xlist[j];
- }
-
- // y coord
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
- if(swap_flag) swap_float(&xlist[0],xlist.size());
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- for(uint j=0;j<clist.size();++j) {
- clist[j].y=xlist[j];
- }
-
- // z coord
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
- if(swap_flag) swap_float(&xlist[0],xlist.size());
- if(gap_flag) ff.read(dummy,sizeof(dummy));
- for(uint j=0;j<clist.size();++j) {
- clist[j].z=xlist[j];
- }
-
cg.AddFrame(clist);
// skip frames (defined by stride)
- if(stride>1) ff.seekg(frame_size*(stride-1),std::ios_base::cur);
+ if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur);
}
- ff.get();
- if(!ff.eof()) {
+ istream.get();
+ if(!istream.eof()) {
LOG_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: "
- << ff.tellg());
+ << istream.tellg());
}
LOG_VERBOSE("Loaded " << cg.GetFrameCount() << " frames with " << cg.GetAtomCount() << " atoms each");
@@ -209,25 +235,89 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
return cg;
}
-} // anon ns
+class DCDCoordSource : public mol::CoordSource {
+public:
+ DCDCoordSource(const mol::AtomHandleList& atoms, const String& filename,
+ uint stride):
+ mol::CoordSource(atoms), filename_(filename),
+ stream_(filename.c_str(), std::ios::binary), loaded_(false), stride_(stride)
+ {
+ this->SetMutable(false);
+ frame_=mol::CoordFramePtr(new mol::CoordFrame(atoms.size()));
+ }
+
+
+ virtual uint GetFrameCount()
+ {
+ if (!frame_count_)
+ const_cast<DCDCoordSource*>(this)->FetchFrame(0);
+ return frame_count_;
+ }
+
+ virtual mol::CoordFramePtr GetFrame(uint frame_id) const {
+ const_cast<DCDCoordSource*>(this)->FetchFrame(frame_id);
+ return frame_;
+ }
-mol::CoordGroupHandle LoadCHARMMTraj(const String& crd_fn,
- const String& trj_fn,
- unsigned int stride)
+ virtual void AddFrame(const std::vector<geom::Vec3>& coords) {}
+ virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords) {}
+private:
+
+ void FetchFrame(uint frame);
+ String filename_;
+ DCDHeader header_;
+ bool skip_flag_;
+ bool swap_flag_;
+ bool gap_flag_;
+ std::ifstream stream_;
+ bool loaded_;
+ uint frame_count_;
+ uint curr_frame_;
+ uint stride_;
+ size_t frame_start_;
+ mol::CoordFramePtr frame_;
+};
+
+
+void DCDCoordSource::FetchFrame(uint frame)
{
- mol::EntityHandle ent=LoadEntity(crd_fn);
- return load_dcd(ent.GetAtomList(),trj_fn,stride);
+ if (!loaded_) {
+ read_dcd_header(stream_, header_, swap_flag_, skip_flag_, gap_flag_);
+ frame_start_=stream_.tellg();
+ loaded_=true;
+ frame_count_=header_.num;
+ }
+ size_t frame_size=calc_frame_size(skip_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_,
+ swap_flag_, xlist, *frame_.get())) {
+ }
}
+typedef boost::shared_ptr<DCDCoordSource> DCDCoordSourcePtr;
+
+
+} // anon ns
+
-mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& e,
- const String& trj_fn,
- unsigned int stride)
+mol::CoordGroupHandle LoadCHARMMTraj(const mol::EntityHandle& ent,
+ const String& trj_fn,
+ unsigned int stride, bool lazy_load)
{
- return load_dcd(e.GetAtomList(),trj_fn,stride);
+ mol::AtomHandleList alist(ent.GetAtomList());
+ std::sort(alist.begin(),alist.end(),less_index);
+ if (lazy_load) {
+ LOG_INFO("Importing CHARMM trajectory with lazy_load=true");
+ DCDCoordSource* source=new DCDCoordSource(alist, trj_fn, stride);
+ return mol::CoordGroupHandle(DCDCoordSourcePtr(source));
+ }
+ LOG_INFO("Importing CHARMM trajectory with lazy_load=false");
+ return load_dcd(alist, trj_fn, stride);
}
-
namespace {
void write_dcd_hdr(std::ofstream& out,
diff --git a/modules/io/src/mol/dcd_io.hh b/modules/io/src/mol/dcd_io.hh
index b20a3fb0f4b7bcbb2ccf6b47b1b206fdeb3b3e57..24d189e634e5f7ec0289eb09680f6debf039a45a 100644
--- a/modules/io/src/mol/dcd_io.hh
+++ b/modules/io/src/mol/dcd_io.hh
@@ -29,22 +29,14 @@
namespace ost { namespace io {
-/*! \brief import a CHARMM trajectory in dcd format
- requires the coordinate and the trajectory file; the format
- of the coordinate file will be automatically deduced from the extension
- the optional stride parameter will cause only every nth frame to be loaded
-*/
-mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const String& coord,
- const String& trj,
- unsigned int stride=1);
-
/*! \brief import a CHARMM trajectory in dcd format with an existing entity
requires the existing entity and the trajectory file - obviously the
atom layout of the entity must match the trajectory file
*/
-mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& e,
- const String& trj,
- unsigned int stride=1);
+mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& ent,
+ const String& trj_filename,
+ unsigned int stride=1,
+ bool lazy_load=false);
/*! \brief export coord group as PDB file and DCD trajectory
diff --git a/modules/mol/alg/molalg.rst b/modules/mol/alg/molalg.rst
index 00735f6def328ed899de4005fb2f803cd23f0206..2daf3ca7597bca27d86b9d86d53e7c67a33b5bb6 100644
--- a/modules/mol/alg/molalg.rst
+++ b/modules/mol/alg/molalg.rst
@@ -22,4 +22,22 @@
naming of the atoms is ambigous. For these residues, the overlap of both
possible solutions to the fixed atoms, that is, everything that is not
ambigous is calculated. The solution that gives higher overlap is then used to
- calculate the actual overlap score.
\ No newline at end of file
+ calculate the actual overlap score.
+
+
+.. function:: SuperposeFrames(frames, sel, from=0, to=-1, ref=-1)
+
+ This function superposes the frames of the given coord group and returns them
+ as a new coord group.
+
+ :param frames: The source coord group.
+ :param from: index of the first frame
+ :param to: index of the last frame plus one. If set to -1, the value is set to
+ the number of frames in the coord group
+ :param sel: An entity view containing the selection of atoms to be used for
+ superposition. If set to an invalid view, all atoms in the coord group are
+ used.
+ :param ref: The index of the reference frame to use for superposition. If set
+ to -1, the each frame is superposed to the previous frame.
+
+ :returns: A newly created coord group containing the superposed frames.
\ 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 240a029a42dd523e8a4dfab72be40fe77bd74fa7..67f4fbf9e13454d105b0f69844297b06a1508309 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -23,7 +23,7 @@
#include <boost/python.hpp>
#include <ost/config.hh>
#include <ost/mol/alg/local_dist_test.hh>
-
+#include <ost/mol/alg/superpose_frames.hh>
using namespace boost::python;
void export_svdSuperPose();
@@ -40,4 +40,7 @@ BOOST_PYTHON_MODULE(_mol_alg)
#endif
def("LocalDistTest", &ost::mol::alg::LocalDistTest);
+ def("SuperposeFrames", &ost::mol::alg::SuperposeFrames,
+ (arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
+ arg("end")=-1, arg("ref")=-1));
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 17dc9acbf2a5d32ac17e62b0e50af18f35122b2f..e613a9c735cff3588598ae612eec5ff0ce805202 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -3,12 +3,14 @@ set(OST_MOL_ALG_HEADERS
module_config.hh
sec_structure_segments.hh
local_dist_test.hh
+ superpose_frames.hh
)
set(OST_MOL_ALG_SOURCES
svd_superpose.cc
sec_structure_segments.cc
local_dist_test.cc
+ superpose_frames.cc
)
set(MOL_ALG_DEPS mol)
diff --git a/modules/mol/alg/src/superpose_frames.cc b/modules/mol/alg/src/superpose_frames.cc
new file mode 100644
index 0000000000000000000000000000000000000000..560293c166dcec0f31986493dcea0518d1e1929f
--- /dev/null
+++ b/modules/mol/alg/src/superpose_frames.cc
@@ -0,0 +1,170 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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 <Eigen/Core>
+#include <Eigen/Array>
+#include <Eigen/SVD>
+#include <Eigen/LU>
+#include <ost/mol/mol.hh>
+#include <ost/mol/alg/superpose_frames.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+
+typedef Eigen::Matrix<Real, 3, 3> EMat3;
+typedef Eigen::Matrix<Real, 1, 3> ECVec3;
+typedef Eigen::Matrix<Real, 3, 1> ERVec3;
+typedef Eigen::Matrix<Real, 3, Eigen::Dynamic> EMat3X;
+typedef Eigen::Matrix<Real, Eigen::Dynamic, 3> EMatX3;
+
+
+inline geom::Vec3 rvec_to_gvec(const ERVec3 &vec) {
+ return *reinterpret_cast<const geom::Vec3*>(&vec);
+}
+
+inline geom::Vec3 cvec_to_gvec(const ECVec3 &vec) {
+ return *reinterpret_cast<const geom::Vec3*>(&vec);
+}
+
+inline geom::Mat3 emat_to_gmat(const EMat3 &mat)
+{
+ return *reinterpret_cast<const geom::Mat3*>(&mat);
+}
+
+inline ERVec3 gvec_to_rvec(const geom::Vec3 &vec)
+{
+ return *reinterpret_cast<const ERVec3*>(&vec);
+}
+
+inline ECVec3 gvec_to_cvec(const geom::Vec3 &vec)
+{
+ return *reinterpret_cast<const ECVec3*>(&vec);
+}
+
+inline EMat3X row_sub(const EMat3X& m, const ERVec3& s)
+{
+ EMat3X r(m);
+ for (int col=0; col<m.cols();++col) {
+ r.col(col)-=s;
+ }
+ return r;
+}
+
+inline EMatX3 col_sub(const EMatX3& m, const ECVec3& s)
+{
+ EMatX3 r(m);
+ for (int row=0; row<m.rows();++row) {
+ r.row(row)-=s;
+ }
+ return r;
+}
+
+CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
+ int begin, int end, int ref)
+{
+ int real_end=end==-1 ? cg.GetFrameCount() : end;
+ CoordFramePtr ref_frame;
+ std::vector<unsigned long> indices;
+ if (!sel.IsValid()) {
+ indices.reserve(cg.GetAtomCount());
+ for (size_t i=0;i<cg.GetAtomCount(); ++i) {
+ indices.push_back(i);
+ }
+ } else {
+ AtomViewList atoms=sel.GetAtomList();
+ indices.reserve(atoms.size());
+ for (AtomViewList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ indices.push_back(i->GetIndex());
+ }
+ }
+ CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+ EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
+ EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
+ EMat3X frame_centered=EMat3X::Zero(3, indices.size());
+ EMat3X frame_mat=EMat3X::Zero(3, indices.size());
+ ECVec3 ref_center;
+ ERVec3 frame_center;
+ if (ref!=-1) {
+ ref_frame=cg.GetFrame(ref);
+ for (size_t i=0; i<indices.size(); ++i) {
+ ref_mat.row(i)=gvec_to_cvec((*ref_frame)[indices[i]]);
+ }
+ ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+ ref_centered=col_sub(ref_mat, ref_center);
+ } else {
+ ref_frame=cg.GetFrame(begin);
+ }
+ for (int i=begin; i<real_end; ++i) {
+ if (ref==i || (ref==-1 && i==begin)) {
+ superposed.AddFrame(*ref_frame);
+ ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
+ continue;
+ }
+ if (ref==-1) {
+ for (size_t j=0; j<indices.size(); ++j) {
+ ref_mat.row(j)=gvec_to_cvec((*ref_frame)[indices[j]]);
+ }
+ ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+ ref_centered=col_sub(ref_mat, ref_center);
+ }
+ CoordFramePtr frame=cg.GetFrame(i);
+ for (size_t j=0; j<indices.size(); ++j) {
+ frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
+ }
+ std::vector<geom::Vec3> frame_data=*frame;
+ frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
+ frame_centered=row_sub(frame_mat, frame_center);
+ //single value decomposition
+ Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
+ EMat3 matrixVT=svd.matrixV().transpose();
+ //determine rotation
+ Real detv=matrixVT.determinant();
+ Real dett=svd.matrixU().determinant();
+ Real det=detv*dett;
+ EMat3 e_rot;
+ if (det<0) {
+ EMat3 tmat=EMat3::Identity();
+ tmat(2,2)=-1;
+ e_rot=(svd.matrixU()*tmat)*matrixVT;
+ }else{
+ e_rot=svd.matrixU()*matrixVT;
+ }
+ // prepare rmsd calculation
+ geom::Vec3 shift=rvec_to_gvec(ref_center);
+ geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
+ geom::Mat3 rot=emat_to_gmat(e_rot);
+ geom::Mat4 mat4_com, mat4_rot, mat4_shift;
+ mat4_rot.PasteRotation(rot);
+ mat4_shift.PasteTranslation(shift);
+ mat4_com.PasteTranslation(com_vec);
+ geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
+ for (std::vector<geom::Vec3>::iterator c=frame_data.begin(),
+ e2=frame_data.end(); c!=e2; ++c) {
+ *c=geom::Vec3(tf*geom::Vec4(*c));
+ }
+ superposed.AddFrame(frame_data);
+ if (ref==-1) {
+ ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
+ }
+ }
+ return superposed;
+}
+
+}}}
diff --git a/modules/mol/alg/src/superpose_frames.hh b/modules/mol/alg/src/superpose_frames.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9416c5a9236d9bcd090eccedc24b8655ee15a2a5
--- /dev/null
+++ b/modules/mol/alg/src/superpose_frames.hh
@@ -0,0 +1,34 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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
+//------------------------------------------------------------------------------
+#ifndef OST_MOL_ALG_SUPERPOSE_FRAMES_HH
+#define OST_MOL_ALG_SUPERPOSE_FRAMES_HH
+
+#include <ost/mol/coord_group.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/alg/module_config.hh>
+namespace ost { namespace mol { namespace alg {
+
+/// \brief returns a superposed version of coord group
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle cg,
+ EntityView sel,
+ int begin=0, int end=-1,
+ int ref=-1);
+}}}
+
+#endif
diff --git a/modules/mol/base/src/coord_group.hh b/modules/mol/base/src/coord_group.hh
index 38fa81d2e5089989c1ecc01300784c694f485217..70844a26b017091ca58cde4b7570fb970d85571b 100644
--- a/modules/mol/base/src/coord_group.hh
+++ b/modules/mol/base/src/coord_group.hh
@@ -33,9 +33,6 @@ namespace ost { namespace mol {
/// \brief coordinate group, for trajectories and such
class DLLEXPORT_OST_MOL CoordGroupHandle {
- friend DLLEXPORT_OST_MOL
- CoordGroupHandle CreateCoordGroup(const std::vector<AtomHandle>&);
-
public:
/// \brief create empty, invalid handle
CoordGroupHandle();
@@ -84,10 +81,10 @@ public:
AtomHandleList GetAtomList() const;
CoordFramePtr GetFrame(uint frame) const;
-private:
- void CheckValidity() const;
CoordGroupHandle(CoordSourcePtr source);
+private:
+ void CheckValidity() const;
CoordSourcePtr source_;
};
diff --git a/modules/mol/base/src/impl/entity_impl.cc b/modules/mol/base/src/impl/entity_impl.cc
index 120c60e475534c1ce4a1080b9fb8364f17705fa2..33d9bef0552fee83ec02d05676d2ca4a6b1526a4 100644
--- a/modules/mol/base/src/impl/entity_impl.cc
+++ b/modules/mol/base/src/impl/entity_impl.cc
@@ -707,7 +707,7 @@ void EntityImpl::UpdateFromXCS()
TraceDirectionality();
}
- if(Logger::Instance().GetLogLevel()>Logger::DEBUG) {
+ if(Logger::Instance().GetVerbosityLevel()>Logger::DEBUG) {
LOG_TRACE("dumping directionality");
for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) {
LOG_TRACE(" " << it->second << ":");
diff --git a/scripts/init_cl.py b/scripts/init_cl.py
index 971b6ea1eff8ae1277a0e2fd93ddf1ac69d27407..fba11d8adb9fa5a6f6b9baca7905b9de55c028e5 100644
--- a/scripts/init_cl.py
+++ b/scripts/init_cl.py
@@ -16,10 +16,9 @@ class OstOptionParser(optparse.OptionParser):
optparse.OptionParser.__init__(self, **kwargs)
def exit(self, status_code, error_message):
print error_message,
- QtGui.QApplication.instance().exit()
sys.exit(-1)
-parser=OstOptionParser(usage=usage,conflict_handler="resolve")
+parser=OstOptionParser(usage=usage,conflict_handler="resolve", prog='ost''')
parser.add_option("-i", "--interactive", action="callback", callback=interactive_flag, help="start interpreter interactively (must be first parameter, ignored otherwise)")
parser.add_option("-h", "--help", action="callback", callback=show_help, help="show this help message and exit")
parser.add_option("-v", "--verbosity_level", action="store", type="int", dest="vlevel", default=2, help="sets the verbosity level [default: %default]")
diff --git a/scripts/ost.in b/scripts/ost.in
index 17b9ba416ce2dbedad14dd288e48128ae1ec2ea9..16dc6738783ce6ef78dcf6c170feb30125a47261 100755
--- a/scripts/ost.in
+++ b/scripts/ost.in
@@ -70,8 +70,10 @@ else
if [ $1 == "-i" ] ;then
interactive="-i"
fi
+ if [ $1 == "-v" ] ;then
+ interactive="-i"
+ fi
fi
-
# finally pass control to python instance
IFS="#"
exec $pyexec $interactive "$DNG_ROOT/@LIBDIR@/openstructure/init_cl.py" $opts