diff --git a/CMakeLists.txt b/CMakeLists.txt
index 00054ead6a5fe2a5bdea7ba0526d34be02046c74..cfebaa1b5068363a4433c737cb9947ccb1e5d320 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -22,6 +22,8 @@ option(ENABLE_GFX "whether graphics support should be enabled"
ON)
option(ENABLE_IMG "whether the image processing module should be compiled"
ON)
+option(USE_NUMPY "whether numpy support is added"
+ OFF)
option(USE_DOUBLE_PRECISION "whether to compile in double precision"
OFF)
option(ENABLE_SPNAV "whether 3DConnexion devices should be supported"
@@ -87,6 +89,11 @@ if (USE_SHADER)
else()
set(_SHADER OFF)
endif()
+if (USE_NUMPY)
+ set(_NUMPY ON)
+else()
+ set(_NUMPY OFF)
+endif()
if (COMPILE_TMTOOLS)
set(_TM_TOOLS ON)
else()
@@ -143,6 +150,10 @@ find_package(PNG REQUIRED)
find_package(Eigen 2.0.0 REQUIRED)
find_package(Python 2.4 REQUIRED)
+if(USE_NUMPY)
+ find_package(Numpy REQUIRED)
+endif()
+
if (ENABLE_IMG)
find_package(FFTW REQUIRED)
find_package(TIFF REQUIRED)
@@ -177,6 +188,10 @@ include_directories(${Boost_INCLUDE_DIRS}
${SPNAV_INCLUDE_DIR}
${PNG_INCLUDE_DIR}
)
+if(USE_NUMPY)
+include_directories(${PYTHON_NUMPY_INCLUDE_DIR})
+endif()
+
if (UNIX)
SET(CMAKE_SKIP_BUILD_RPATH FALSE)
SET(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
@@ -209,6 +224,7 @@ message(STATUS
" OpenGL support (-DENABLE_GFX) : ${_OPENGL}\n"
" Image Processing support (-DENABLE_IMG) : ${_IMG}\n"
" Shader support (-DUSE_SHADER) : ${_SHADER}\n"
+ " Numpy support (-DUSE_NUMPY) : ${_NUMPY}\n"
" Optimize (-DOPTIMIZE) : ${_OPT}\n"
" Profiling support (-DPROFILE) : ${_PROFILE}\n"
" Double Precision (-DUSE_DOUBLE_PRECISION) : ${_DOUBLE_PREC}\n"
diff --git a/cmake_support/FindNumpy.cmake b/cmake_support/FindNumpy.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..6348b9c91a934e4c01edac011cc861adc9fd059a
--- /dev/null
+++ b/cmake_support/FindNumpy.cmake
@@ -0,0 +1,32 @@
+if (PYTHON_NUMPY_INCLUDE_DIR)
+ # in cache already
+ set (PYTHON_NUMPY_FIND_QUIETLY TRUE)
+endif (PYTHON_NUMPY_INCLUDE_DIR)
+
+#INCLUDE(FindPython)
+
+IF(PYTHON_BINARY)
+ EXEC_PROGRAM ("${PYTHON_BINARY}"
+ ARGS "-c 'import numpy; print numpy.get_include()'"
+ OUTPUT_VARIABLE PYTHON_NUMPY_INCLUDE_DIR
+ RETURN_VALUE PYTHON_NUMPY_NOT_FOUND)
+
+ if (PYTHON_NUMPY_INCLUDE_DIR)
+ set (PYTHON_NUMPY_FOUND TRUE)
+ set (PYTHON_NUMPY_INCLUDE_DIR ${PYTHON_NUMPY_INCLUDE_DIR} CACHE STRING "Numpy include path")
+ else (PYTHON_NUMPY_INCLUDE_DIR)
+ set(PYTHON_NUMPY_FOUND FALSE)
+ endif (PYTHON_NUMPY_INCLUDE_DIR)
+ENDIF(PYTHON_BINARY)
+
+if (PYTHON_NUMPY_FOUND)
+ if (NOT PYTHON_NUMPY_FIND_QUIETLY)
+ message (STATUS "Numpy headers found")
+ endif (NOT PYTHON_NUMPY_FIND_QUIETLY)
+else (PYTHON_NUMPY_FOUND)
+ if (PYTHON_NUMPY_FIND_REQUIRED)
+ message (FATAL_ERROR "Numpy headers missing")
+ endif (PYTHON_NUMPY_FIND_REQUIRED)
+endif (PYTHON_NUMPY_FOUND)
+
+MARK_AS_ADVANCED (PYTHON_NUMPY_INCLUDE_DIR)
\ No newline at end of file
diff --git a/deployment/linux/create_bundle.py b/deployment/linux/create_bundle.py
index cf0ce5f6103cb48a99632a732d6672c88895c070..dbb88bb70feaa46329922ba30e3efefaea0b89fb 100644
--- a/deployment/linux/create_bundle.py
+++ b/deployment/linux/create_bundle.py
@@ -8,11 +8,11 @@ import datetime
# custom parameters
-if len(sys.argv) < 2:
- print 'usage: create_bundle.py additional_label'
+if len(sys.argv) < 3:
+ print 'usage: create_bundle.py additional_label system_python_version'
sys.exit()
-system_python_version='python2.6'
+system_python_version=sys.argv[2]
system_python_bin='/usr/bin/'+system_python_version
system_python_libs='/usr/lib/'+system_python_version
second_system_python_libs_flag=True
diff --git a/modules/config/CMakeLists.txt b/modules/config/CMakeLists.txt
index ddf07651188d815e51ae669972981c4aeab469c1..65c9ca6bebf69a4def544c5c58dfc50be04ee98d 100644
--- a/modules/config/CMakeLists.txt
+++ b/modules/config/CMakeLists.txt
@@ -5,14 +5,20 @@ dllexport.hh
version.hh
)
-
set(animations_enabled 0)
+
if (USE_SHADER)
set(shader_support 1)
else()
set(shader_support 0)
endif()
+if (USE_NUMPY)
+ set(numpy_support 1)
+else()
+ set(numpy_support 0)
+endif()
+
if (PROFILE)
set(profiling_enabled 1)
else()
diff --git a/modules/config/config.hh.in b/modules/config/config.hh.in
index 0095327dd737efa71141f72144cf9d8b5bfed213..88653c1e8df5bc5f75bf99dd6554c5a25c5ffaa6 100644
--- a/modules/config/config.hh.in
+++ b/modules/config/config.hh.in
@@ -31,4 +31,5 @@
#define OST_STATIC_PROPERTY_WORKAROUND @static_props@
#define OST_SPNAV_ENABLED @spnav_enabled@
#define OST_FFT_USE_THREADS @fftw_use_threads@
+#define OST_NUMPY_SUPPORT_ENABLED @numpy_support@
#endif
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index ec9a6cb75d2b269bc229ef0a2311bdc9fe15b6da..db2b8bb7097dc76d0b74c7906a0ea4a63bda5929 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -205,7 +205,7 @@ public:
Vec3List(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
Vec3List(const Vec3List& rhs) : base_type(rhs) { }
-
+ Vec3List(const base_type& rhs) : base_type(rhs) { }
Vec3List& operator=(const Vec3List& rhs)
{
base_type::operator=(rhs);
diff --git a/modules/gui/doc/dialogs.rst b/modules/gui/doc/dialogs.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a2bc2a7f239fa296a00bc0658bd2760ea6380a1b
--- /dev/null
+++ b/modules/gui/doc/dialogs.rst
@@ -0,0 +1,25 @@
+Dialogs
+================================================================================
+.. currentmodule:: ost.gui
+
+OpenStructure provides several :class:`Dialogs` which can be called from its
+many menus in the top bar.
+
+
+.. image:: images/100208_Dialogs.png
+
+Scene Menu
+--------------------------------------------------------------------------------
+
+
+Superpose
+^^^^^^^^^
+
+ Structurally superpose two entities.
+
+ .. image:: images/100624_superpose_dialog.png
+
+ The ``Superpose`` entry in the ``Scene`` menu is only accessible if two
+ instances of the graphical :class:`~ost.gfx.Entity` are selected.
+
+ .. autoclass:: ost.gui.dng.superpositiondialog.SuperpositionDialog
diff --git a/modules/gui/doc/gui.rst b/modules/gui/doc/gui.rst
index 675bf6b139294abb9507955d9e4b1978ffb075b7..19a407ac087b2ea40dd5c3803a599a68e3730d4d 100644
--- a/modules/gui/doc/gui.rst
+++ b/modules/gui/doc/gui.rst
@@ -19,4 +19,4 @@ C++ Qt widgets as well as with PyQt widgets. Learn more about :doc:`python_cpp`.
tools
python_shell
sequence_viewer
-
+ dialogs
diff --git a/modules/gui/doc/images/100208_Dialogs.png b/modules/gui/doc/images/100208_Dialogs.png
new file mode 100644
index 0000000000000000000000000000000000000000..f42f8809a3848f01390db42a71dbcc96bba53a01
Binary files /dev/null and b/modules/gui/doc/images/100208_Dialogs.png differ
diff --git a/modules/gui/doc/images/100624_superpose_dialog.png b/modules/gui/doc/images/100624_superpose_dialog.png
new file mode 100644
index 0000000000000000000000000000000000000000..c01e657c23b0d60f6d08af94d632e95b0b74f6d3
Binary files /dev/null and b/modules/gui/doc/images/100624_superpose_dialog.png differ
diff --git a/modules/gui/pymod/CMakeLists.txt b/modules/gui/pymod/CMakeLists.txt
index dc632913113930804afd77130237c76850ad8edf..9a623a9a1877f68e6aa7ff1f39229489acee84bd 100644
--- a/modules/gui/pymod/CMakeLists.txt
+++ b/modules/gui/pymod/CMakeLists.txt
@@ -90,6 +90,7 @@ set(OST_GUI_PYMOD_MODULES
set(OST_GUI_PYMOD_DNG_MODULES
__init__.py
termuse.py
+ superpositiondialog.py
init.py
menu.py
)
diff --git a/modules/gui/pymod/__init__.py b/modules/gui/pymod/__init__.py
index c543655248a41ba9543326e84e15c7a4b6a249d8..a63b91820c28a4eeb9b03a7c0eb916b99fcc84b9 100644
--- a/modules/gui/pymod/__init__.py
+++ b/modules/gui/pymod/__init__.py
@@ -184,7 +184,27 @@ class OneOf:
if isinstance(node, cl):
return True
return False
-
+
+class TwoOf:
+ def __init__(self, *classes):
+ self.classes=classes
+ def __call__(self):
+ sel=SceneSelection.Instance()
+ act_count=sel.GetActiveNodeCount()
+ if act_count<2:
+ return False
+ found=0
+ for i in range(0, act_count):
+ node=sel.GetActiveNode(i)
+ for cl in self.classes:
+ if isinstance(node, cl):
+ found += 1
+ if found > 2:
+ return False
+ if found == 2:
+ return True
+ return False
+
class ManyOf:
def __init__(self, *classes):
self.classes=classes
diff --git a/modules/gui/pymod/dng/init.py b/modules/gui/pymod/dng/init.py
index 661ddadc40d4c2d02c2d90588a5bc3d625437908..d46af73d8e1a93678a576698e06bdd69b552b84c 100644
--- a/modules/gui/pymod/dng/init.py
+++ b/modules/gui/pymod/dng/init.py
@@ -21,6 +21,7 @@ from ost.gui.init_spacenav import _InitSpaceNav
from ost.gui.init_context_menu import _InitContextMenu
from ost.gui.init_splash import _InitSplash
from ost.gui.dng import termuse
+from ost.gui.dng import superpositiondialog
import ost.gui.dng.menu
from PyQt4.QtGui import *
def _my_exit(code):
diff --git a/modules/gui/pymod/dng/menu.py b/modules/gui/pymod/dng/menu.py
index c2b672c2b28e7a11a9855eec4aff08c722845016..de943d3225a2c25ea3f6a42b0aefbb6969e3f9ad 100644
--- a/modules/gui/pymod/dng/menu.py
+++ b/modules/gui/pymod/dng/menu.py
@@ -1,9 +1,12 @@
from PyQt4.QtCore import *
from PyQt4.QtGui import *
-from ost import gui, gfx, io
+from ost import *
+from ost import gui
from ost.gui.scene.loader_manager_widget import LoaderManagerWidget
from ost.gui.init_splash import _InitSplash
from ost.gui.dng import termuse
+from ost.gui.dng import superpositiondialog
+
import sys
class FileMenu(QMenu):
def __init__(self, parent=None):
@@ -134,6 +137,8 @@ class SceneMenu(QMenu):
enabled=gui.ManyOf(gfx.GfxObj))
gui.AddMenuAction(self, 'Fit To Screen', self._FitToScreen,
enabled=gui.OneOf(gfx.Entity))
+ gui.AddMenuAction(self, 'Superpose', self._SuperposeDialog,
+ enabled=gui.TwoOf(gfx.Entity))
gui.AddMenuAction(self, 'Save Snapshot', self._ExportScene)
gui.AddMenuAction(self, 'Scene Clipping', self._ClipScene, shortcut='Ctrl+Shift+C')
@@ -169,7 +174,31 @@ class SceneMenu(QMenu):
def _FitToScreen(self):
sel=gui.SceneSelection.Instance()
gfx.FitToScreen(sel.GetActiveNode(0))
-
+
+ def _SuperposeDialog(self):
+ sel=gui.SceneSelection.Instance()
+ act_count=sel.GetActiveNodeCount()
+ # we now that there have to be 2 gfx.Entities, because of using TwoOf to
+ # enable menu entry!
+ i = 0;
+ gfx_ent_1 = sel.GetActiveNode(i)
+ while not isinstance(gfx_ent_1, gfx.Entity):
+ i += 1
+ gfx_ent_1 = sel.GetActiveNode(i)
+ i += 1
+ gfx_ent_2 = sel.GetActiveNode(i)
+ while not isinstance(gfx_ent_2, gfx.Entity):
+ i += 1
+ gfx_ent_2 = sel.GetActiveNode(i)
+ sd = superpositiondialog.SuperpositionDialog(gfx_ent_1, gfx_ent_2)
+ if sd.reference == 0:
+ gfx_ent_1.UpdatePositions()
+ gfx.Scene().CenterOn(gfx_ent_1)
+ else:
+ gfx_ent_2.UpdatePositions()
+ gfx.Scene().CenterOn(gfx_ent_2)
+ LogScript('RMSD: %.3f'%sd.rmsd)
+
class WindowMenu(QMenu):
def __init__(self, parent=None):
QMenu.__init__(self, parent)
diff --git a/modules/gui/pymod/dng/superpositiondialog.py b/modules/gui/pymod/dng/superpositiondialog.py
new file mode 100644
index 0000000000000000000000000000000000000000..f426bc3307ac8af3412fb3b373e322a96f411662
--- /dev/null
+++ b/modules/gui/pymod/dng/superpositiondialog.py
@@ -0,0 +1,282 @@
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+#
+# Authors: Stefan Bienert
+#
+from PyQt4.QtCore import *
+from PyQt4.QtGui import *
+from ost.mol.alg import Superpose
+from ost import mol
+
+class ChainComboBox(QComboBox):
+ def __init__(self, ent, parent=None):
+ # class variables
+ self.all_chains = 'All'
+ QComboBox.__init__(self, parent)
+ self.entity = ent
+ self.addItem(self.all_chains)
+ for chain in self.entity.chains:
+ self.addItem(chain.name)
+ if self.count()>0:
+ self.setCurrentIndex(0)
+
+ def SetItems(self, ent):
+ self.clear()
+ self.entity = ent
+ self.addItem(self.all_chains)
+ for chain in self.entity.chains:
+ self.addItem(chain.name)
+ if self.count()>0:
+ self.setCurrentIndex(0)
+ return
+
+ def _GetSelectedChain(self):
+ if self.currentIndex() == -1:
+ return mol.EntityHandle()
+ elif self.currentText() == self.all_chains:
+ return self.entity
+ return self.entity.Select('cname=' + str(self.currentText()))
+
+ def _SetSelectedChain(self, chain):
+ if hasattr(chain, 'name'):
+ name = chain.name
+ else:
+ name = str(chain)
+ for i in range(self.count()):
+ if self.itemText(i) == name:
+ self.setCurrentIndex(i)
+ break
+ selected_chain = property(_GetSelectedChain, _SetSelectedChain)
+
+class SuperpositionDialog(QDialog):
+ """
+ Provides a graphical user interface to structurally superpose two entities.
+ Uses function :func:`~ost.mol.alg.Superpose`. The RMSD of two superposed
+ molecules will be stored in attribute ``rmsd``. An index for the selected
+ reference molecule will be stored in attribute ``reference``.
+
+ :param ent_one: The first entity
+ :type ent_one: :class:`~ost.mol.EntityView`, :class:`~ost.mol.EntityHandle`
+ or :class:`~ost.gfx.Entity`
+ :param ent_two: The second entity
+ :type ent_two: :class:`~ost.mol.EntityView`, :class:`~ost.mol.EntityHandle`
+ or :class:`~ost.gfx.Entity`
+
+ **Example Usage:**
+
+ .. code-block:: python
+
+ e1=io.LoadPDB('examples/code_fragments/entity/pdb1ake.ent')
+ e2=io.LoadPDB('examples/code_fragments/entity/pdb4ake.ent')
+
+ sd = superpositiondialog.SuperpositionDialog(e1, e2)
+
+ g1=gfx.Entity('G1', e1)
+ g2=gfx.Entity('G2', e2)
+ scene.Add(g1)
+ scene.Add(g2)
+
+ if sd.reference == 0:
+ scene.CenterOn(g1)
+ else:
+ scene.CenterOn(g2)
+
+ LogScript('RMSD: %.3f'%sd.rmsd)
+ """
+
+ def __init__(self, ent_one, ent_two, parent=None):
+ # class variables
+ self.rmsd = 0.0
+ self._mmethod_dict = {'number': 'number',
+ 'index': 'index',
+ 'local alignment': 'local-aln',
+ 'global alignment': 'global-aln'}
+ QDialog.__init__(self, parent)
+ self.setWindowTitle('Superpose structures')
+ if not isinstance(ent_one, mol.EntityHandle) and \
+ not isinstance(ent_one, mol.EntityView):
+ n_one = ent_one.GetName()
+ self.ent_one = ent_one.GetView()
+ else:
+ if isinstance(ent_one, mol.EntityHandle):
+ n_one = ent_one.GetName()
+ elif isinstance(ent_one, mol.EntityView):
+ n_one = ent_one.GetHandle().GetName()
+ self.ent_one = ent_one
+ if len(n_one) == 0:
+ n_one = '1'
+ if not isinstance(ent_two, mol.EntityHandle) and \
+ not isinstance(ent_two, mol.EntityView):
+ n_two = ent_two.GetName()
+ self.ent_two = ent_two.GetView()
+ else:
+ if isinstance(ent_two, mol.EntityHandle):
+ n_two = ent_two.GetName()
+ elif isinstance(ent_two, mol.EntityView):
+ n_two = ent_two.GetHandle().GetName()
+ self.ent_two = ent_two
+ if len(n_two) == 0:
+ n_two = '2'
+ if n_one == n_two:
+ n_one = n_one + ' 1'
+ n_two = n_two + ' 2'
+ layout = QGridLayout(self)
+ # select reference molecule
+ self.reference = 0;
+ self._reference = self._ReferenceSelection(n_one, n_two)
+ grow = 0
+ layout.addWidget(QLabel("reference"), grow, 0)
+ layout.addWidget(self._reference, grow, 1)
+ grow += 1
+ # chains
+ self._chain_one = ChainComboBox(self.ent_one, self)
+ self._chain_two = ChainComboBox(self.ent_two, self)
+ layout.addWidget(QLabel("reference chain"), grow, 0)
+ layout.addWidget(self._chain_one, grow, 1)
+ grow += 1
+ layout.addWidget(QLabel("chain"), grow, 0)
+ layout.addWidget(self._chain_two, grow, 1)
+ grow += 1
+ # link chain and reference selection
+ QObject.connect(self._reference,
+ SIGNAL('currentIndexChanged(int)'),
+ self._ChangeChainSelection)
+ # match methods
+ self._methods = self._MatchMethods()
+ layout.addWidget(QLabel('match residues by'), grow, 0)
+ grow += 1
+ layout.addWidget(self._methods)
+ # atoms
+ self._atoms = self._FetchAtoms(self._methods.size(),
+ self.ent_one,
+ self.ent_two)
+ self._atmselectbx, self._atmselectgrp = self._AtomSelectionBox()
+ layout.addWidget(self._atmselectbx, grow, 1)
+ grow += 1
+ # buttons
+ ok_button = QPushButton("Superpose")
+ QObject.connect(ok_button, SIGNAL('clicked()'), self.accept)
+ cancel_button = QPushButton("Cancel")
+ hbox_layout = QHBoxLayout()
+ hbox_layout.addStretch(1)
+ layout.addLayout(hbox_layout, grow, 0, 1, 2)
+ grow += 1
+ QObject.connect(cancel_button, SIGNAL('clicked()'), self.reject)
+ QObject.connect(self, SIGNAL('accepted()'), self._Superpose)
+ hbox_layout.addWidget(cancel_button, 0)
+ hbox_layout.addWidget(ok_button, 0)
+ ok_button.setDefault(True)
+ self.exec_()
+
+ def _Superpose(self):
+ view_one = self._chain_one.selected_chain
+ view_two = self._chain_two.selected_chain
+ atoms = self._GetAtomSelection()
+ sp = Superpose(view_one, view_two,
+ self._mmethod_dict[str(self._methods.currentText())],
+ atoms)
+ self.rmsd = sp.rmsd
+
+ def _toggle_atoms(self, checked):
+ if checked:
+ self._atoms.setEnabled(True)
+ else:
+ self._atoms.setEnabled(False)
+
+ def _AtomSelectionBox(self):
+ bt1 = QRadioButton('All')
+ bt2 = QRadioButton('Backbone')
+ bt3 = QRadioButton('CA')
+ self.cstmbtntxt = 'Custom'
+ custom_rbutton = QRadioButton(self.cstmbtntxt)
+ group = QButtonGroup()
+ group.addButton(bt1)
+ group.addButton(bt2)
+ group.addButton(bt3)
+ group.addButton(custom_rbutton)
+ bt1.setChecked(True)
+ vbox_layout = QVBoxLayout()
+ vbox_layout.addWidget(bt1)
+ vbox_layout.addWidget(bt2)
+ vbox_layout.addWidget(bt3)
+ vbox_layout.addWidget(custom_rbutton)
+ vbox_layout.addWidget(self._atoms)
+ QObject.connect(custom_rbutton, SIGNAL('toggled(bool)'), self._toggle_atoms)
+ box = QGroupBox("atom selection")
+ box.setLayout(vbox_layout)
+ return box, group
+
+ def _GetAtomSelection(self):
+ checkedbtn = self._atmselectgrp.checkedButton()
+ if str(checkedbtn.text()) != self.cstmbtntxt:
+ return str(checkedbtn.text())
+ slctn_model = self._atoms.selectionModel()
+ dt_model = slctn_model.model()
+ atms = list()
+ for idx in slctn_model.selectedRows():
+ slctn = dt_model.data(idx, Qt.DisplayRole).toString()
+ atms.append(str(slctn))
+ return atms
+
+ def _FetchAtoms(self, dim, ent_a, ent_b):
+ # fetch list of atoms: only those which are in both entities are considered
+ atm_dict = {}
+ for atm in ent_a.GetAtomList():
+ atm_dict[atm.name] = 0
+ for atm in ent_b.GetAtomList():
+ if atm.name in atm_dict:
+ atm_dict[atm.name] = 1
+ atmlst = QStringList()
+ for atm in sorted(atm_dict.keys()):
+ if atm_dict[atm]:
+ atmlst.append(atm)
+ elems = QStringListModel(atmlst)
+ atoms = QListView(self)
+ dim.setHeight(3*dim.height())
+ atoms.setFixedSize(dim)
+ atoms.setModel(elems)
+ atoms.setSelectionMode(QAbstractItemView.MultiSelection)
+ atoms.setEditTriggers(QAbstractItemView.NoEditTriggers)
+ atoms.setEnabled(False)
+ return atoms
+
+ def _ReferenceSelection(self, name_a, name_b):
+ cbox = QComboBox()
+ cbox.addItem(name_a)
+ cbox.addItem(name_b)
+ if cbox.count() > 0:
+ cbox.setCurrentIndex(0)
+ return cbox
+
+ def _ChangeChainSelection(self, index):
+ if index == 0:
+ self._chain_one.SetItems(self.ent_one)
+ self._chain_two.SetItems(self.ent_two)
+ self.reference = 0;
+ elif index == 1:
+ self._chain_one.SetItems(self.ent_two)
+ self._chain_two.SetItems(self.ent_one)
+ self.reference = 1;
+ return
+
+ def _MatchMethods(self):
+ methods=QComboBox(self)
+ for method in sorted(self._mmethod_dict):
+ methods.addItem(method)
+ return methods
diff --git a/modules/img/alg/pymod/wrap_alg.cc b/modules/img/alg/pymod/wrap_alg.cc
index f0faf8c13e1b73f1b2d0266a28978c4cb1c6ba37..c141fe11310076374194677f601e98a1929eb446 100644
--- a/modules/img/alg/pymod/wrap_alg.cc
+++ b/modules/img/alg/pymod/wrap_alg.cc
@@ -187,12 +187,12 @@ BOOST_PYTHON_MODULE(_ost_img_alg)
.def(self_ns::str(self))
;
class_<StatMinMax, bases<NonModAlgorithm> >("StatMinMax", init<>() )
- .def("GetMinimum",&Stat::GetMinimum)
- .def("GetMinimumPosition",&Stat::GetMinimumPosition)
- .def("SetMinimum",&Stat::SetMinimum)
- .def("GetMaximum",&Stat::GetMaximum)
- .def("GetMaximumPosition",&Stat::GetMaximumPosition)
- .def("SetMaximum",&Stat::SetMaximum)
+ .def("GetMinimum",&StatMinMax::GetMinimum)
+ .def("GetMinimumPosition",&StatMinMax::GetMinimumPosition)
+ .def("SetMinimum",&StatMinMax::SetMinimum)
+ .def("GetMaximum",&StatMinMax::GetMaximum)
+ .def("GetMaximumPosition",&StatMinMax::GetMaximumPosition)
+ .def("SetMaximum",&StatMinMax::SetMaximum)
.def(self_ns::str(self))
;
diff --git a/modules/img/alg/tests/test_stat.cc b/modules/img/alg/tests/test_stat.cc
index b34689ec9e9a802feb3750234d63a6b613cbd619..6f74a044a0d8f600ae95debcac1984581d25e2d7 100644
--- a/modules/img/alg/tests/test_stat.cc
+++ b/modules/img/alg/tests/test_stat.cc
@@ -49,14 +49,8 @@ void test() {
Stat stat;
im.Apply(stat);
-
- std::ostringstream msg;
-
- msg << "expected 5.0 as mean but found " << stat.GetMean();
- BOOST_CHECK_MESSAGE(stat.GetMean()==5.0,msg.str());
- msg.str("");
- msg << "expected 2.73861 as stdev but found " << stat.GetStandardDeviation();
- BOOST_CHECK_MESSAGE((stat.GetStandardDeviation()-2.73861)<1e-5,msg.str());
+ BOOST_CHECK_CLOSE(stat.GetMean(),5.0,0.0001);
+ BOOST_CHECK_CLOSE(stat.GetStandardDeviation(),2.58198889747,0.0001);
}
diff --git a/modules/index.rst b/modules/index.rst
index 8cc338a82b9480f516d927e48b3894b80b61e8d0..afc2142cfb81f0b06d9c83755207b29523132d7a 100644
--- a/modules/index.rst
+++ b/modules/index.rst
@@ -41,7 +41,9 @@ For Starters
Molecules
--------------------------------------------------------------------------------
-**Overview**: :doc:`molecules intro <intro-01>` | :doc:`mol overview <mol/base/mol>` | :doc:`graphical entity<gfx/entity>` | :doc:`entity <mol/base/entity>` | :doc:`queries <mol/base/query>`
+**Overview**: :doc:`molecules intro <intro-01>` | :doc:`mol overview <mol/base/mol>` | :doc:`graphical entity<gfx/entity>` | :doc:`entity <mol/base/entity>` | :doc:`queries <mol/base/query>` | :doc:`algorithms <mol/alg/molalg>`
+
+**Trajectories**: :doc:`basics <mol/base/traj>` | :ref:`analysis <traj-analysis>`
**Input/Output**: :ref:`loading and saving molecules <mol-io>`
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index 3fdb02f55727260e8679c0ac30e922d10e8d78c6..e5f48a93d4486e2de8c358ebb310dfb589f3d9e6 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/modules/io/src/mol/dcd_io.cc
@@ -184,7 +184,7 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
}
-mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
+mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted!
const String& trj_fn,
unsigned int stride)
{
@@ -196,10 +196,6 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
}
Profile profile_load("LoadCHARMMTraj");
- mol::AtomHandleList alist(alist2);
- std::sort(alist.begin(),alist.end(),less_index);
-
-
DCDHeader header;
bool swap_flag=false, skip_flag=false, gap_flag=false;
read_dcd_header(istream, header, swap_flag, skip_flag, gap_flag);
diff --git a/modules/io/src/mol/sdf_writer.cc b/modules/io/src/mol/sdf_writer.cc
index 9f3154e9b99f91a47354a7dcb2f9ccf31cd1b758..40c0660153853bc5edbfeb8e94384194a8daf63b 100644
--- a/modules/io/src/mol/sdf_writer.cc
+++ b/modules/io/src/mol/sdf_writer.cc
@@ -27,6 +27,7 @@
#include <ost/mol/chain_view.hh>
#include <ost/mol/bond_handle.hh>
#include <boost/regex.hpp>
+#include <boost/bind.hpp>
namespace ost { namespace io {
@@ -38,7 +39,7 @@ namespace {
public:
SDFAtomWriter(std::ostream& ostream, std::map<long, int>& atom_indices)
: ostr_(ostream), atom_indices_(atom_indices), counter_(0) {
- atom_indices_.clear();
+ atom_indices_.clear();
}
private:
public:
@@ -60,23 +61,61 @@ namespace {
class SDFBondWriter : public mol::EntityViewVisitor {
public:
- SDFBondWriter(std::ostream& ostream, std::map<long, int>& atom_indices)
+ SDFBondWriter(std::ostream& ostream,
+ const std::map<long, int>& atom_indices)
: ostr_(ostream), atom_indices_(atom_indices), counter_(0) {
}
private:
+ // compare two atoms according to their indices (used for sorting)
+ bool CompareAtomIdx(const mol::AtomView& first,
+ const mol::AtomView& second) {
+ std::map<long, int>::const_iterator aidx_first(
+ atom_indices_.find(first.GetHashCode()));
+ std::map<long, int>::const_iterator aidx_second(
+ atom_indices_.find(second.GetHashCode()));
+
+ if(aidx_first==atom_indices_.end() || aidx_second==atom_indices_.end()) {
+ throw IOException("Cannot write bond: atom idx not found for sorting");
+ }
+ return (aidx_first->second < aidx_second->second);
+ }
+
public:
virtual bool VisitAtom(const mol::AtomView& atom) {
- counter_++;
+ ++counter_; // current atom index
+
+ // get all neighboring atoms and sort them according to their atom index
mol::AtomViewList atoms = atom.GetBondPartners();
- mol::AtomViewList::iterator atom_iter = atoms.begin();
- for(; atom_iter != atoms.end(); ++atom_iter) {
- int atom_index = atom_indices_.find((*atom_iter).GetHashCode())->second;
- if(atom_index > counter_) {
- int type = 1;
- mol::BondHandle bond = atom.GetHandle().FindBondToAtom(atom_iter->GetHandle());
- if(bond.IsValid()) type = bond.GetBondOrder();
+ std::sort(atoms.begin(), atoms.end(), bind(&SDFBondWriter::CompareAtomIdx,
+ this, _1, _2));
+
+ // iterate all neighboring atoms and print bonds to all atoms with index
+ // larger than current atom index
+ for(mol::AtomViewList::iterator atom_iter = atoms.begin();
+ atom_iter != atoms.end(); ++atom_iter) {
+ std::map<long, int>::const_iterator aidx(
+ atom_indices_.find((*atom_iter).GetHashCode()));
+
+ // check if index was found
+ if(aidx==atom_indices_.end()) {
+ throw IOException("Cannot write bond between " +
+ atom.GetQualifiedName() + " and " +
+ atom_iter->GetQualifiedName() +
+ ": atom index not found");
+ }
+
+ // only print bonds to atoms with larger index than current index
+ if(aidx->second > counter_) {
+ mol::BondHandle bond(atom.GetHandle().FindBondToAtom(
+ atom_iter->GetHandle()));
+ if(!bond.IsValid()) {
+ throw IOException("Bond is invalid between " +
+ atom.GetQualifiedName() + " and " +
+ atom_iter->GetQualifiedName());
+ }
+ int type = bond.GetBondOrder();
ostr_ << format("%3i") % counter_
- << format("%3i") % atom_index
+ << format("%3i") % aidx->second
<< format("%3i") % type
<< " 0 0 0"
<< std::endl;
@@ -87,17 +126,17 @@ namespace {
private:
std::ostream& ostr_;
- std::map<long, int>& atom_indices_;
+ const std::map<long, int>& atom_indices_;
int counter_;
};
}
SDFWriter::SDFWriter(std::ostream& ostream)
- : outfile_(), ostr_(ostream), counter_(0) {
+ : outfile_(), ostr_(ostream), counter_(0), atom_indices_() {
}
SDFWriter::SDFWriter(const String& filename)
- : outfile_(filename.c_str()), ostr_(outfile_), counter_(0) {
+ : outfile_(filename.c_str()), ostr_(outfile_), counter_(0), atom_indices_() {
}
SDFWriter::SDFWriter(const boost::filesystem::path& filename):
@@ -106,7 +145,7 @@ SDFWriter::SDFWriter(const boost::filesystem::path& filename):
#else
outfile_(filename.file_string().c_str()),
#endif
- ostr_(outfile_), counter_(0) {
+ ostr_(outfile_), counter_(0), atom_indices_() {
}
void SDFWriter::Write(const mol::EntityView& ent) {
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 578d9d1d1a20c924416e5d6eeb66bc14d822ca2d..4730eec073af926ec47ef043930e75f1685515fb 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -82,3 +82,137 @@ The following function detects steric clashes in atomic structures. Two atoms ar
:returns: The filtered :class:`~ost.mol.EntityView`
+.. _traj-analysis:
+
+Trajectory Analysis
+--------------------------------------------------------------------------------
+
+This is a set of functions used for basic trajectory analysis such as extracting
+positions, distances, angles and RMSDs. The organization is such that most
+functions have their counterpart at the individual :class:`frame level
+<ost.mol.CoordFrame>` so that they can also be called on one frame instead
+of the whole trajectory.
+
+All these functions have a "stride" argument that defaults to stride=1, which is
+used to skip frames in the analysis.
+
+
+.. function:: AnalyzeAtomPos(traj, atom1, stride=1)
+
+ This function extracts the position of an atom from a trajectory. It returns
+ a vector containing the position of the atom for each analyzed frame.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+.. function:: AnalyzeCenterOfMassPos(traj, sele, stride=1)
+
+ This function extracts the position of the center-of-mass of a selection
+ (:class:`~ost.mol.EntityView`) from a trajectory and returns it as a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param sele: The selection from which the center of mass is computed
+ :type sele: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+.. function:: AnalyzeDistanceBetwAtoms(traj, atom1, atom2, stride=1)
+
+ This function extracts the distance between two atoms from a trajectory
+ and returns it as a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+.. function:: AnalyzeAngle(traj, atom1, atom2, atom3, stride=1)
+
+ This function extracts the angle between three atoms from a trajectory
+ and returns it as a vector. The second atom is taken as being the central
+ atom, so that the angle is between the vectors (atom1.pos-atom2.pos) and
+ (atom3.pos-atom2.pos).
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param atom3: The third :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+.. function:: AnalyzeDihedralAngle(traj, atom1, atom2, atom3, atom4, stride=1)
+
+ This function extracts the dihedral angle between four atoms from a trajectory
+ and returns it as a vector. The angle is between the planes containing the
+ first three and the last three atoms.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param atom1: The first :class:`~ost.mol.AtomHandle`.
+ :param atom2: The second :class:`~ost.mol.AtomHandle`.
+ :param atom3: The third :class:`~ost.mol.AtomHandle`.
+ :param atom4: The fourth :class:`~ost.mol.AtomHandle`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+.. function:: AnalyzeDistanceBetwCenterOfMass(traj, sele1, sele2, stride=1)
+
+ This function extracts the distance between the center-of-mass of two
+ selections (:class:`~ost.mol.EntityView`) from a trajectory and returns it as
+ a vector.
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param sele1: The selection from which the first center of mass is computed
+ :type sele1: :class:`~ost.mol.EntityView`.
+ :param sele2: The selection from which the second center of mass is computed
+ :type sele2: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+.. function:: AnalyzeRMSD(traj, reference_view, sele_view)
+
+ This function extracts the rmsd between two :class:`~ost.mol.EntityView` and
+ returns it as a vector. The views don't have to be from the same entity. The
+ reference positions are taken directly from the reference_view, evaluated only
+ once. The positions from the sele_view are evaluated for each frame.
+ If you want to compare to frame i of the trajectory t, first use
+ t.CopyFrame(i) for example:
+
+ .. code-block:: python
+
+ eh=io.LoadPDB(...)
+ t=io.LoadCHARMMTraj(eh,...)
+ sele=eh.Select(...)
+ t.CopyFrame(0)
+ mol.alg.AnalyzeRMSD(t,sele,sele)
+
+ :param traj: The trajectory to be analyzed.
+ :type traj: :class:`~ost.mol.CoordGroupHandle`
+ :param reference_view: The selection used as reference structure
+ :type reference_view: :class:`~ost.mol.EntityView`.
+ :param sele_view: The selection compared to the reference_view
+ :type sele_view: :class:`~ost.mol.EntityView`.
+ :param stride: Size of the increment of the frame's index between two
+ consecutive frames analyzed.
+
+
+
+
+
+
+
+
+
diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index 2fbf108cbe9133fcab98cd55b30fb99c74dabd2a..0ffb3f65183cb217c5fcfe6b1093a884ea42bca9 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -2,6 +2,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
wrap_mol_alg.cc
export_svd_superpose.cc
export_clash.cc
+ export_trajectory_analysis.cc
)
set(OST_MOL_ALG_PYMOD_MODULES
diff --git a/modules/mol/alg/pymod/export_superpose_frames.cc b/modules/mol/alg/pymod/export_superpose_frames.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e1d55652052388dbf7721f6d030fbe6ef0cbe79f
--- /dev/null
+++ b/modules/mol/alg/pymod/export_superpose_frames.cc
@@ -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
+//------------------------------------------------------------------------------
+
+/*
+ * Author Marco Biasini
+ */
+#include <boost/python.hpp>
+using namespace boost::python;
+#include <ost/mol/alg/superpose_frames.hh>
+
+using namespace ost;
+using namespace ost::mol::alg;
+
+void export_SuperposeFrames()
+{
+ def("SuperposeFrames", &SuperposeFrames);
+}
+
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..61e153bb3f060f61937bb75d8ebb4458d6edcb1b
--- /dev/null
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -0,0 +1,46 @@
+//------------------------------------------------------------------------------
+// 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/trajectory_analysis.hh>
+
+using namespace ost;
+using namespace ost::mol::alg;
+
+//"thin wrappers" for default parameters
+/*
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractapos, ExtractAtomPosition, 2, 3)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmpos, ExtractCMPosition, 2, 3)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractdist, ExtractDistance, 3, 4)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractang, ExtractAngle, 4, 5)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractdih, ExtractDihedral, 5, 6)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmdist, ExtractCMDistance, 3, 4)
+BOOST_PYTHON_FUNCTION_OVERLOADS(extractrmsd, ExtractRMSD, 3, 4)
+*/
+void export_TrajectoryAnalysis()
+{
+ def("AnalyzeAtomPos",&AnalyzeAtomPos, (arg("traj"), arg("atom"), arg("stride")=1));
+ def("AnalyzeCenterOfMassPos",&AnalyzeCenterOfMassPos, (arg("traj"), arg("selection"), arg("stride")=1));
+ def("AnalyzeDistanceBetwAtoms",&AnalyzeDistanceBetwAtoms, (arg("traj"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeAngle",&AnalyzeAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeDihedralAngle",&AnalyzeDihedralAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+ def("AnalyzeDistanceBetwCenterOfMass",&AnalyzeDistanceBetwCenterOfMass, (arg("traj"), arg("selection"), arg("selection"), arg("stride")=1));
+ def("AnalyzeRMSD",&AnalyzeRMSD, (arg("traj"), arg("reference_view"), arg("selection"), arg("stride")=1));
+}
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 82803c983ac9eed3ecc0364326310476be57ece9..48e19d84fc7b35e599454f809c07b27abcf1c461 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -26,7 +26,8 @@ using namespace boost::python;
using namespace ost;
void export_svdSuperPose();
-
+void export_TrajectoryAnalysis();
+void export_Clash();
#if OST_IMG_ENABLED
void export_entity_to_density();
#endif
@@ -42,6 +43,7 @@ mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterC
BOOST_PYTHON_MODULE(_ost_mol_alg)
{
export_svdSuperPose();
+ export_TrajectoryAnalysis();
#if OST_IMG_ENABLED
export_entity_to_density();
#endif
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index ae27e33e6f0160d5556d0028cc3d75fe32e16c62..1bd52057e93d59a5de0ec09edab397b5bc812261 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -7,6 +7,7 @@ set(OST_MOL_ALG_HEADERS
filter_clashes.hh
construct_cbeta.hh
clash_score.hh
+ trajectory_analysis.hh
)
set(OST_MOL_ALG_SOURCES
@@ -17,6 +18,7 @@ set(OST_MOL_ALG_SOURCES
superpose_frames.cc
filter_clashes.cc
construct_cbeta.cc
+ trajectory_analysis.cc
)
set(MOL_ALG_DEPS ost_mol ost_seq)
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index bfb435c2f2f06d97dba79180591d057af380fe62..0fbca1696e5bc3ef1196f6d2a4863a1b0d72dcef 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -132,6 +132,7 @@ EntityView FilterClashes(const EntityView& ent, Real tolerance,
filtered.AddAtom(atom);
}
}
+ continue;
}
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS);
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c8660fc99d1ae118daac98fe802c5fe8ad29564a
--- /dev/null
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -0,0 +1,157 @@
+//------------------------------------------------------------------------------
+// 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 <stdexcept>
+#include <ost/base.hh>
+#include <ost/geom/vec3.hh>
+#include <ost/base.hh>
+#include <ost/geom/geom.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/coord_group.hh>
+#include "trajectory_analysis.hh"
+
+
+namespace ost { namespace mol { namespace alg {
+
+geom::Vec3List AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
+//This function extracts the position of an atom from a trajectory and returns it as a vector of geom::Vec3
+{
+ CheckHandleValidity(traj);
+ geom::Vec3List pos;
+ pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ pos.push_back(frame->GetAtomPos(i1));
+ }
+ return pos;
+}
+
+std::vector<Real> AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ unsigned int stride)
+//This function extracts the distance between two atoms from a trajectory and returns it as a vector
+{
+ CheckHandleValidity(traj);
+ std::vector<Real> dist;
+ dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex();
+ int i2=a2.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ dist.push_back(frame->GetDistanceBetwAtoms(i1,i2));
+ }
+ return dist;
+}
+
+std::vector<Real> AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ const AtomHandle& a3, unsigned int stride)
+//This function extracts the angle between three atoms from a trajectory and returns it as a vector
+{
+ CheckHandleValidity(traj);
+ std::vector<Real> ang;
+ ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ ang.push_back(frame->GetAngle(i1,i2,i3));
+ }
+ return ang;
+}
+
+std::vector<Real> AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
+ const AtomHandle& a3, const AtomHandle& a4, unsigned int stride)
+//This function extracts the dihedral angle between four atoms from a trajectory and returns it as a vector
+{
+ CheckHandleValidity(traj);
+ std::vector<Real> ang;
+ ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex(),i4=a4.GetIndex();
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ ang.push_back(frame->GetDihedralAngle(i1,i2,i3,i4));
+ }
+ return ang;
+}
+
+geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride)
+//This function extracts the position of the CenterOfMass of a selection (entity view) from a trajectory
+//and returns it as a vector.
+ {
+ CheckHandleValidity(traj);
+ geom::Vec3List pos;
+ pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> indices;
+ std::vector<Real> masses;
+ GetIndicesAndMasses(Sele, indices, masses);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ pos.push_back(frame->GetCenterOfMassPos(indices,masses));
+ }
+ return pos;
+}
+
+std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& Sele1,
+ const EntityView& Sele2, unsigned int stride)
+//This function extracts the distance between the CenterOfMass of two selection (entity views) from a trajectory
+//and returns it as a vector.
+ {
+ CheckHandleValidity(traj);
+ std::vector<Real> dist;
+ dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> indices1,indices2;
+ std::vector<Real> masses1,masses2;
+ GetIndicesAndMasses(Sele1, indices1, masses1);
+ GetIndicesAndMasses(Sele2, indices2, masses2);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ dist.push_back(frame->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2));
+ }
+ return dist;
+}
+
+std::vector<Real> AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view,
+ const EntityView& sele_view, unsigned int stride)
+// This function extracts the rmsd between two entity views and returns it as a vector
+// The views don't have to be from the same entity
+// If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
+// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,Sele,Sele)
+ {
+ CheckHandleValidity(traj);
+ int count_ref=reference_view.GetAtomCount();
+ int count_sele=sele_view.GetAtomCount();
+ if (count_ref!=count_sele){
+ throw std::runtime_error("atom counts of the two views are not equal");
+ }
+ std::vector<Real> rmsd;
+ rmsd.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ std::vector<unsigned long> sele_indices;
+ std::vector<geom::Vec3> ref_pos;
+ GetIndices(reference_view, sele_indices);
+ GetPositions(sele_view, ref_pos);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ rmsd.push_back(frame->GetRMSD(ref_pos,sele_indices));
+ }
+ return rmsd;
+}
+
+}}} //ns
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d2947e63b52a499e84569695e3bfe9f7b6af51a2
--- /dev/null
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -0,0 +1,46 @@
+//------------------------------------------------------------------------------
+// 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_TRAJECTORY_ANALYSIS_HH
+#define OST_TRAJECTORY_ANALYSIS_HH
+
+#include <ost/mol/alg/module_config.hh>
+
+#include <ost/base.hh>
+#include <ost/geom/geom.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/coord_group.hh>
+
+
+namespace ost { namespace mol { namespace alg {
+
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
+ geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& sele,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& sele1, const EntityView& sele2,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4,unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view, const EntityView& sele,unsigned int stride=1);
+
+
+}}}//ns
+#endif
diff --git a/modules/mol/alg/tests/CMakeLists.txt b/modules/mol/alg/tests/CMakeLists.txt
index 09d1934f76b97f775c01381c741d1a9c33a530f5..1008a3578e384b05160bf5ea4a127d1be693b349 100644
--- a/modules/mol/alg/tests/CMakeLists.txt
+++ b/modules/mol/alg/tests/CMakeLists.txt
@@ -4,5 +4,5 @@ set(OST_MOL_ALG_UNIT_TESTS
test_convenient_superpose.py
)
-ost_unittest(mol_alg "${OST_MOL_ALG_UNIT_TESTS}")
+ost_unittest(MODULE mol_alg SOURCES "${OST_MOL_ALG_UNIT_TESTS}")
diff --git a/modules/mol/base/doc/entity.rst b/modules/mol/base/doc/entity.rst
index 4520df521aa86b3f51c8e61bd59a03864903a76b..1ee9e6bdae35cce7abfe51d84e0315a354bc25da 100644
--- a/modules/mol/base/doc/entity.rst
+++ b/modules/mol/base/doc/entity.rst
@@ -583,7 +583,18 @@ The Handle Classes
.. attribute:: pos
- The atom's position in global coordinates. Also available as :meth:`GetPos`.
+ The atom's position in global coordinates, transformed by the entity
+ transformation. Also available as :meth:`GetPos`.
+
+ Read-only.
+
+ :type: :class:`~ost.geom.Vec3`
+
+ .. attribute:: original_pos
+
+ The atom's untransformed position in global coordinates. Also available as
+ :meth:`GetOriginalPos`.
+
Read-only.
:type: :class:`~ost.geom.Vec3`
diff --git a/modules/mol/base/doc/mol.rst b/modules/mol/base/doc/mol.rst
index 1ac7e3f2691b0a28cd299f7176ffa43e8bb4451b..976c2c8aa1ced240e5c461fb7f46745c60d0b36a 100644
--- a/modules/mol/base/doc/mol.rst
+++ b/modules/mol/base/doc/mol.rst
@@ -12,4 +12,5 @@ The mol module implements data structures to work with molecular datasets. At it
entity
editors
query
- surface
\ No newline at end of file
+ surface
+ traj
\ No newline at end of file
diff --git a/modules/mol/base/doc/traj.rst b/modules/mol/base/doc/traj.rst
new file mode 100644
index 0000000000000000000000000000000000000000..41730915f3efd9c94fc6fe59598f01f9a4125812
--- /dev/null
+++ b/modules/mol/base/doc/traj.rst
@@ -0,0 +1,220 @@
+Trajectories
+================================================================================
+
+.. currentmodule:: ost.mol
+
+
+.. function:: CreateCoordGroup(atoms)
+
+ :param atoms: List of atoms. All atoms must be from the same entity.
+ :type atoms: :class:`AtomHandleList`
+
+ :rtype: :class:`CoordGroupHandle`
+ :returns: A coord group with zero frames
+
+.. class:: CoordGroupHandle
+
+ A collection of coordinate frames, e.g. an MD trajectory. Create with
+ :func:`CreateCoordGroup`.
+
+
+ .. attribute:: entity
+
+ The attached entity.
+
+ :type: :class:`EntityHandle`
+
+ .. attribute:: atoms
+
+ The atoms of this coord group. The order of atoms is the same as the
+ positions in the coord frames. All atoms are from the same entity.
+
+ :type: :class:`AtomHandleList`
+ .. method:: AddFrames(frames)
+
+ Combine two trajectories by appending the frames of the second to the first
+ :param frames: a valid coord group
+ :type frames: :class:`CoordGroupHandle`
+
+ .. method:: Capture()
+
+ Record the atom positions of the entity attached to this coord group in a
+ new coordinate frame. Note that the atom positions transformed by the entity
+ transform will be stored. Only available for mutable coord groups.
+
+ :see: :attr:`AtomHandle.pos`
+
+ .. method:: CaptureInto(frame_index)
+
+ Same as :meth:`Capture`, but doesn't create a new frame and stores the
+ coordinates directly into frame with index *frame_index*. Only available for
+ mutable coord groups.
+
+ :param frame_index: index of the frame
+ :type frame_index: int
+
+ .. method:: CopyFrame(frame_index)
+
+ Copies the coordinates of frame with index *frame_index* to the attached
+ entity.
+
+ :param frame_index: index of the frame
+ :type frame_index: int
+
+ .. method:: Filter(view)
+
+ Returns a new trajectory containing only coordinates of the atoms in view.
+ Useful to remove water and other solvent molecules from a trajectory to
+ save memory.
+
+ :param view: a valid entity view
+ :type view: :class:`EntityView`
+ :rtype: :class:`CoordGroupHandle`
+
+ .. method:: GetAtomCount()
+
+ Returns the number of atoms in the coord group
+
+ :rtype: int
+
+
+ .. method:: GetAtomList()
+
+ Returns the atoms of the coord group in the same order they appear in the
+ coordinate frames.
+
+ :rtype: :class:`AtomHandleList`
+
+ .. method:: GetAtomPos(frame_index, atom)
+
+ Get position of *atom* in frame with index *frame_index*.
+
+
+ :param frame_index: frame index
+ :type frame_index: int
+ :param atom: A valid atom
+ :type atom: :class:`AtomHandle`
+ :rtype: :class:`Vec3`
+
+ .. method:: GetEntity()
+
+ Returns the attached entity
+
+ :rtype: :class:`EntityHandle`
+
+ .. method:: GetFrameCount()
+
+ Returns the number of frames of this coord group
+
+ :rtype: int
+
+ .. method:: IsValid()
+
+ Whether this coord group is valid
+
+ :rtype: bool
+
+ .. method:: SetAtomPos(frame_index, atom, pos)
+
+ Set position of *atom* in frame with index *frame_index* to *pos*
+
+ :param frame_index: index of the frame
+ :type frame_index: int
+ :param atom: a valid atom
+ :type atom: :class:`AtomHandle`
+ :param pos: new position of the atom
+ :type pos: :class:`Vec3`
+
+ .. method:: SetFramePositions(frame_index, positions)
+
+ Set the frame positions of frame with index *frame_index*. Order and count
+ of positions must match :attr:`atoms`.
+
+ :param frame_index: index of frame
+ :type frame_index: int
+ :param positions: list of positions
+ :type positions: :class:`~ost.geom.Vec3List`
+
+
+.. class:: CoordFrame
+
+ A single frame of coordinates in a :class:`CoordGroupHandle`.
+
+ .. method:: GetAngle(atom1, atom2, atom3)
+
+ :param atom1: first atom
+ :type atom1: :class:`AtomHandle`
+ :param atom2: second (central) atom
+ :type atom2: :class:`AtomHandle`
+ :param atom3: third atom
+ :type atom3: :class:`AtomHandle`
+
+ :returns: the angle in radians between the 3 atoms
+ :rtype: float
+
+ .. method:: GetAtomPos(atom)
+
+ Returns the position of the atom in the coord frame
+
+ :param atom: A valid atom handle
+ :type atom: :class:`AtomHandle`
+
+
+
+ :rtype: :class:`Vec3`
+
+ .. method:: GetCenterOfMassPos(view)
+
+
+ :param view: A valid entity view
+ :type view: :class:`EntityView`
+ :rtype: :class:`Vec3`
+
+ .. method:: GetDihedralAngle(atom1, atom2, atom3, atom4)
+
+ Get dihedral angle of the four atoms.
+
+ :param atom1: First atom. Must be valid
+ :type atom1: :class:`AtomHandle`
+ :param atom2: Second atom. Must be valid
+ :type atom2: :class:`AtomHandle`
+ :param atom3: Third atom. Must be valid
+ :type atom3: :class:`AtomHandle`
+ :param atom3: Fourth atom. Must be valid
+ :type atom3: :class:`AtomHandle`
+
+ :rtype: float
+
+ .. method:: GetDistanceBetwAtoms(atom1, atom2)
+
+ Get distance in (Angstroem) between *atom1* and *atom2* in coord frame.
+
+ :param atom1: First atom. Must be valid
+ :type atom1: :class:`AtomHandle`
+ :param atom2: Second atom. Must be valid
+ :type atom2: :class:`AtomHandle`
+ :rtype: float
+
+ .. method:: GetDistanceBetwCenterOfMass(view1, view2)
+
+ Get distance between center of mass of the first selection and the second.
+
+ :param view1: First view. Must be valid
+ :type view1: :class:`EntityView`
+ :param view2: Second view. Must be valid
+ :type view2: :class:`EntityView`
+ :rtype: float
+
+ .. method:: GetRMSD(view1, view2)
+
+ Get RMSD between two views in the coord frame. The coordinates of the views
+ are taken as is without superposition. The two views must have the same
+ number of atoms. Atoms are matches as they appear in
+ :attr:`EntityView.atoms`.
+
+ :param view1: First view. Must be valid
+ :type view1: :class:`EntityView`
+ :param view2: Second view. Must be valid
+ :type view2: :class:`EntityView`
+ :rtype: float
+
diff --git a/modules/mol/base/pymod/CMakeLists.txt b/modules/mol/base/pymod/CMakeLists.txt
index 29c336dca643d9bbafac80ad9fb5e3959e6ef48c..87369f260e9e199fbcd107304f729f616d0d45c6 100644
--- a/modules/mol/base/pymod/CMakeLists.txt
+++ b/modules/mol/base/pymod/CMakeLists.txt
@@ -5,6 +5,7 @@ export_bond.cc
export_chain.cc
export_chain_view.cc
export_coord_group.cc
+export_coord_frame.cc
export_editors.cc
export_entity.cc
export_entity_view.cc
diff --git a/modules/mol/base/pymod/__init__.py b/modules/mol/base/pymod/__init__.py
index 75c503acf6f3ec1d787879f494bece9bf00e4359..8eb77fa6d3055f6259d7d1aa4c31ccffb9da7098 100644
--- a/modules/mol/base/pymod/__init__.py
+++ b/modules/mol/base/pymod/__init__.py
@@ -32,4 +32,4 @@ def MergeCoordGroups(*coord_groups):
cg=CreateCoordGroup(coord_groups[0].atoms)
for coord_group in coord_groups:
cg.AddFrames(coord_group)
- return cg
\ No newline at end of file
+ return cg
diff --git a/modules/mol/base/pymod/export_atom.cc b/modules/mol/base/pymod/export_atom.cc
index 42e18e0ce5f720e439429536192f10576608c5bd..f232d498cc3f52d61419d0c441c59e530470f737 100644
--- a/modules/mol/base/pymod/export_atom.cc
+++ b/modules/mol/base/pymod/export_atom.cc
@@ -53,6 +53,9 @@ void export_Atom()
.add_property("pos",
make_function(&AtomBase::GetPos,
return_value_policy<copy_const_reference>()))
+ .add_property("original_pos",
+ make_function(&AtomBase::GetOriginalPos,
+ return_value_policy<copy_const_reference>()))
.add_property("name",
make_function(&AtomBase::GetName,
return_value_policy<copy_const_reference>()),
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
new file mode 100644
index 0000000000000000000000000000000000000000..16665e0bb1ea3759f7cbcc62da1e2a4931df6aa4
--- /dev/null
+++ b/modules/mol/base/pymod/export_coord_frame.cc
@@ -0,0 +1,50 @@
+//------------------------------------------------------------------------------
+// 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/coord_frame.hh>
+#include <ost/geom/vec3.hh>
+#include <ost/mol/entity_handle.hh>
+
+using namespace ost;
+using namespace ost::mol;
+
+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_rmsd)(const mol::EntityView& Reference_View, const mol::EntityView& sele_View) = &CoordFrame::GetRMSD;
+
+
+void export_CoordFrame()
+{
+ class_<CoordFrame>("CoordFrame",no_init)
+ .def("GetAtomPos", get_atom_pos)
+ .def("GetDistanceBetwAtoms", get_dist_atom)
+ .def("GetAngle", get_angle)
+ .def("GetDihedralAngle", get_dihedral)
+ .def("GetCenterOfMassPos", get_cm)
+ .def("GetDistanceBetwCenterOfMass", get_dist_cm)
+ .def("GetRMSD",get_rmsd)
+ ;
+ def("CreateCoordFrame",CreateCoordFrame);
+}
diff --git a/modules/mol/base/pymod/export_coord_group.cc b/modules/mol/base/pymod/export_coord_group.cc
index ebd0f1f33daa0fc703d92f4e5351626741b52f61..794db5363c6ae050262e51c29608977ef83d569a 100644
--- a/modules/mol/base/pymod/export_coord_group.cc
+++ b/modules/mol/base/pymod/export_coord_group.cc
@@ -53,6 +53,7 @@ void export_CoordGroup()
.def("GetAtomPos",&CoordGroupHandle::GetAtomPos)
.def("CopyFrame",&CoordGroupHandle::CopyFrame)
.add_property("atoms", &CoordGroupHandle::GetAtomList)
+ .add_property("entity", &CoordGroupHandle::GetEntity)
.def("IsValid", &CoordGroupHandle::IsValid)
.def("Capture", capture1)
.def("Capture", capture2)
diff --git a/modules/mol/base/pymod/export_editors.cc b/modules/mol/base/pymod/export_editors.cc
index 0b469fec3909ddcf45ccf9bd4f4e4e3b5d0b4c18..34d3a7fe8960cba5019702873fca197995f22eb0 100644
--- a/modules/mol/base/pymod/export_editors.cc
+++ b/modules/mol/base/pymod/export_editors.cc
@@ -26,6 +26,10 @@ using namespace boost::python;
using namespace ost;
using namespace ost::mol;
+#if OST_NUMPY_SUPPORT_ENABLED
+#include <numpy/arrayobject.h>
+#endif
+
namespace {
BondHandle (EditorBase::*connect_a)(const AtomHandle&,
@@ -53,11 +57,143 @@ void (ICSEditor::*rotate_torsion_a)(TorsionHandle, Real)=&ICSEditor::RotateTorsi
void (ICSEditor::*rotate_torsion_b)(const AtomHandle&, const AtomHandle&,
const AtomHandle&, const AtomHandle&,
Real)=&ICSEditor::RotateTorsionAngle;
-
+
+#if OST_NUMPY_SUPPORT_ENABLED
+template<typename T, bool O>
+void set_pos2_nc_t(XCSEditor& e, const AtomHandleList& alist, PyArrayObject* na)
+{
+ size_t count=0;
+ for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait,++count) {
+ if(O) {
+ e.SetAtomOriginalPos(*ait,geom::Vec3(static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,0))),
+ static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,1))),
+ static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,2)))));
+ } else {
+ e.SetAtomTransformedPos(*ait,geom::Vec3(static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,0))),
+ static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,1))),
+ static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,2)))));
+ }
+ }
+}
+
+template<bool O>
+void set_pos2_t(XCSEditor& e, const AtomHandleList& alist, object pyobj)
+{
+ size_t acount = alist.size();
+
+ if(!PyArray_Check(pyobj.ptr())) {
+ throw std::runtime_error("expected a numpy array");
+ return;
+ }
+ PyArrayObject* na=reinterpret_cast<PyArrayObject*>(pyobj.ptr());
+
+ if(PyArray_NDIM(na)!=2 || PyArray_DIM(na,0)!=int(acount) || PyArray_DIM(na,1)!=3) {
+ throw std::runtime_error("excpected a numpy array of shape (NAtoms, 3)");
+ return;
+ }
+
+ if(PyArray_ISCONTIGUOUS(na)) {
+ if(PyArray_TYPE(na)==NPY_FLOAT) {
+ if(O) {
+ e.SetAtomOriginalPos(alist,reinterpret_cast<float*>(PyArray_DATA(na)));
+ } else {
+ e.SetAtomTransformedPos(alist,reinterpret_cast<float*>(PyArray_DATA(na)));
+ }
+ } else if(PyArray_TYPE(na)==NPY_DOUBLE) {
+ if(O) {
+ e.SetAtomOriginalPos(alist,reinterpret_cast<double*>(PyArray_DATA(na)));
+ } else {
+ e.SetAtomTransformedPos(alist,reinterpret_cast<double*>(PyArray_DATA(na)));
+ }
+ } else {
+ throw std::runtime_error("expected a numpy array of type float or double");
+ return;
+ }
+ } else {
+ // non-contiguous
+#if 0
+ throw std::runtime_error("expected contiguous numpy array");
+#else
+ if(PyArray_TYPE(na)==NPY_FLOAT) {
+ set_pos2_nc_t<float,O>(e,alist,na);
+ } else if(PyArray_TYPE(na)==NPY_DOUBLE) {
+ set_pos2_nc_t<double,O>(e,alist,na);
+ } else {
+ throw std::runtime_error("expected a numpy array of type float or double");
+ return;
+ }
+#endif
+ }
+}
+#endif
+
+void set_pos(XCSEditor& e, object o1, object o2, bool trans)
+{
+ extract<AtomHandle> eah(o1);
+ extract<geom::Vec3> ev3(o2);
+ if(eah.check() && ev3.check()) {
+ if(trans) {
+ e.SetAtomTransformedPos(eah(),ev3());
+ } else {
+ e.SetAtomOriginalPos(eah(),ev3());
+ }
+ return;
+ }
+
+#if OST_NUMPY_SUPPORT_ENABLED
+
+ extract<AtomHandleList> eal(o1);
+ if(eal.check()) {
+ if(trans) {
+ set_pos2_t<false>(e,eal(),o2);
+ } else {
+ set_pos2_t<true>(e,eal(),o2);
+ }
+ return;
+ }
+
+ std::map<unsigned long,AtomHandle> amap;
+ EntityHandle eh=e.GetEntity();
+ for(AtomHandleIter ait=eh.AtomsBegin(), aite=eh.AtomsEnd(); ait!=aite; ++ait) {
+ amap[(*ait).GetIndex()]=*ait;
+ }
+
+ AtomHandleList alist;
+ for(int i=0;i<len(o1);++i) {
+ int gid = extract<int>(o1[i]);
+ std::map<unsigned long,AtomHandle>::iterator ait=amap.find(static_cast<unsigned long>(gid));
+ alist.push_back(ait==amap.end() ? AtomHandle() : ait->second);
+ }
+
+ if(trans) {
+ set_pos2_t<false>(e,alist,o2);
+ } else {
+ set_pos2_t<true>(e,alist,o2);
+ }
+
+#else
+ throw std::runtime_error("SetAtom*Pos(...,ndarray) not available, because numpy support not compiled in");
+#endif
+}
+
+void set_o_pos(XCSEditor& e, object o1, object o2)
+{
+ set_pos(e,o1,o2,false);
+}
+
+void set_t_pos(XCSEditor& e, object o1, object o2)
+{
+ set_pos(e,o1,o2,true);
+}
+
}
void export_Editors()
-{
+{
+#if OST_NUMPY_SUPPORT_ENABLED
+ import_array();
+#endif
+
class_<EditorBase>("EditorBase", no_init)
.def("InsertChain", &EditorBase::InsertChain)
.def("InsertAtom", &EditorBase::InsertAtom,
@@ -84,8 +220,9 @@ void export_Editors()
;
class_<XCSEditor, bases<EditorBase> >("XCSEditor", no_init)
- .def("SetAtomPos", &XCSEditor::SetAtomPos)
- .def("SetAtomOriginalPos", &XCSEditor::SetAtomOriginalPos)
+ .def("SetAtomPos", set_t_pos)
+ .def("SetAtomTransformedPos", set_t_pos)
+ .def("SetAtomOriginalPos", set_o_pos)
.def("ApplyTransform", &XCSEditor::ApplyTransform)
.def("SetTransform", &XCSEditor::SetTransform)
.def("UpdateICS", &XCSEditor::UpdateICS)
diff --git a/modules/mol/base/pymod/export_entity.cc b/modules/mol/base/pymod/export_entity.cc
index 163918c75642b0e1da67336dcad2619e30c84da8..6ae3727031bf27cd0ced5adcb638d20f0d7b83e2 100644
--- a/modules/mol/base/pymod/export_entity.cc
+++ b/modules/mol/base/pymod/export_entity.cc
@@ -30,6 +30,10 @@ using namespace ost::mol;
#include <ost/export_helper/generic_property_def.hh>
+#if OST_NUMPY_SUPPORT_ENABLED
+#include <numpy/arrayobject.h>
+#endif
+
namespace {
EntityHandle create1() { return CreateEntity(); }
@@ -64,11 +68,55 @@ ICSEditor depr_request_ics_editor(EntityHandle e, EditMode m)
return e.EditICS(m);
}
+bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
+{
+ return a1.GetIndex()<a2.GetIndex();
+}
+
+
+PyObject* get_pos2(EntityHandle& entity, bool id_sorted)
+{
+#if OST_NUMPY_SUPPORT_ENABLED
+ npy_intp dims[]={entity.GetAtomCount(),3};
+ PyObject* na = PyArray_SimpleNew(2,dims,NPY_FLOAT);
+ npy_float* nad = reinterpret_cast<npy_float*>(PyArray_DATA(na));
+ if(id_sorted) {
+ AtomHandleList alist = entity.GetAtomList();
+ std::sort(alist.begin(),alist.end(),less_index);
+ for(AtomHandleList::const_iterator it=alist.begin();it!=alist.end();++it,nad+=3) {
+ geom::Vec3 pos=(*it).GetPos();
+ nad[0]=static_cast<npy_float>(pos[0]);
+ nad[1]=static_cast<npy_float>(pos[1]);
+ nad[2]=static_cast<npy_float>(pos[2]);
+ }
+ } else {
+ for(AtomHandleIter it=entity.AtomsBegin();it!=entity.AtomsEnd();++it,nad+=3) {
+ geom::Vec3 pos=(*it).GetPos();
+ nad[0]=static_cast<npy_float>(pos[0]);
+ nad[1]=static_cast<npy_float>(pos[1]);
+ nad[2]=static_cast<npy_float>(pos[2]);
+ }
+ }
+ return na;
+#else
+ throw std::runtime_error("GetPositions disabled, since numpy support is not compiled in");
+ return 0;
+#endif
+}
+PyObject* get_pos1(EntityHandle& entity)
+{
+ return get_pos2(entity,true);
}
+} // ns
+
void export_Entity()
{
+#if OST_NUMPY_SUPPORT_ENABLED
+ import_array();
+#endif
+
class_<EntityBase> ent_base("EntityBase", no_init);
ent_base
.def(self_ns::str(self))
@@ -139,6 +187,11 @@ void export_Entity()
.def("IsTransformationIdentity",&EntityHandle::IsTransformationIdentity)
.def(self==self)
.def(self!=self)
+#if OST_NUMPY_SUPPORT_ENABLED
+ .def("GetPositions",get_pos1)
+ .def("GetPositions",get_pos2)
+ .add_property("positions",get_pos1)
+#endif
;
def("CreateEntity",create1);
diff --git a/modules/mol/base/pymod/wrap_mol.cc b/modules/mol/base/pymod/wrap_mol.cc
index 0ab6b45681dfc9fd1d6f3c272407c56dff3e4cee..f29d1377f4bc0ccf14289b50df67e994208e99e5 100644
--- a/modules/mol/base/pymod/wrap_mol.cc
+++ b/modules/mol/base/pymod/wrap_mol.cc
@@ -39,6 +39,7 @@ void export_AtomView();
void export_ResidueView();
void export_Editors();
void export_CoordGroup();
+void export_CoordFrame();
void export_PropertyID();
void export_BoundingBox();
void export_QueryViewWrapper();
@@ -65,6 +66,7 @@ BOOST_PYTHON_MODULE(_ost_mol)
export_EntityView();
export_Editors();
export_CoordGroup();
+ export_CoordFrame();
export_PropertyID();
export_BoundingBox();
export_QueryViewWrapper();
diff --git a/modules/mol/base/src/atom_base.cc b/modules/mol/base/src/atom_base.cc
index f25d40135a9770c0d9fc1fbad7c5641ed98c5e36..642ea5aa93baad015633ab13caedbb2f7e3f704d 100644
--- a/modules/mol/base/src/atom_base.cc
+++ b/modules/mol/base/src/atom_base.cc
@@ -42,25 +42,25 @@ const GenericPropContainerImpl* AtomBase::GpImpl() const
const String& AtomBase::GetName() const
{
this->CheckValidity();
- return impl_->GetName();
+ return impl_->Name();
}
void AtomBase::SetName(const String& atom_name)
{
this->CheckValidity();
- return impl_->SetName(atom_name);
+ impl_->Name()=atom_name;
}
const geom::Vec3& AtomBase::GetPos() const
{
this->CheckValidity();
- return impl_->GetPos();
+ return impl_->TransformedPos();
}
const geom::Vec3& AtomBase::GetOriginalPos() const
{
this->CheckValidity();
- return impl_->GetOriginalPos();
+ return impl_->OriginalPos();
}
geom::Vec3 AtomBase::GetAltPos(const String& alt_group) const
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..80105c128fb3805765c0ccdd7c2f6f5095242ea6 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -0,0 +1,183 @@
+//------------------------------------------------------------------------------
+// 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/invalid_handle.hh>
+#include <ost/integrity_error.hh>
+#include <ost/log.hh>
+#include <ost/mol/in_mem_coord_source.hh>
+#include <ost/mol/view_op.hh>
+#include <ost/mol/mol.hh>
+#include "coord_frame.hh"
+
+namespace ost { namespace mol {
+
+ CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos)
+ {
+ return CoordFrame(atom_pos);
+ }
+
+ geom::Vec3 CoordFrame::GetAtomPos(const AtomHandle& atom){
+ //Returns the position of the atom in this frame
+ return this->GetAtomPos(atom.GetIndex());
+ }
+ geom::Vec3 CoordFrame::GetAtomPos(int i1){
+ //Returns the position in this frame of the atom with index i1
+ return (*this)[i1];
+ }
+
+ Real CoordFrame::GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2){
+ //Return the distance in this frame between two atoms
+ 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
+ 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)
+ 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
+ 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
+ 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 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
+ 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 v;
+ for (unsigned int i=0,e=indices.size();i!=e; i++) {
+ v+=masses[i]*(*this)[indices[i]];
+ }
+ return v;
+ }
+
+ Real CoordFrame::GetDistanceBetwCenterOfMass(const EntityView& sele1,const EntityView& sele2){
+ //Returns the distance between the centers of mass of the two EntityViews
+ std::vector<unsigned long> indices1,indices2;
+ std::vector<Real> masses1,masses2;
+ GetIndicesAndMasses(sele1,indices1,masses1);
+ GetIndicesAndMasses(sele2,indices2,masses2);
+ 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
+ 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 rmsd=0.0,val;
+ for (unsigned int i1=0; i1!=indices_sele.size(); ++i1) {
+ geom::Vec3 av_sele=(*this)[indices_sele[i1]];
+ geom::Vec3 av_ref=ref_pos[i1];
+ val=geom::Length2(av_ref-av_sele);
+ rmsd+=val;
+ }
+ 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
+ int count_ref=reference_view.GetAtomCount();
+ int count_sele=sele_view.GetAtomCount();
+ if (count_ref!=count_sele){
+ throw std::runtime_error("atom counts of the two views are not equal");
+ }
+ std::vector<unsigned long> indices_sele;
+ std::vector<geom::Vec3> ref_pos;
+ GetIndices(sele_view,indices_sele);
+ GetPositions(reference_view,ref_pos);
+ return this->GetRMSD(ref_pos,indices_sele);
+ }
+
+}}//ns
\ No newline at end of file
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index c55b855a45ee8508aacece0145a04678dcd68ef0..c856e583c180bf9dbb73cbbf6ac179f26dbdd466 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -20,19 +20,55 @@
#define OST_MOL_COORD_FRAME_HH
/*
- Author: Marco Biasini
+ Author: Marco Biasini and Niklaus Johner
*/
#include <boost/shared_ptr.hpp>
#include <ost/geom/geom.hh>
#include <ost/mol/module_config.hh>
+#include <ost/mol/entity_view.hh>
+#include "atom_handle.hh"
namespace ost { namespace mol {
-
-typedef std::vector<geom::Vec3> CoordFrame;
+class DLLEXPORT_OST_MOL CoordFrame : public geom::Vec3List{
+public:
+ typedef geom::Vec3List 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) { }
+
+ 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);
+ 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);
+};
+
+ void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
+ void GetMasses(const EntityView& sele, std::vector<Real>& masses);
+ void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
+ void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
+
typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
typedef std::vector<CoordFramePtr> CoordFrameList;
+// factory method
+// create a frame froma Vec3List containing the positions of the atoms
+ DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos);
+
}}
#endif
diff --git a/modules/mol/base/src/coord_group.cc b/modules/mol/base/src/coord_group.cc
index a68985d6ceab558e821b6d6f05dfb8eae19cc00b..f92c58e22b44f9dc4cd81753e288acefd164be9a 100644
--- a/modules/mol/base/src/coord_group.cc
+++ b/modules/mol/base/src/coord_group.cc
@@ -24,6 +24,7 @@
#include <ost/mol/mol.hh>
#include "coord_group.hh"
+
namespace ost { namespace mol {
CoordGroupHandle CreateCoordGroup(const AtomHandleList& atomlist)
@@ -83,8 +84,9 @@ CoordGroupHandle::operator bool() const
return this->IsValid();
}
-void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
-{
+//void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
+ void CoordGroupHandle::AddFrame(const geom::Vec3List& clist)
+ {
this->CheckValidity();
if (source_->IsMutable()) {
source_->AddFrame(clist);
@@ -203,5 +205,6 @@ CoordGroupHandle CoordGroupHandle::Filter(const EntityView& selected) const
}
return filtered_cg;
}
-
+
+
}} // ns
diff --git a/modules/mol/base/src/coord_group.hh b/modules/mol/base/src/coord_group.hh
index a3e582b15c12bfcc787d21a30b743cc2c04c1488..2b6d6fe7f4c40a68e323151f299ffeb184943e45 100644
--- a/modules/mol/base/src/coord_group.hh
+++ b/modules/mol/base/src/coord_group.hh
@@ -25,7 +25,6 @@
#include <vector>
#include <boost/shared_array.hpp>
-
#include "atom_handle.hh"
#include "coord_source.hh"
@@ -65,7 +64,8 @@ public:
void Capture(uint frame);
/// \brief add frame
- void AddFrame(const std::vector<geom::Vec3>& clist);
+ //void AddFrame(const std::vector<geom::Vec3>& clist);
+ void AddFrame(const geom::Vec3List& clist);
void AddFrames(const CoordGroupHandle& cg);
/// \brief set an indidivial atom position in the given frame
@@ -88,9 +88,10 @@ public:
CoordGroupHandle Filter(const EntityView& selected) const;
CoordGroupHandle(CoordSourcePtr source);
+
+
private:
void CheckValidity() const;
-
CoordSourcePtr source_;
};
diff --git a/modules/mol/base/src/editor_base.cc b/modules/mol/base/src/editor_base.cc
index c5eed6d60790041bff56e16e0751fe63da9db0f6..7718234673d7e35c33aee25b113a33e1c2dc2d54 100644
--- a/modules/mol/base/src/editor_base.cc
+++ b/modules/mol/base/src/editor_base.cc
@@ -149,7 +149,7 @@ void EditorBase::ReorderAllResidues()
void EditorBase::RenameAtom(AtomHandle atom, const String& new_name)
{
CheckHandleValidity(atom);
- atom.Impl()->SetName(new_name);
+ atom.Impl()->Name()=new_name;
}
BondHandle EditorBase::Connect(const AtomHandle& first,
@@ -202,11 +202,6 @@ TorsionHandle EditorBase::AddTorsion(const String& name, const AtomHandle& a1,
}
-EditMode EditorBase::GetMode() const
-{
- return mode_;
-}
-
void EditorBase::UpdateTrace()
{
if (mode_==UNBUFFERED_EDIT) {
diff --git a/modules/mol/base/src/editor_base.hh b/modules/mol/base/src/editor_base.hh
index e9409ecd31b4a8d53e25455045f91192856ffbec..48bf067bd19ea8c1957d60fe3ba3c439f86b4917 100644
--- a/modules/mol/base/src/editor_base.hh
+++ b/modules/mol/base/src/editor_base.hh
@@ -192,8 +192,11 @@ public:
void ReorderAllResidues();
/// \brief Get edit mode of editor
- EditMode GetMode() const;
+ EditMode GetMode() const {return mode_;}
+ /// \ brief return entity this editor works on
+ EntityHandle GetEntity() const {return ent_;}
+
/// \brief change the name of the atom to the new name
void RenameAtom(AtomHandle atom, const String& new_name);
protected:
diff --git a/modules/mol/base/src/entity_handle.hh b/modules/mol/base/src/entity_handle.hh
index 559a2cf109fd1a1089d38ead2b60bc4781249584..fff9f6efc306b77df8ca9aa6f6b804f1480329c0 100644
--- a/modules/mol/base/src/entity_handle.hh
+++ b/modules/mol/base/src/entity_handle.hh
@@ -137,7 +137,7 @@ public:
///
/// The two iterators returned by ResiduesBegin() and ResiduesEnd() form a
/// half-closed range. It is cheaper to cache the iterator returned by
- /// ResiduesEnd() than to call ResiduesEnd() after very loop, i.e. like
+ /// ResiduesEnd() than to call ResiduesEnd() after every loop, i.e. like
/// \code
/// // e is an instance of EntityHandle
/// for (ResidueHandleIter i=e.ResiduesBegin(), x=e.ResiduesEnd(); i!=x; ++i) {
diff --git a/modules/mol/base/src/impl/atom_impl.cc b/modules/mol/base/src/impl/atom_impl.cc
index bf66239bb46741d3de82d4dbd2a1384717f59a46..4f7f059b818eda692de6447a9b3d2b5d547037e7 100644
--- a/modules/mol/base/src/impl/atom_impl.cc
+++ b/modules/mol/base/src/impl/atom_impl.cc
@@ -94,11 +94,11 @@ void AtomImpl::TraceDirectionality(FragmentImplP frag, ConnectorImplP conn,
#if !defined(NDEBUG)
if (conn->GetFirst()==shared_from_this()) {
LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber()
- << "." << GetName() << " [" << conn->GetSecond()->GetQualifiedName()
+ << "." << Name() << " [" << conn->GetSecond()->GetQualifiedName()
<< " ]");
} else {
LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber()
- << "." << GetName() << " [" << conn->GetFirst()->GetQualifiedName()
+ << "." << Name() << " [" << conn->GetFirst()->GetQualifiedName()
<< " ]");
}
@@ -106,7 +106,7 @@ void AtomImpl::TraceDirectionality(FragmentImplP frag, ConnectorImplP conn,
#endif
} else {
LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber()
- << "." << GetName() << " [ ]");
+ << "." << Name() << " [ ]");
}
// presence of a primary connector indicates ring closure
@@ -195,7 +195,7 @@ void AtomImpl::UpdateFromXCS()
// stack before calling UpdateFromICS() on the next atom.
{
// Derive direction and length of connector from atom positions.
- geom::Vec3 global_d=((*i)->GetSecond()->GetOriginalPos()-this->GetOriginalPos());
+ geom::Vec3 global_d=((*i)->GetSecond()->OriginalPos()-this->OriginalPos());
// Set direction and length of connector. Direction is relative to
// local coordinate system of this atom.
// Note the order of arguments for the matrix multiplication. This is the
@@ -224,7 +224,7 @@ std::ostream& operator<<(std::ostream& o, const AtomImplPtr ap)
{
o << ap->GetResidue()->GetChain()->GetName() << ".";
o << ap->GetResidue()->GetKey() << ap->GetResidue()->GetNumber() << ".";
- o << ap->GetName();
+ o << ap->Name();
return o;
}
@@ -257,7 +257,7 @@ ConnectorImplP GetConnector(const AtomImplPtr& a, const AtomImplPtr& b) {
}
String AtomImpl::GetQualifiedName() const {
- return this->GetResidue()->GetQualifiedName()+"."+this->GetName();
+ return this->GetResidue()->GetQualifiedName()+"."+this->Name();
}
void AtomImpl::DeleteAllConnectors() {
diff --git a/modules/mol/base/src/impl/atom_impl.hh b/modules/mol/base/src/impl/atom_impl.hh
index 963f5ab60a3b0a1c5a8ccf714f0d5860e275799c..b342c526ae56bc5acdfbef3942255c7ff0a9f69c 100644
--- a/modules/mol/base/src/impl/atom_impl.hh
+++ b/modules/mol/base/src/impl/atom_impl.hh
@@ -56,21 +56,21 @@ public:
~AtomImpl();
void Apply(EntityVisitor& h);
- const String& GetName() const {return name_;}
-
- void SetName(const String& atom_name) {
- name_=atom_name;
- }
+ // for efficiency reasons, the simple setter/getter methods are
+ // replaced by direct access - this is the impl layer after all
+ const String& Name() const {return name_;}
+ String& Name() {return name_;}
- const geom::Vec3& GetPos() const {return tf_pos_;}
+ // DEPRECATED
+ const String& GetName() const {return name_;}
- const geom::Vec3& GetOriginalPos() const {return pos_;}
+ const geom::Vec3& TransformedPos() const {return tf_pos_;}
+ geom::Vec3& TransformedPos() {return tf_pos_;}
- void SetTransformedPos(const geom::Vec3& pos) { tf_pos_=pos; }
+ const geom::Vec3& OriginalPos() const {return pos_;}
+ geom::Vec3& OriginalPos() {return pos_;}
- void SetOriginalPos(const geom::Vec3& pos) { pos_=pos; }
-
ResidueImplPtr GetResidue() const;
void SetPrimaryConnector(const ConnectorImplP& bp) {
diff --git a/modules/mol/base/src/impl/chain_impl.cc b/modules/mol/base/src/impl/chain_impl.cc
index e0f51ada14d03f8b6fee8457c71aa8536471caf7..421f98dd4a43bafb12537e1649b51aba49abb9d3 100644
--- a/modules/mol/base/src/impl/chain_impl.cc
+++ b/modules/mol/base/src/impl/chain_impl.cc
@@ -420,8 +420,8 @@ geom::AlignedCuboid ChainImpl::GetBounds() const
ResidueImplPtr r=*i;
for (AtomImplList::iterator j=r->GetAtomList().begin();
j!=r->GetAtomList().end(); ++j) {
- mmin=geom::Min(mmin, (*j)->GetPos());
- mmax=geom::Max(mmax, (*j)->GetPos());
+ mmin=geom::Min(mmin, (*j)->TransformedPos());
+ mmax=geom::Max(mmax, (*j)->TransformedPos());
atoms=true;
}
}
@@ -441,7 +441,7 @@ geom::Vec3 ChainImpl::GetCenterOfAtoms() const
ResidueImplPtr r=*i;
for (AtomImplList::iterator j=r->GetAtomList().begin();
j!=r->GetAtomList().end(); ++j) {
- sum+=(*j)->GetPos();
+ sum+=(*j)->TransformedPos();
}
}
sum/=this->GetAtomCount();
@@ -459,7 +459,7 @@ geom::Vec3 ChainImpl::GetCenterOfMass() const
ResidueImplPtr r=*i;
for (AtomImplList::iterator j=r->GetAtomList().begin();
j!=r->GetAtomList().end(); ++j) {
- center+=(*j)->GetPos() * (*j)->GetMass();
+ center+=(*j)->TransformedPos() * (*j)->GetMass();
}
}
center/=mass;
diff --git a/modules/mol/base/src/impl/connector_impl.cc b/modules/mol/base/src/impl/connector_impl.cc
index dd4b285da0903649afb63cda9c9c4bab815e804b..8bafe009b6a110edecd96cadc3492ed348dc1f19 100644
--- a/modules/mol/base/src/impl/connector_impl.cc
+++ b/modules/mol/base/src/impl/connector_impl.cc
@@ -77,14 +77,14 @@ geom::Mat3 find_rotation(const geom::Vec3& d) {
geom::Vec3 ConnectorImpl::GetPos() const
{
- return (GetFirst()->GetPos()+GetSecond()->GetPos())*0.5;
+ return (GetFirst()->TransformedPos()+GetSecond()->TransformedPos())*0.5;
}
Real ConnectorImpl::GetLength() const
{
Real length;
if (this->GetFirst()->GetEntity()->HasICS()==false) {
- length=geom::Distance(this->GetFirst()->GetOriginalPos(),this->GetSecond()->GetOriginalPos());
+ length=geom::Distance(this->GetFirst()->OriginalPos(),this->GetSecond()->OriginalPos());
} else {
length=len_;
}
@@ -110,8 +110,8 @@ bool ConnectorImpl::IsConnectorOf(const AtomImplPtr& a,
geom::Vec3 ConnectorImpl::GetOriginalPos() const
{
- return (this->GetFirst()->GetOriginalPos()+
- this->GetSecond()->GetOriginalPos())*0.5;
+ return (this->GetFirst()->OriginalPos()+
+ this->GetSecond()->OriginalPos())*0.5;
}
void ConnectorImpl::Switch()
diff --git a/modules/mol/base/src/impl/dihedral.cc b/modules/mol/base/src/impl/dihedral.cc
index 3580829b198937ec8523481f041ed48b97f38faf..cac771039b96ae38e478edb7c9a596e337bd869c 100644
--- a/modules/mol/base/src/impl/dihedral.cc
+++ b/modules/mol/base/src/impl/dihedral.cc
@@ -46,8 +46,8 @@ Dihedral::Dihedral(const AtomImplList& atoms)
Real Dihedral::GetAngleXCS() const {
- return geom::DihedralAngle(atoms_[0]->GetPos(), atoms_[1]->GetPos(),
- atoms_[2]->GetPos(), atoms_[3]->GetPos());
+ return geom::DihedralAngle(atoms_[0]->TransformedPos(), atoms_[1]->TransformedPos(),
+ atoms_[2]->TransformedPos(), atoms_[3]->TransformedPos());
}
Real Dihedral::GetAngleICS() const {
diff --git a/modules/mol/base/src/impl/entity_impl.cc b/modules/mol/base/src/impl/entity_impl.cc
index 26709b43130b4119502ccc41a3d7d2af8814d9c2..1441af7cc3a079e5d08589951aff5e8f7be14e76 100644
--- a/modules/mol/base/src/impl/entity_impl.cc
+++ b/modules/mol/base/src/impl/entity_impl.cc
@@ -270,10 +270,10 @@ geom::AlignedCuboid EntityImpl::GetBounds() const
if (this->GetAtomCount()>0) {
geom::Vec3 mmin, mmax;
AtomImplMap::const_iterator it=atom_map_.begin();
- mmin=mmax=it->second->GetPos();
+ mmin=mmax=it->second->TransformedPos();
for (++it; it!=atom_map_.end();++it) {
- mmin=geom::Min(mmin,it->second->GetPos());
- mmax=geom::Max(mmax,it->second->GetPos());
+ mmin=geom::Min(mmin,it->second->TransformedPos());
+ mmax=geom::Max(mmax,it->second->TransformedPos());
}
return geom::AlignedCuboid(mmin, mmax);
} else {
@@ -287,7 +287,7 @@ geom::Vec3 EntityImpl::GetCenterOfAtoms() const {
if (this->GetAtomCount()>0) {
for (AtomImplMap::const_iterator it = atom_map_.begin();
it!=atom_map_.end();++it) {
- center+=it->second->GetPos();
+ center+=it->second->TransformedPos();
}
center/=static_cast<Real>(atom_map_.size());
}
@@ -299,7 +299,7 @@ geom::Vec3 EntityImpl::GetCenterOfMass() const {
Real mass = this->GetMass();
if (this->GetAtomCount()>0 && mass>0) {
for(AtomImplMap::const_iterator it = atom_map_.begin();it!=atom_map_.end();++it) {
- center+=it->second->GetPos()*it->second->GetMass();
+ center+=it->second->TransformedPos()*it->second->GetMass();
}
center/=mass;
}
@@ -326,12 +326,12 @@ AtomImplPtr EntityImpl::CreateAtom(const ResidueImplPtr& rp,
#else
AtomImplPtr ap(new AtomImpl(shared_from_this(), rp, name, pos, ele, next_index_++));
#endif
- if (identity_transf_ == false) {
+ if (!identity_transf_) {
geom::Vec3 transformed_pos = geom::Vec3(transformation_matrix_*geom::Vec4(pos));
- ap->SetTransformedPos(transformed_pos);
+ ap->TransformedPos()=transformed_pos;
atom_organizer_.Add(ap,transformed_pos);
} else {
- ap->SetTransformedPos(pos);
+ ap->TransformedPos()=pos;
atom_organizer_.Add(ap,pos);
}
atom_map_.insert(AtomImplMap::value_type(ap.get(),ap));
@@ -498,8 +498,8 @@ Real EntityImpl::GetAngleXCS(const AtomImplPtr& a1, const AtomImplPtr& a2,
ConnectorImplP c23=GetConnector(a2, a3);
if (c12 && c12->GetFirst() && c12->GetSecond()) {
if (c23 && c23->GetFirst() && c23->GetSecond()) {
- return Angle(a2->GetPos()-a1->GetPos(),
- a2->GetPos()-a3->GetPos());
+ return Angle(a2->TransformedPos()-a1->TransformedPos(),
+ a2->TransformedPos()-a3->TransformedPos());
} else {
AtomHandle ah2(a2), ah3(a3);
throw NotConnectedError(ah2, ah3);
@@ -765,12 +765,7 @@ void EntityImpl::SetTransform(const geom::Mat4 transfmat)
{
transformation_matrix_=transfmat;
inverse_transformation_matrix_=Invert(transformation_matrix_);
- geom::Mat4 identity = geom::Mat4();
- if (transformation_matrix_ == identity) {
- identity_transf_ = true;
- } else {
- identity_transf_ = false;
- }
+ identity_transf_ = (transformation_matrix_==geom::Mat4());
this->UpdateTransformedPos();
this->MarkOrganizerDirty();
}
@@ -1049,7 +1044,7 @@ void EntityImpl::UpdateOrganizer()
atom_organizer_.Clear();
for (AtomImplMap::const_iterator i=atom_map_.begin(),
e=atom_map_.end(); i!=e; ++i) {
- atom_organizer_.Add(i->second, i->second->GetPos());
+ atom_organizer_.Add(i->second, i->second->TransformedPos());
}
dirty_flags_&=~DirtyOrganizer;
}
@@ -1189,7 +1184,7 @@ void EntityImpl::RenameChain(ChainImplPtr chain, const String& new_name)
void EntityImpl::UpdateTransformedPos(){
for(AtomImplMap::iterator it = atom_map_.begin();it!=atom_map_.end();++it) {
- it->second->SetTransformedPos(geom::Vec3(transformation_matrix_*geom::Vec4(it->second->GetOriginalPos())));
+ it->second->TransformedPos()=geom::Vec3(transformation_matrix_*geom::Vec4(it->second->OriginalPos()));
}
}
diff --git a/modules/mol/base/src/impl/query_ast.cc b/modules/mol/base/src/impl/query_ast.cc
index 697ff173f1271b36e576d37acff20e9e50df5d58..ea96a485877b9c01226f4acdfa94e69c526fcb99 100644
--- a/modules/mol/base/src/impl/query_ast.cc
+++ b/modules/mol/base/src/impl/query_ast.cc
@@ -120,6 +120,7 @@ StringOrRegexParam::StringOrRegexParam():
StringOrRegexParam::StringOrRegexParam(const String& s):
is_regex_(false),r_(),s_(s)
{
+ String special("[]{}()");
for(String::const_iterator it=s.begin();it!=s.end();++it) {
if((*it)=='?' || (*it)=='*') {
is_regex_=true;
@@ -130,12 +131,16 @@ StringOrRegexParam::StringOrRegexParam(const String& s):
if(is_regex_) {
std::ostringstream e;
for(String::const_iterator it=s.begin();it!=s.end();++it) {
- if((*it)=='?') {
- e << ".";
- } else if((*it)=='*') {
- e << ".*";
+
+ if((*it)=='?' && (it==s.begin() || (*(it-1))!='\\')) {
+ e << ".";
+ } else if((*it)=='*' && (it==s.begin() || (*(it-1))!='\\')) {
+ e << ".*";
} else {
- e << *it;
+ if (special.find(*it)!=String::npos) {
+ e << '\\';
+ }
+ e << *it;
}
}
//std::cerr << "assembling regex [" << e.str() << "]... ";
diff --git a/modules/mol/base/src/impl/query_impl.cc b/modules/mol/base/src/impl/query_impl.cc
index cc5f80240bd1ddaff070ea42d518891c84419772..a70cfa813b39d63d9337cfb6135846d63796f8fb 100644
--- a/modules/mol/base/src/impl/query_impl.cc
+++ b/modules/mol/base/src/impl/query_impl.cc
@@ -148,7 +148,7 @@ QueryToken QueryLexer::LexNumericToken() {
}
bool is_ident_or_str(char c) {
- static String allowed_chars("_*?");
+ static String allowed_chars("_*?\\");
return isalnum(c) || allowed_chars.find_first_of(c)!=String::npos;
}
@@ -158,7 +158,8 @@ QueryToken QueryLexer::LexIdentOrStringToken() {
bool force_string=false;
while (current_<query_string_.length() &&
is_ident_or_str(query_string_[current_])) {
- if (query_string_[current_]=='*' || query_string_[current_]=='?') {
+ if (query_string_[current_]=='*' || query_string_[current_]=='?' ||
+ query_string_[current_]=='\\') {
force_string=true;
}
current_++;
diff --git a/modules/mol/base/src/impl/residue_impl.cc b/modules/mol/base/src/impl/residue_impl.cc
index cf4729239627e27f7f21477d7bbe7e9775ae2689..d8ce7f6d3f97b6c0bfe9e452d0c2760d54bbd01c 100644
--- a/modules/mol/base/src/impl/residue_impl.cc
+++ b/modules/mol/base/src/impl/residue_impl.cc
@@ -60,8 +60,8 @@ AtomImplPtr ResidueImpl::InsertAtom(const String& name,
AtomImplPtr ResidueImpl::InsertAtom(const AtomImplPtr& atom)
{
- AtomImplPtr dst_atom=this->InsertAtom(atom->GetName(),
- atom->GetPos(),
+ AtomImplPtr dst_atom=this->InsertAtom(atom->Name(),
+ atom->TransformedPos(),
atom->GetElement());
dst_atom->Assign(*atom.get());
@@ -166,12 +166,12 @@ AtomImplPtr ResidueImpl::GetCentralAtom() const
if (chem_class_.IsNucleotideLinking()) {
for (AtomImplList::const_iterator it=atom_list_.begin();
it!=atom_list_.end();++it) {
- if((*it)->GetName()=="P") return *it;
+ if((*it)->Name()=="P") return *it;
}
} else if (chem_class_.IsPeptideLinking()) {
for (AtomImplList::const_iterator it=atom_list_.begin();
it!=atom_list_.end();++it) {
- if((*it)->GetName()=="CA") return *it;
+ if((*it)->Name()=="CA") return *it;
}
}
@@ -217,14 +217,14 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const
AtomImplPtr a1 = FindAtom("C");
AtomImplPtr a2 = FindAtom("O");
if(a1 && a2) {
- nrvo = geom::Normalize(a2->GetPos()-a1->GetPos());
+ nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos());
} else {
a1 = FindAtom("CB");
a2 = FindAtom("CA");
if(a1 && a2) {
- nrvo = geom::Normalize(a2->GetPos()-a1->GetPos());
+ nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos());
} else {
- geom::Vec3 v0=GetCentralAtom()->GetPos();
+ geom::Vec3 v0=GetCentralAtom()->TransformedPos();
nrvo=geom::Cross(geom::Normalize(v0),
geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1])));
LOG_VERBOSE("warning: could not find atoms for proper central normal calculation");
@@ -235,9 +235,9 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const
AtomImplPtr a2 = FindAtom("OP1");
AtomImplPtr a3 = FindAtom("OP2");
if(a1 && a2 && a3) {
- nrvo = geom::Normalize(a1->GetPos()-(a2->GetPos()+a3->GetPos())*.5);
+ nrvo = geom::Normalize(a1->TransformedPos()-(a2->TransformedPos()+a3->TransformedPos())*.5);
} else {
- geom::Vec3 v0=GetCentralAtom()->GetPos();
+ geom::Vec3 v0=GetCentralAtom()->TransformedPos();
nrvo=geom::Cross(geom::Normalize(v0),
geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1])));
LOG_VERBOSE("warning: could not find atoms for proper central normal calculation");
@@ -251,7 +251,7 @@ AtomImplPtr ResidueImpl::FindAtom(const String& aname) const
{
for (AtomImplList::const_iterator it=atom_list_.begin();
it!=atom_list_.end();++it) {
- if ((*it)->GetName()==aname) {
+ if ((*it)->Name()==aname) {
return *it;
}
}
@@ -366,18 +366,25 @@ void ResidueImpl::DeleteAtom(const AtomImplPtr& atom) {
}
}
+namespace {
+ struct aname_matcher {
+ aname_matcher(const String& n): aname(n) {}
+ bool operator()(AtomImplPtr& a) {return a->Name()==aname;}
+ String aname;
+ };
+}
+
void ResidueImpl::DeleteAtoms(const String& atom_name) {
AtomImplList::iterator i=atom_list_.begin();
for (; i!=atom_list_.end(); ++i) {
- if ((*i)->GetName()==atom_name) {
+ if ((*i)->Name()==atom_name) {
(*i)->DeleteAllTorsions();
(*i)->DeleteAllConnectors();
this->GetEntity()->DeleteAtom(*i);
}
}
AtomImplList::iterator new_end;
- new_end=std::remove_if(atom_list_.begin(), atom_list_.end(),
- bind(&AtomImpl::GetName, _1)==atom_name);
+ new_end=std::remove_if(atom_list_.begin(), atom_list_.end(), aname_matcher(atom_name));
atom_list_.erase(new_end, atom_list_.end());
}
@@ -452,10 +459,10 @@ geom::AlignedCuboid ResidueImpl::GetBounds() const
if (atom_list_.size()>0) {
AtomImplList::const_iterator i=atom_list_.begin();
- mmin=mmax=(*i)->GetPos();
+ mmin=mmax=(*i)->TransformedPos();
for (++i; i!=atom_list_.end(); ++i) {
- mmax=geom::Max(mmax,(*i)->GetPos());
- mmin=geom::Min(mmin,(*i)->GetPos());
+ mmax=geom::Max(mmax,(*i)->TransformedPos());
+ mmin=geom::Min(mmin,(*i)->TransformedPos());
}
return geom::AlignedCuboid(mmin, mmax);
} else {
@@ -469,7 +476,7 @@ geom::Vec3 ResidueImpl::GetCenterOfAtoms() const
if (!atom_list_.empty()) {
for (AtomImplList::const_iterator i=atom_list_.begin();
i!=atom_list_.end(); ++i) {
- sum+=(*i)->GetPos();
+ sum+=(*i)->TransformedPos();
}
sum/=atom_list_.size();
}
@@ -483,7 +490,7 @@ geom::Vec3 ResidueImpl::GetCenterOfMass() const
if (this->GetAtomCount() > 0 && mass > 0) {
for (AtomImplList::const_iterator i=atom_list_.begin();
i!=atom_list_.end(); ++i) {
- center+=(*i)->GetPos()*(*i)->GetMass();
+ center+=(*i)->TransformedPos()*(*i)->GetMass();
}
}
return center/mass;
@@ -566,7 +573,7 @@ bool ResidueImpl::SwitchAtomPos(const String& group) {
for (; k!=gr.atoms.end(); ++k) {
AtomGroupEntry& entry=*k;
assert(!entry.atom.expired());
- entry.pos=entry.atom.lock()->GetOriginalPos();
+ entry.pos=entry.atom.lock()->OriginalPos();
}
}
AtomGroup& agr=i->second;
@@ -575,11 +582,11 @@ bool ResidueImpl::SwitchAtomPos(const String& group) {
AtomGroupEntry& entry=*j;
assert(!entry.atom.expired());
- entry.atom.lock()->SetOriginalPos(entry.pos);
+ entry.atom.lock()->OriginalPos()=entry.pos;
EntityHandle ent = entry.atom.lock()->GetEntity();
geom::Mat4 transf_matrix = ent.GetTransformationMatrix();
geom::Vec3 transf_pos = geom::Vec3(transf_matrix*geom::Vec4(entry.pos));
- entry.atom.lock()->SetTransformedPos(transf_pos);
+ entry.atom.lock()->TransformedPos()=transf_pos;
}
curr_group_=group;
return true;
diff --git a/modules/mol/base/src/query_state.cc b/modules/mol/base/src/query_state.cc
index b90767a7e20de61fad90328283940bbc9e0da2f2..d8d4da028c2f727bde33cb8305a7bf8ae09a0b08 100644
--- a/modules/mol/base/src/query_state.cc
+++ b/modules/mol/base/src/query_state.cc
@@ -319,7 +319,7 @@ boost::logic::tribool QueryState::EvalAtom(const AtomImplPtr& a) {
int int_value;
switch (ss.sel_id) {
case Prop::ANAME:
- str_value = a->GetName();
+ str_value = a->Name();
s_[*i] = cmp_string(ss.comp_op,str_value,
boost::get<StringOrRegexParam>(ss.param));
break;
@@ -328,17 +328,17 @@ boost::logic::tribool QueryState::EvalAtom(const AtomImplPtr& a) {
s_[*i]=cmp_num<int>(ss.comp_op, int_value,boost::get<int>(ss.param));
break;
case Prop::AX:
- float_value=(a->GetPos())[0];
+ float_value=(a->TransformedPos())[0];
s_[*i]=cmp_num<Real>(ss.comp_op, float_value,
boost::get<float>(ss.param));
break;
case Prop::AY:
- float_value=(a->GetPos())[1];
+ float_value=(a->TransformedPos())[1];
s_[*i]=cmp_num<Real>(ss.comp_op, float_value,
boost::get<float>(ss.param));
break;
case Prop::AZ:
- float_value=(a->GetPos())[2];
+ float_value=(a->TransformedPos())[2];
s_[*i]=cmp_num<Real>(ss.comp_op, float_value,
boost::get<float>(ss.param));
break;
@@ -358,7 +358,7 @@ boost::logic::tribool QueryState::EvalAtom(const AtomImplPtr& a) {
boost::get<float>(ss.param));
break;
case Prop::WITHIN:
- s_[*i]= this->do_within(a->GetPos(),
+ s_[*i]= this->do_within(a->TransformedPos(),
boost::get<WithinParam>(ss.param),
ss.comp_op);
break;
diff --git a/modules/mol/base/src/xcs_editor.cc b/modules/mol/base/src/xcs_editor.cc
index f5a35fa66dda31b0230f5ee06db98b2b7688eb25..0c2454b187c8d29902b8380e4caf3b7969db3fe5 100644
--- a/modules/mol/base/src/xcs_editor.cc
+++ b/modules/mol/base/src/xcs_editor.cc
@@ -70,33 +70,121 @@ XCSEditor& XCSEditor::operator=(const XCSEditor& rhs)
return *this;
}
-void XCSEditor::SetAtomPos(const AtomHandle& atom,
- const geom::Vec3& position)
+void XCSEditor::SetAtomTransformedPos(const AtomHandle& atom,
+ const geom::Vec3& position)
{
CheckHandleValidity(atom);
- atom.Impl()->SetTransformedPos(position);
- geom::Mat4 inv_transformation_matrix = ent_.Impl()->GetInvTransfMatrix();
- geom::Vec3 original_pos = geom::Vec3(inv_transformation_matrix*geom::Vec4(position));
- atom.Impl()->SetOriginalPos(original_pos);
+ atom.Impl()->TransformedPos()=position;
+ if(ent_.Impl()->IsTransfIdentity()) {
+ atom.Impl()->OriginalPos()=position;
+ } else {
+ atom.Impl()->OriginalPos() = geom::Vec3(ent_.Impl()->GetInvTransfMatrix()*geom::Vec4(position));
+ }
ent_.Impl()->MarkICSDirty();
ent_.Impl()->MarkOrganizerDirty();
this->Update();
}
+namespace {
+ template<typename T>
+ void set_transformed_pos(impl::EntityImpl* ent, const AtomHandleList& alist, T *positions)
+ {
+ bool has_tf=ent->IsTransfIdentity();
+ for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait) {
+ if(ait->IsValid()) {
+ ait->Impl()->TransformedPos()[0]=static_cast<Real>(positions[0]);
+ ait->Impl()->TransformedPos()[1]=static_cast<Real>(positions[1]);
+ ait->Impl()->TransformedPos()[2]=static_cast<Real>(positions[2]);
+ if(has_tf) {
+ ait->Impl()->OriginalPos()=ait->Impl()->TransformedPos();
+ } else {
+ ait->Impl()->OriginalPos() = geom::Vec3(ent->GetInvTransfMatrix()*geom::Vec4(ait->Impl()->TransformedPos()));
+ }
+ }
+ positions+=3;
+ }
+ ent->MarkICSDirty();
+ ent->MarkOrganizerDirty();
+ }
+} // anon ns
+
+void XCSEditor::SetAtomTransformedPos(const AtomHandleList& alist, float *positions)
+{
+ set_transformed_pos<float>(ent_.Impl().get(),alist,positions);
+ this->Update();
+}
+
+void XCSEditor::SetAtomTransformedPos(const AtomHandleList& alist, double *positions)
+{
+ set_transformed_pos<double>(ent_.Impl().get(),alist,positions);
+ this->Update();
+}
void XCSEditor::SetAtomOriginalPos(const AtomHandle& atom,
- const geom::Vec3& position)
+ const geom::Vec3& position)
{
CheckHandleValidity(atom);
- atom.Impl()->SetOriginalPos(position);
- geom::Mat4 transformation_matrix = atom.GetEntity().GetTransformationMatrix();
- geom::Vec3 transformed_pos = geom::Vec3(transformation_matrix*geom::Vec4(position));
- atom.Impl()->SetTransformedPos(transformed_pos);
+ atom.Impl()->OriginalPos()=position;
+ if(ent_.Impl()->IsTransfIdentity()) {
+ atom.Impl()->TransformedPos()=position;
+ } else {
+ atom.Impl()->TransformedPos() = geom::Vec3(ent_.Impl()->GetTransfMatrix()*geom::Vec4(position));
+ }
ent_.Impl()->MarkICSDirty();
ent_.Impl()->MarkOrganizerDirty();
this->Update();
}
+namespace {
+ template<typename T>
+ void set_original_pos(impl::EntityImpl* ent, const AtomHandleList& alist, T *positions)
+ {
+ bool has_tf=ent->IsTransfIdentity();
+ for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait) {
+ if(ait->IsValid()) {
+ ait->Impl()->OriginalPos()[0]=static_cast<Real>(positions[0]);
+ ait->Impl()->OriginalPos()[1]=static_cast<Real>(positions[1]);
+ ait->Impl()->OriginalPos()[2]=static_cast<Real>(positions[2]);
+ if(has_tf) {
+ ait->Impl()->TransformedPos()=ait->Impl()->OriginalPos();
+ } else {
+ ait->Impl()->TransformedPos() = geom::Vec3(ent->GetTransfMatrix()*geom::Vec4(ait->Impl()->OriginalPos()));
+ }
+ }
+ positions+=3;
+ }
+ ent->MarkICSDirty();
+ ent->MarkOrganizerDirty();
+ }
+} // anon ns
+
+void XCSEditor::SetAtomOriginalPos(const AtomHandleList& alist, float *positions)
+{
+ set_original_pos<float>(ent_.Impl().get(),alist,positions);
+ this->Update();
+}
+
+void XCSEditor::SetAtomOriginalPos(const AtomHandleList& alist, double *positions)
+{
+ set_original_pos<double>(ent_.Impl().get(),alist,positions);
+ this->Update();
+}
+
+void XCSEditor::SetAtomPos(const AtomHandle& atom, const geom::Vec3& position)
+{
+ this->SetAtomTransformedPos(atom,position);
+}
+
+void XCSEditor::SetAtomPos(const AtomHandleList& alist, float *positions)
+{
+ this->SetAtomTransformedPos(alist,positions);
+}
+
+void XCSEditor::SetAtomPos(const AtomHandleList& alist, double *positions)
+{
+ this->SetAtomTransformedPos(alist,positions);
+}
+
void XCSEditor::ApplyTransform(const geom::Mat4& transform)
{
ent_.Impl()->ApplyTransform(transform);
diff --git a/modules/mol/base/src/xcs_editor.hh b/modules/mol/base/src/xcs_editor.hh
index 4b89083c6f82b5952a00270683a82bfddeada13c..73d96a652829932c469a6a0382093de7fd53b5f4 100644
--- a/modules/mol/base/src/xcs_editor.hh
+++ b/modules/mol/base/src/xcs_editor.hh
@@ -51,11 +51,47 @@ public:
void SetAtomOriginalPos(const AtomHandle& atom,
const geom::Vec3& position);
+ /// \brief numpy float interface
+ ///
+ /// the passed in float array must have a length of 3*alist.size()
+ void SetAtomOriginalPos(const AtomHandleList& alist,
+ float *positions);
+
+ /// \brief numpy double interface
+ ///
+ /// the passed in double array must have a length of 3*alist.size()
+ void SetAtomOriginalPos(const AtomHandleList& alist,
+ double *positions);
+
/// \brief set transformed position of atom
///
/// This function also updates the original coordinates
+ void SetAtomTransformedPos(const AtomHandle& atom,
+ const geom::Vec3& position);
+
+ /// \brief numpy float interface
+ ///
+ /// the passed in float array must have a length of 3*alist.size()
+ void SetAtomTransformedPos(const AtomHandleList& alist,
+ float *positions);
+
+ /// \brief numpy double interface
+ ///
+ /// the passed in double array must have a length of 3*alist.size()
+ void SetAtomTransformedPos(const AtomHandleList& alist,
+ double *positions);
+
+ /// \brief same as SetAtomTransformedPos(AtomHandle, geom::Vec3)
void SetAtomPos(const AtomHandle& atom,
- const geom::Vec3& position);
+ const geom::Vec3& position);
+
+ /// \brief same as SetTransformedPos(AtomHandleList,float*)
+ void SetAtomPos(const AtomHandleList& alist,
+ float *positions);
+
+ /// \brief same as SetTransformedPos(AtomHandleList,double*)
+ void SetAtomPos(const AtomHandleList& alist,
+ double *positions);
/// \brief apply additional transformation to all atoms
///
@@ -69,6 +105,7 @@ public:
/// \brief immediately update internal coordinate system
void UpdateICS();
+
protected:
XCSEditor(const EntityHandle& ent, EditMode mode);
diff --git a/modules/mol/base/tests/CMakeLists.txt b/modules/mol/base/tests/CMakeLists.txt
index e5d1c868ca25f6b2d138e176b37ea71b806458f8..b9771e0c30ec443cba99f7cc35ee8cf6ca773a9a 100644
--- a/modules/mol/base/tests/CMakeLists.txt
+++ b/modules/mol/base/tests/CMakeLists.txt
@@ -15,7 +15,7 @@ set(OST_MOL_BASE_UNIT_TESTS
tests.cc
)
-ost_unittest(mol "${OST_MOL_BASE_UNIT_TESTS}")
+ost_unittest(MODULE mol SOURCES "${OST_MOL_BASE_UNIT_TESTS}")
# for valgrind debugging
# executable(NAME test_query_standalone SOURCES test_query_standalone.cc DEPENDS_ON mol)
diff --git a/modules/mol/base/tests/test_numpy.py b/modules/mol/base/tests/test_numpy.py
new file mode 100644
index 0000000000000000000000000000000000000000..1255425f85ad0bfed89b238f7264d7ec394b367d
--- /dev/null
+++ b/modules/mol/base/tests/test_numpy.py
@@ -0,0 +1,61 @@
+import sys,unittest
+sys.path.append("../../../../stage/lib/openstructure")
+sys.path.append("../../../../stage/lib64/openstructure")
+from ost import *
+import numpy
+
+def v2v(v):
+ return geom.Vec3(float(v[0]),float(v[1]),float(v[2]))
+
+def dd(v1,v2):
+ return geom.Distance(v1,v2)<1e-8
+
+class TestNumpy(unittest.TestCase):
+ def setUp(self):
+ pass
+
+ def test_(self):
+ entity=mol.CreateEntity()
+ ed=entity.EditXCS()
+ ch=ed.InsertChain("X")
+ re=ed.AppendResidue(ch,"ALA")
+ a0=ed.InsertAtom(re,"A",geom.Vec3(0,0,0))
+ self.assertEqual(a0.GetIndex(),0)
+ a1=ed.InsertAtom(re,"B",geom.Vec3(1,0,0))
+ self.assertEqual(a1.GetIndex(),1)
+ a2=ed.InsertAtom(re,"C",geom.Vec3(2,0,0))
+ self.assertEqual(a2.GetIndex(),2)
+ a3=ed.InsertAtom(re,"D",geom.Vec3(3,0,0))
+ self.assertEqual(a3.GetIndex(),3)
+
+ self.assertTrue(dd(a0.pos,geom.Vec3(0,0,0)))
+ self.assertTrue(dd(a1.pos,geom.Vec3(1,0,0)))
+ self.assertTrue(dd(a2.pos,geom.Vec3(2,0,0)))
+ self.assertTrue(dd(a3.pos,geom.Vec3(3,0,0)))
+
+ ed.SetAtomTransformedPos(entity.GetAtomList(),
+ numpy.array([[0,1,0],[0,2,0],[0,3,0],[0,4,0]], dtype=numpy.float32))
+
+ self.assertTrue(dd(a0.pos,geom.Vec3(0,1,0)))
+ self.assertTrue(dd(a1.pos,geom.Vec3(0,2,0)))
+ self.assertTrue(dd(a2.pos,geom.Vec3(0,3,0)))
+ self.assertTrue(dd(a3.pos,geom.Vec3(0,4,0)))
+
+ na=entity.positions
+
+ self.assertTrue(dd(v2v(na[0]),geom.Vec3(0,1,0)))
+ self.assertTrue(dd(v2v(na[1]),geom.Vec3(0,2,0)))
+ self.assertTrue(dd(v2v(na[2]),geom.Vec3(0,3,0)))
+ self.assertTrue(dd(v2v(na[3]),geom.Vec3(0,4,0)))
+
+ ed.SetAtomTransformedPos([3,99,2],
+ numpy.array([[0,0,-3],[-1,-1,-1],[0,0,-2]], dtype=numpy.float32))
+
+ self.assertTrue(dd(a0.pos,geom.Vec3(0,1,0)))
+ self.assertTrue(dd(a1.pos,geom.Vec3(0,2,0)))
+ self.assertTrue(dd(a2.pos,geom.Vec3(0,0,-2)))
+ self.assertTrue(dd(a3.pos,geom.Vec3(0,0,-3)))
+
+if __name__== '__main__':
+ unittest.main()
+
diff --git a/modules/seq/alg/src/merge_pairwise_alignments.cc b/modules/seq/alg/src/merge_pairwise_alignments.cc
index 5636f32c97af7b96297663f75d6456abf954a4d4..d5da17127fc3e5233552a721d12acdaa1f1eb0f5 100644
--- a/modules/seq/alg/src/merge_pairwise_alignments.cc
+++ b/modules/seq/alg/src/merge_pairwise_alignments.cc
@@ -55,7 +55,8 @@ void update_shifts(const AlignmentHandle& aln,
res_index=-1;
}
else if (i > 0) {
- res_index=s1.GetResidueIndex(i-1);
+ res_index=s1.GetResidueIndex(i-1)-s1.GetOffset();
+// res_index=s1.GetResidueIndex(i-1);
}
ShiftMap::iterator p=shifts.find(res_index);
if (p!=shifts.end()) {
@@ -111,7 +112,8 @@ SequenceHandle realign_sequence(const AlignmentHandle& aln,
p=shifts.find(-1);
}
else if (i>0 && s1.GetOneLetterCode(i-1)!='-') {
- p=shifts.find(s1.GetResidueIndex(i-1));
+ p=shifts.find(s1.GetResidueIndex(i-1)-s1.GetOffset());
+// p=shifts.find(s1.GetResidueIndex(i-1));
}
if (p!=shifts.end()) {
int d=p->second-shift;
@@ -139,11 +141,16 @@ AlignmentHandle MergePairwiseAlignments(const AlignmentList& pairwise_alns,
update_shifts(*i, ref_seq, shifts);
}
+
AlignmentHandle merged=CreateAlignment();
merged.AddSequence(shift_reference(ref_seq, shifts));
+ size_t ref_len=merged.GetSequence(0).GetLength();
for (AlignmentList::const_iterator i=pairwise_alns.begin(),
e=pairwise_alns.end(); i!=e; ++i) {
SequenceHandle new_seq=realign_sequence(*i, shifts);
+ for (size_t j=new_seq.GetLength(); j<ref_len; ++j) {
+ new_seq.Append('-');
+ }
merged.AddSequence(new_seq);
}
return merged;
diff --git a/modules/seq/alg/tests/CMakeLists.txt b/modules/seq/alg/tests/CMakeLists.txt
index f6d8ebe1493d5198882f57f101956baa8e6ff830..9933524408a7fd8bcf4107c72217fe6632b7f3d3 100644
--- a/modules/seq/alg/tests/CMakeLists.txt
+++ b/modules/seq/alg/tests/CMakeLists.txt
@@ -7,5 +7,5 @@ set(OST_SEQ_ALG_UNIT_TESTS
test_global_align.py
)
-ost_unittest(seq_alg "${OST_SEQ_ALG_UNIT_TESTS}")
+ost_unittest(MODULE seq_alg SOURCES "${OST_SEQ_ALG_UNIT_TESTS}")
diff --git a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc b/modules/seq/alg/tests/test_merge_pairwise_alignments.cc
index 3f67ac5d161840199ac352ffb09d5bf7c883603d..f89715adeccc627198420c1fc80826502db617a6 100644
--- a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc
+++ b/modules/seq/alg/tests/test_merge_pairwise_alignments.cc
@@ -107,4 +107,30 @@ BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_three)
BOOST_CHECK_EQUAL(seqs[1].GetString(), "xyabcdefghijk");
BOOST_CHECK_EQUAL(seqs[2].GetString(), "-zabcdefghijk");
}
+
+BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_four)
+{
+ SequenceHandle ref=CreateSequence("REF", "abcdefghijk");
+ SequenceHandle s1=CreateSequence("S1", "abcdefghijk--");
+ SequenceHandle s2=CreateSequence("S2", "abcdefghijkxy");
+ AlignmentHandle aln1=CreateAlignment();
+ aln1.AddSequence(s1);
+ aln1.AddSequence(s2);
+
+ SequenceHandle s3=CreateSequence("S1", "-abcdefghijk");
+ SequenceHandle s4=CreateSequence("S2", "zabcdefghijk");
+
+ AlignmentHandle aln2=CreateAlignment();
+ aln2.AddSequence(s3);
+ aln2.AddSequence(s4);
+ AlignmentList l;
+ l.push_back(aln1);
+ l.push_back(aln2);
+ AlignmentHandle m=alg::MergePairwiseAlignments(l, ref);
+ ConstSequenceList seqs=m.GetSequences();
+ BOOST_CHECK_EQUAL(seqs[0].GetString(), "-abcdefghijk--");
+ BOOST_CHECK_EQUAL(seqs[1].GetString(), "-abcdefghijkxy");
+ BOOST_CHECK_EQUAL(seqs[2].GetString(), "zabcdefghijk--");
+}
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/modules/seq/base/pymod/export_sequence.cc b/modules/seq/base/pymod/export_sequence.cc
index df18efd0db8a778a6779d2082a674732a214eb59..0cb3f4df3d62de4a2b6311b8effa132c714c300f 100644
--- a/modules/seq/base/pymod/export_sequence.cc
+++ b/modules/seq/base/pymod/export_sequence.cc
@@ -341,7 +341,7 @@ void export_sequence()
.def("GetLength", &AlignmentHandle::GetLength)
.def("__len__", &AlignmentHandle::GetLength)
.def("GetSequences", &AlignmentHandle::GetSequences)
- .def("GetCoverage", &AlignmentHandle::GetSequences)
+ .def("GetCoverage", &AlignmentHandle::GetCoverage)
.def("AttachView", attach_view_a)
.def("AttachView", attach_view_b)
.def("Cut", &AlignmentHandle::Cut)
diff --git a/modules/seq/base/tests/CMakeLists.txt b/modules/seq/base/tests/CMakeLists.txt
index 528cb936ddae91eabe14d1b06c6aa743de7e29ac..645dba066e39c51c1fc0b9211228e8490ed85a07 100644
--- a/modules/seq/base/tests/CMakeLists.txt
+++ b/modules/seq/base/tests/CMakeLists.txt
@@ -7,5 +7,5 @@ set(OST_SEQ_UNIT_TESTS
tests.cc
)
-ost_unittest(seq "${OST_SEQ_UNIT_TESTS}")
+ost_unittest(MODULE seq SOURCES "${OST_SEQ_UNIT_TESTS}")