Add atom container align and rmsd methods (Zube 2914).

moldesign/_notebooks/Example 6. Align Two Molecules.ipynb

@@ -0,0 +1,264 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<span style=\"float:right\"><a href=\"http://moldesign.bionano.autodesk.com/\" target=\"_blank\" title=\"About\">About</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href=\"https://forum.bionano.autodesk.com/c/Molecular-Design-Toolkit\" target=\"_blank\" title=\"Forum\">Forum</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href=\"https://github.com/autodesk/molecular-design-toolkit/issues\" target=\"_blank\" title=\"Issues\">Issues</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href=\"http://bionano.autodesk.com/MolecularDesignToolkit/explore.html\" target=\"_blank\" title=\"Tutorials\">Tutorials</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href=\"http://autodesk.github.io/molecular-design-toolkit/\" target=\"_blank\" title=\"Documentation\">Documentation</a></span>\n",
+    "</span>\n",
+    "![Molecular Design Toolkit](img/Top.png)\n",
+    "\n",
+    "\n",
+    "<center><h1>Example 6: Align Two Molecules </h1> </center>\n",
+    "---\n",
+    "\n",
+    "\n",
+    "This notebook reads in two structures and calculates the transformation matrix aligning the C-alpha atoms from a subset of residues one structure to the other.\n",
+    "\n",
+    " - _Author_: [Dave Parker](https://github.com/ktbolt), Autodesk Research<br>\n",
+    " - _Created on_: August 23, 2017\n",
+    " - _Tags_: align, fit\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 70,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import moldesign as mdt\n",
+    "from moldesign import units as u\n",
+    "from moldesign.molecules.atomcollections import AtomList\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Read in and display the source molecule."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 71,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "2395475165f342f5b5d94da1027b4e26"
+      }
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "src_mol = mdt.from_pdb(\"1yu8\")\n",
+    "src_mol.draw()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Extract the source molecule C-alpha atoms for residues 64-72."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "src_residues = src_mol.chains[\"X\"].residues\n",
+    "src_residue_ids = set(range(64,73))\n",
+    "src_atoms = [res.atoms[\"CA\"] for res in src_residues if res.pdbindex in src_residue_ids]\n",
+    "src_atoms_list = AtomList(src_atoms)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Read in and display the destination molecule."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 76,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "WARNING: Residue GLU502 (index 271, chain A) is missing expected atoms. Attempting to infer chain end\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "4948f6a1dc234e02b35d375f52ffa115"
+      }
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "dst_mol = mdt.from_pdb(\"3ac2\")\n",
+    "dst_mol.draw()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Extract the destination molecule C-alpha atoms for residues 446-454."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 77,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "dst_residues = dst_mol.chains[\"A\"].residues\n",
+    "dst_residue_ids = set(range(446,455))\n",
+    "dst_atoms = [res.atoms[\"CA\"] for res in residues if res.pdbindex in residue_ids]\n",
+    "dst_atoms_list = AtomList(dst_atoms)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Calculate the transformation aligning the source atoms to destination atoms."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 88,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "rmsd = 0.412441\n",
+      "transformation = \n",
+      "[[  1.00000000e+00   4.56737751e-16   2.21965480e-16  -7.10542736e-15]\n",
+      " [ -3.52269708e-16   1.00000000e+00  -6.83377612e-17   8.88178420e-15]\n",
+      " [ -4.55009138e-16   5.86301682e-16   1.00000000e+00   8.88178420e-15]\n",
+      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00000000e+00]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "rmsd, xform = src_atoms_list.align(dst_atoms_list)\n",
+    "print(\"rmsd = %f\" % rmsd)\n",
+    "print(\"transformation = \\n\" + np.array_str(xform))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Transform all of the source atoms."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 79,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "src_atoms = src_mol.atoms\n",
+    "coords = np.array([0.0, 0.0, 0.0, 1.0],dtype=float)\n",
+    "\n",
+    "for i in range(0,len(src_atoms)):\n",
+    "    coords[0:3] = src_atoms[i].position.value_in(u.angstrom)\n",
+    "    xcoords = xform.dot(coords)\n",
+    "    src_atoms[i].x = xcoords[0] * u.angstrom\n",
+    "    src_atoms[i].y = xcoords[1] * u.angstrom\n",
+    "    src_atoms[i].z = xcoords[2] * u.angstrom\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Merge the source and destination molecules into s single molecule and display it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 80,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "WARNING: atom indices modified due to name clashes\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "888fa1f3e38b498ea3c11fc37be561f7"
+      }
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "merged_mol = mdt.Molecule([src_mol, dst_mol])\n",
+    "merged_mol.draw()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

moldesign/_tests/test_alignments.py

@@ -6,10 +6,11 @@
 import moldesign as mdt
 from moldesign import geom
 from moldesign.mathutils import normalized
+from moldesign.molecules.atomcollections import AtomList
 from moldesign import units as u
 
 from .molecule_fixtures import *
-from .helpers import assert_almost_equal
+from .helpers import assert_almost_equal, get_data_path
 
 
 __PYTEST_MARK__ = 'internal'  # mark all tests in this module with this label (see ./conftest.py)
@@ -146,3 +147,34 @@ def test_pmi_reorientation_on_benzene(benzene):
         np.testing.assert_allclose(newdistmat.defunits_value(),
                                    original_distmat)
 
+# Test the alignmnet of the atoms from two molecules.
+def test_align_molecules():
+
+    # Get the CA atoms for chain X residues 64-72.
+    mol1 = mdt.read(get_data_path('1yu8.pdb'))
+    residues = mol1.chains["X"].residues
+    residue_ids = set(range(64,73))
+    src_atoms = [res.atoms["CA"] for res in residues if res.pdbindex in residue_ids]
+
+    # Get the CA atoms for chain A residues 446-454.
+    mol2 = mdt.read(get_data_path('3ac2.pdb'))
+    residues = mol2.chains["A"].residues
+    residue_ids = set(range(446,455))
+    dst_atoms = [res.atoms["CA"] for res in residues if res.pdbindex in residue_ids]
+
+    # Align the two atom lists.
+    src_atoms_list = AtomList(src_atoms)
+    dst_atoms_list = AtomList(dst_atoms)
+    rmsd_error, xform = src_atoms_list.align(dst_atoms_list)
+
+    # Compare results to expected values.
+    rtol = 1e-5
+    expected_rmsd_error = 0.412441
+    expected_xform = np.array([(4.42719613e-01,  8.84833988e-01, -1.45148745e-01,  1.13863353e+01),   \
+                               (1.30739446e-02, -1.68229930e-01, -9.85661079e-01,  2.42427092e+01),   \
+                               (-8.96564786e-01,  4.34473825e-01, -8.60469602e-02,  2.24877465e+01),  \
+                               (0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.00000000e+00)])
+    np.testing.assert_allclose(rmsd_error, expected_rmsd_error, rtol)
+    np.testing.assert_allclose(xform, expected_xform, rtol)
+
+

moldesign/mathutils/__init__.py

@@ -1,3 +1,4 @@
 from .vectormath import *
 from .eigen import *
-from .grids import *
+from .grids import *
+from .align import *