{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9229fb0449cb4bd185eb93727e851d2c", "version_major": 2, "version_minor": 0 }, "text/plain": [] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from simtk.openmm import app, LangevinIntegrator, XmlSerializer\n", "from simtk import unit\n", "from openforcefield.typing.engines.smirnoff import ForceField as OFFForceField\n", "from openforcefield.topology import Molecule as OFFMolecule\n", "from openforcefield.topology import Topology as OFFTopology\n", "from openforcefield.topology import TopologyAtom\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: Fun with VirtualSites. Define a twofold sulfur lone pair virtual site and run a short simulation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "vsite_offxml = '''\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "'''\n", "ff = OFFForceField('openff-1.2.0.offxml', vsite_offxml)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1RTV/o38CcQ7iDhKiCCeEGFIChaL+AoilNr0XZUKGLxbWsbqtPBdjkzceysYqdT39jpbxY/3rFjqF2r8V6oo0Iro1AvgFZroNwF8cJFQCBAIEBIIDnvH7umyE0CSU4uz2f1D3Jycs5Dga/77L3PPgyKogAhhNBEmdFdAEIIGTaMUYQQmhSMUYQQmhSMUYQQmhSMUYQQmhSMUYQQmhSMUYQQmhSMUYQQmhSMUYT0EZvNPnHihOrlmTNn5s2bR2M9aAwYowghNCkYowghNCkYowghNCkMXJoEIT3EZrNra2utrKzIS7lc7uXlVVlZSW9VaETYGkVITx08eLDoKR6PR3c5aFRMugswGDEVFQ+lUtXLY/PnB9ja0lgPMnpOTk7e3t7ka2dnZ3qLQWPAGFXDRzNmrHdyIl9bmGFDHiEEgDGqFiaDYYnpiXTr8uXLHR0dMpmM7kLQqHCIabxiKire8PDYgNdWSLd8fX3r6upqamp8fX3prgWNDFujauDV1v5PfT0A+FlbH507l+5ykEmwt7cHgO7ubroLQaPCGFXD76dNi2CxAIDJYNBdCzIVdnZ2gDGq3zBG1eDAZLpbWtJdBTItpDXa09NDdyFoVDhggpBew4t6/YcxqoYHfX2CJ0+qenvpLgSZEHJRj61RfYYxOl5pAQGSgYH/19BQir/QSIewNar/MEbV4GBuDgAShYLuQpAJwdao/sMYVYO9uTkAdGOMIh3C1qj+wxhVA7ZGke4ZU2t0RtkMYa+Q7io0D2NUDSRGuwcG6C4EmRADnTd6Vnx2WdUy2yJbtxK3lx+8/GPPj3RXpEUYo2rAi3qke4Y4bzRVlPpm7ZtvubxVHVhdHVj9psubgjYB3UVpEU6/VwNe1CPdM7i+UalSym3gfj7tc44rh2zZytq6hbWF3qq0ClujarBnMgFjFOmWwfWNFvQWiBXiWKfYwRsZYMz3T2OMqgFbo0j3DK412qZoszOzm2I+he5CdAdjVA0Yo0j3DKg12q3sftL/xMXcpUfZ06Xoorsc3cEYVYO1mZkFgyFXKuW4SCvSFUNpjeZ25y66uyiuJm6R7SKWOeubjm8Gv0uBMf/JYIyqxx7nPCHd0v8JTx2Kjrdq31p9b3W1rLpT0SlVSnnTeH9s+ONR0dGm/iaxQnxWfHZ33W66y9QiHKlXj4O5ecfAgEShcLawoLsWZBL0fMLT953fv1v37uP+x9Zm1typ3P0e+y0ZlgmuCS5Ml380/yPxcaKdmd1Su6UfenxId6VahA8RUU9UXNyNq1f/e/780qVL6a4FmYQffvghMjLS0tLy66+/jo2NZejNkuEtAy1/fPzH4+3HAWCF3YqjvkfnW8+nuyh64EW9emStreInTzo7O+kuBBk/sVj8zjvvrFu3DgDkcnlcXNzixYsvX75Md10AAOkd6YEVgcfbj9ua2fKm8XL9c002QwFjVF0sFgsAxGIx3YUgDWOz2SdOnFC9PHPmzLx582isJzMzMygo6OjRo0wm889//vPJkyf9/PwKCwtffPHF8PDwGzdu0FVYU3/T7x7+LuZRjGhA9Bv73xTNL+JO5ZozzOmqRx9gjKoHYxRpW1NT09atWzdt2vT48eOwsLCioqJDhw7FxcVVVlby+Xw3N7cbN26sXLkyJiamurpal4VRQKWKUudVzDsvPs8yZ/F9+Nf8r82xmqPLGvQTxqh6SIziRT3SBoqiUlNT582bd/bsWUdHx+Tk5Nzc3ICAAPKupaUlh8OpqqricrnW1tbp6emBgYEJCQlNTU06qO2h7OG66nUJdQldiq4ox6jSgFKOK8e4700aP4xR9Tg6OgLGKNKC6urqtWvXJiQkdHV1RUVFlZaW7tmzx8xs6F+ok5MTj8errq7mcDgkdmfPnr1v3z7t/U4OKAc+b/6cfZf9g+SHqRZT0/zSMmdlelt4a+l0BolC6khJSQGA3//+93QXgjQsMDDQ3t7e5SkHB4e5c+fq5tRyuZzH41lZWQGAh4eHQCAY5wfv3r0bHR1Nxu5dXFx4PJ5UKtVsbaWlpUuXLV14ayEUQPTD6Nb+Vs0e3zhgjKrn2LFjALB9+3a6C0EaFhgYmJKSUv/U4cOHSYxevHixoKBAe+fNz88nl+0MBiM+Pr6trU3dI9y+fXv16tWkVTR9+nQ+nz8wMDD5wmQyWVJSkqWlJQAsf235xc6Lkz+mscIYVU9GRgYAvPzyy3QXgjQsMDDw+PHjqpenT5+eO3euRCLx9PRkMBjR0dG1tbWaPaNYLE5MTCSX7bNnz/7hhx/G3r+xsXGMd7Ozs0NCQkiYBgQEpKWlTaa2n3/+eeHChSTcORxOZ2fnZI5m9DBG1XP9+nUACA8Pp7sQpGGjxeif/vQncrltb2//6aefauqqOSMjY/r06QBgYWHB5XL7+vrG3r+qqsrS0jI6Ovr+/fuj7aNQKNLS0mbOnEnCdPny5bm5ueoW1tvby+Vyzc3NAWDmzJnPDXdEYYyqq7i4GADYbDbdhSANGzFGyde1tbXx8fEkm7y9vQUCgVKpnPCJmpqaVEdbvnx5WVnZeD4lEAhImltaWu7Zs6elpWW0PWUyGZ/Pd3d3J6eIjIwsKSkZZ215eXlz584FADMzMw6H093dPc4PmjiMUfXU1taSvyW6C0G6dvXqVdVV8wsvvHDz5k11j6BUKgUCgbOzMwDY2tryeDy1OjHr6uo4HA5pJ9rZ2XG53DGutSUSCY/Hc3BwIJkYHR396NGjMQ7e2dmp6mFgs9m3b98ef2EIY1Q9Dx8+BAAbG5uGhga6a0G6plAoBAKBh4cH6TRUq8O0urp6zZo1JIVffvnlCfe0lpeXR0dHk+O4urryeLwxOgRaWlq4XC4ZJrK1teVyuR0dHcN3u3jxoo+Pz/h7GNAQGKPP9+DBA4FAwOFwAgICGAyGmZkZaRF4enpGRUUlJSVlZGS0tuJEEFPR3d2dlJRkbW2tyiaJRDLG/mQ+E9l/6tSp45/PNIabN2+uXLmShKmvry+fz1coFKPtXFVVpZoX5ezszOPxent7yVsdHR0czi9PTAoJCSksLJx8bSYIY3QEnZ2dly9fPnDgwPr168l8exVbW1vSRTVlyjPPSGAwGP7+/tu3bz915Ah18yb19NcUGau6urr4+HiSTWN0mN68eTMwMFA1n0kkEmmwhuzs7AULFpDfQDabnZmZOcbOP/30k6o57O3tzefzz5075+XlRa6u1O1hQINhjP6itq8vq77+3XffXbBgAWlsqvj4+Gzbti0lJeXOnTv9/f1kGd2cnJyGhoaMjIxNmzbZ2NjY2tqSnXeuXEkBUEwmFRBAxcdTyclUXh4lkz1zslOnqOBgysqKcnen3nmHUn+qINITt27dWrZsGfnRL1my5MaNG6q3enp6VEPes2bNys7O1kYBZHR+xowZqgEloVA4xv5ZWVmqHl7yb8Dq1avHGP1H42G6MdqvVFb09JxqbuY+eLCuqChUKFxeUGA7ZQoAMJnMgIAADocjEAgePnw45IPktvovv/ySvCRDunK5vLCw8MiRI1f376eCgihzcwrg1//s7KiVK6m9e6nqaurIEYrFok6fpsRiqqKCWr+eWriQwt4og0UGjjw9PVUdpjU1Nd999x3pbWQymYmJidoe8pbJZMnJyU5OTqoa7t27N9rOCoXixIkTLBbL2tr673//+2RmHSDCOGM0urw8VChU/Vfe00O2N8vl2e3t/1NX93/u3l1aUDB4n3XFxXvv3z9y6lR+fv7YcwPJVJL9+/eTl4Nnxvyqu5vKy6OSk6n4eCoggGIwfsnTW7coJyfqq69+3bOnh/LwoPh8DX77SPe6urr+8pe/kA5QMqRDBvSLi4t1VkN7ezuXy7WxsSGDRRwOZ4wZ+2RiU2Vlpc7KM2JGG6MXRCKZQkH+a+vv/8vDhxtKSgbn5pKCgm3l5bza2u9Fonp1GoPTpk0DgG3btpGXI8foEK2t1MWL1McfU1evUgBUV9cz73I41Nat6n2HSC/V1dVFR0c7OTnZ2dnR1dtYX1/P4XCYTKZqXpRYLB6+G+mxHeesVTQ243yISExFxRseHhucncnLAYpaXVTUp1TampvPsbEJsbcPtrcPsbefYj6RtWZnzpz56NEjJpNJRp/kcrmXl1dlZeW4Pnz+PLz+Ogx5PNmBA3D9OqSkgI8PPDuihQxOTU2Nn5+fr69vTU0NjWWUl5fv37+f3LscFhaWn58/ZIeQkJDi4uKioqLg4GA6CjQqJvFIOyaD8amf33Qrqxk2NpNfGdDCwgIAWCzWzz//DAAZGRlk2adxcXWFnh6QSMDB4deNTU3g6gqvvw4PH0JcHCQmQmDgpMtE9CBjSkqlkt4yAgMDL1y4cPv2bS6X+9577w3fgTRXB/AZt5pgtOuN8mpr1xYXry0ufruqCgBWsVgzNZGhAGBhYcFkMtva2pydnb29vZ2ftnnHJTQUWCxIT/91S28vZGTA8uXg5gY9PZCaCkFBsGEDXLoExnihYPRIjOpJPC1duvTatWtbtmwZ/hbGqAYZbYz+ftq00/Pnn54//7OnKzVokKurK0VRE7lqs7GBTz+FvXvhm2+gsxMqK2HLFnB3h927IScHKishMRFsbSErC9avB39/OHQI8IElBoXEk0KhoLuQX434MFGMUQ0y2hh1YDLdLS3dLS218UB5Nzc3AMjJyWlpaVH7w7t3w7/+BQcPgrs7rFwJ06bBlStgbQ0A4O8P//u/0NAAycng6wv378O+feDrCwkJUFWl6W/iGfr2QDfDRVqjw2N01apVs2bNevz4MR1FjQBjVIOMNka1p6ysbMWKFQCwZ8+eY8eOxcbGjnd8SWX7diguBpkMWlvh6FFwcXnmXUdH2LMHHj6EjAyIjISuLkhNhYAAWLcOMjPxSl/PjRajDQ0NDx8+lMlkdBQ1AoxRDcIYnQg/Pz/yRX9/v7bOYWYGGzdCdjYUFMCOHWBhATk58MorH7z88ldffdXX16et86LJGS2eRotXumCMapBxxmhaQMAGtUZ+1KSKUblcrr2z/GLRIhAIoLYWPv74ydq1yVlZb7/99vTp0z/88MOGhoZJHru9vf3ChQtGOemNLqPFJcaoETPOGNU2XbRGh5g6FT76yPn779PS0lasWCESiQ4ePOjj47Nx48acnBy1jtTb25uTk7Nv377w8PCpU6e++uqr5eXlALBr1y7Xp1Sr/iB1jRaX+hZb+laPQcMYnQjVShC6aI0OQh4jcePGDaFQGB8fb25u/t13361bty40NDQ1NVUqlY72QZlMdv369aSkpPDwcBaLtW7dukOHDt24cYPJZEZERPT29gLAwYMHi57i8Xg6/LZGZYgDX3hRb4JMYvq9xrm5uVlZWclkMolEQksBoaGhx44d++yzz44cOfLFF18UFhYmJCR89NFHb7zxxnvvveft7Q0ACoWiqKgoJycnPz8/Nze3q6uLfNbc3Dw0NDQyMjIsLGz16tUOT28EcHJyIh8EAPUmw6JBzMzMGAwGWe9j8EwjjFEjhjE6QS4uLo2NjW1tbTTW4OHhceDAgX379p08eTIlJaWkpOTQoUPJyckREREWFha5ubmdnZ1kTwaDsWDBgjVr1qxZs2bVqlVDFktFmmVubj4wMKBQKEhUEfoWW/pWj0HDGJ0gNze3xsbG9vZ2ugsBa2vrnTt37ty5Mz8/PyUl5dy5c+Xl5fX19QAwc+ZM0upcu3YtWVEF6cCIMYqtUSOGMTpBU6dOBQCxPt1iFB4eHh4enpubu2rVKhaLVVJSQh7hOx5lZWWDX8bGxsbGxmqhRrXt2rXr/fffJ1+TVWDorWc8RkxMjFEjhkNME0S6FFtbW+kuZKiVK1c6OTmJxWILLdy+pXt6OPD1XCMmlF7daw/6V49BwxidiPT09MuXLzs4ONTX1+/bt4/25XwIMohEukEBoKioiO6KNIAMfBGGMvA1YsNT3+61x9aoBmGMqqehoWHDhg0xMTESicTX19fMzOzQoUOvvvqqajCHLiEhISwWi/TVLly4EIwlRg0RXtSbGoxRNaSnpwcHB2dlZbFYLD6fX1pamp2d7e7unpmZuXjx4tLSUhprs7OzoyiKLIFKYpR8begoijp9+vS3335LdyFqGOOiHmPUKGGMjktzc/PmzZtjYmLa2tpeeumlsrIycp9PRESEUChcsmTJ/fv3V6xY8c0339BV4eAWKHn0oxHEaFlZ2WuvvRYXF7d9+3YAmMgqMHQY46Jef2JL3zoZDBrG6POlp6ez2exz5845Ojry+fyLFy8Onjw0ffr03Nzct956q7u7OzY2NiEhQXd3iA4yuAUaEBBgbW19//591ZR7w0X+ZxrWcBle1JsajNGxtLa2RkdHx8TEiESiF198UdUIHcLa2vqrr77i8/mWlpapqamRkZFPnjzRcamDW6BMJjMwMJCiKHr7GTTCEGMUL+pNDcboqL7//vvg4OBvv/12ypQpfD4/KytLda/kiDgczpUrVzw9PXNzcxcvXnzr1i2dlQoAbDbbwsKiqqqK3CBvNN2jZNUCw4pRg7io17dYN2gYoyMQi8UJCQlRUVFNTU2RkZGlpaUcDmfEJzEMERYWVlxcHBER0dDQsHr16i+//FIH1RJWVlbz589XKBRkIr3RdI8aYmsUL+pNDcboUFlZWWw2OzU11dbWlsfjXbp0ycfHZ/wfd3Nzu3z5MpfLlclkHA5nx44dYyy8pFmDW6BGM+fJEGMUL+pNDcborzo7OxMSEjZs2NDQ0BAWFlZUVMTlcs3M1P5fxGQyeTzeiRMnbG1tjx8/Hh4eXltbq42ChxjcAg0ODjYzMysrK9PxUn4aZ4gxahAX9fpWj0HDGP3FpUuXgoKCUlNTbWxseDxebm7unDlzJnPA7du337hxw8/Pr7CwcPHixeourjwBg1ugdnZ2c+bMkcvld+/e1fZ5tcpoYpSu1qhSqezo6Bi+HWNUgzBGAQCSkpJeeuml+vp60rk5sUbocCEhIXfu3HnxxRdFItH69esPHTqk1cd1LFy4kMFglJSUkL8N4xhlMtwY1YeL+tLS0rCwsLi4uOFvYYxqkMnF6IgLqq9bt440Qq9fvz7JRugQLi4uWVlZPB6Poqh9+/bFxcX19PRo8PiDTZkyxc/PTyqV3rt3D55e4xt69yiJUUtLS7oLUcOIM9s/+eSThoaGN998Uzc19PX1/fWvfw0NDb1161ZZWdnwx3ZhjGoQLpQHABAeHl5TU0OePq9xDAaDy+UGBATs2LHjzJkzRUVF//nPf+bPn6+Nc23fvl0sFltbWwO2RukzYsPT1dVVZwXk5+e/8847lZWVDAaDw+H84x//GL5QN8aoBplca3Q0WspQlY0bN96+fTswMLCysnLZsmXnzp3Txln+9re/paSkzJw5EwAcHR0BoKCgQK8WRVWX4cYoLQklFov37NmzatWqysrKOXPmXLlyhc/nj/iwgzt37gDAjz/+qPtbRYwPxqju+Pv737p1a+vWrV1dXVu2bNHGCnuNjY3p6ekJCQl+fn7Lli0DAEtLSx8fn4SEhKqqKs2eSzcMMUbpul09MzMzKCgoJSXF3Nycy+WWlpauXr16+G41NTUbNmw4cuSIg4PD3bt3Z82atW/fvhGHodB4USYmMDDQ3t7e5SkHB4e5c+fqsgClUpmcnEz+0jZs2NDe3j7JA9bV1QkEgjfeeEP1vFLC3d1906ZNJEwBwNzcfMuWLfn5+Rr5LnQmKysLANavX093IePV19fn7+8PAHv37h0YGNDNSRsbG7ds2UJ+0GFhYeXl5SPuplAo+Hw+WXGcxWIdOHAgOjqa3Ffi4ODA5XI7Ozt1U7CRMcUYTUlJqX/q8OHDOo5R4urVq+7u7gAwe/bskpKSsXcODAw8fvy46uXp06dnzZqVlpaWmJgYGho6ODodHBwiIyN5PJ5QKCQPp6QoqqioiMPh2NjYkH0WLVokEAjkcrkWvz3NycjIAICNGzfSXci4VFRUkJE9Ozs7AJg3b55AIFAoFNo7o1KpVF22Ozo6Jicnj3a64uLipUuXkt+BqKiox48fk+0//fRTVFQU2e7q6srj8aRSqfYKNkqmGKNDIomWGKUoqq6u7oUXXgAAa2vrr7/+eow9Sc0tLS0ZGRlcLtfPz29wdNrb26uic4y/2Obm5qSkJNVAh4eHR1JSUltbmxa+M006e/YsAGzevJnuQp5PIBCQ9JwxY8Ynn3yimvLBZrPT0tK0ccbS0tLly5erkrG+vn7E3Xp7e5OSkshsBy8vr7Nnzw7fJz8//ze/+Q051PTp0/l8fn9/vzZqNkoYo7TFKEVRUql0586d5HeXw+EMbyF2dXVlZ2e7urr6+fkNvqnfzMyMRGdeXp5a7cq+vj6BQMBms1URzOFwKioqNPptadKZM2cA4LXXXqO7kLG0tLRs2rSJ/C+Nj4+XSCQURcnlcj6fr1pTcfny5VeuXNHUGeVyOY/Hs7KyAgBPT8/09PTR9rx+/fq8efPI7wyHw+nq6hrjsNnZ2WSCBwDMnTtX201po4ExSmeMEmSFPQBYuXJlU1OTRCL57rvvdu7cOXv27MHRaWtrGxkZmZSUtH//fn9//0meNC8vLyoqihyfhHJGRoZGvh3NOn78OAC8/vrrdBcyquzsbPK8UkdHx5MnTw55VyaT8fl88hxZAIiMjBQKhZM8Y15eHpkwx2Aw4uPjR7uk6OjoUC2pExQU9OOPP47n4EqlMi0tTQdNaWNicjGqn/Lz88mfoq2tLTzL39/f3d3dxsZGG8NiVVVViYmJqm7TkJAQPp+vV11jqampehujfX19qhveIiIiRrumpihKIpHweDwyBY3BYERFRRUXF0/gjGKxODExkZyRzGcabc+MjAzSELa2tk5KSpLJZGqdSC6XCwQCVffRsmXLNNiUNj4Yo/qipaXFy8vL3t6e/KU5OzsvWbIkISHhypUr2h4Wa25uPnDgABnyAoC9ezc2NR3s76e521QoFCYmJtra2jo6OlpbWycmJj569IjekgZTjSYxmcykpKTxDMqLRCIul0v+0TIzM4uOjr5///74z1hYWOjp6QkAVlZWBw4cGC0ZGxoafve735Ef5cqVK+/evTv+UwwxvCl9586dCR/NiGGM6hHSIP3vf/87pENKNx0RUqn06NGjISHBly5NFwqhoMDq0aN4qXTkqTPaU1FR8eGHHw6ZvEUwmczo6Gja52yRwXFy3eDn53fz5k21Pt7c3Mzlckm3poWFBYfDaWxsHM8He3p6Zs6cuWLFinHOZxpj1F4tXV1dBw4cIJMBGAzG9u3bKb3sHKMRxqi+qKioIMMFw9/S8a+sRJJXXR0lFDKEQhAKze7dixSLMyhKqb0zUhQlEon4fH5YWJgqNL29vclKr/Hx8YWFhRwOh9zkCgCLFi2iq/OhpaVl48aNQ0aTJqCmpobD4ZD7nWxtbblc7nhmENfW1o6WjCUlJSPOZ9KUtrY20pT+wx/+QGGMPgtjVF988cUXABAXFzf8LVp+Zfv67tXVJRYW2gqFIBRCeXlwaytfodBwcvX0SI4fP75+/XoSKKQZ9fbbb1+7do08O8DV1bWlpYXs/OTJEx6PR9rs8HTOVmtrq2ZLGsPly5fJZTWLxTp9+vTkD1heXj75CfBSqfS585k0paGhgfwPxxgdDGNUX8TExABAamoq3YU8Y2BA3NycXFLiTcK0uHhqQ0NSf//kk0shkeTV1HCKimZZWJgDgKWlZVRUlEAg6OnpoSiqqanJyckJAIanFZmzFRQURMLUysoqPj6+tLR00iWNRSqVqkaT1qxZM8Zo0gTcvn17whPgc3NzyXwmsgqJzm5DwhgdDGNULyiVStKRf+/ePbprGYFSKWtvT7t79wUSpqTbtLe3bAKH6um5U1f3fnGxBzmUUMj45JNt//73v4fM2iG3Nr700ktjHCovLy86OlrVjA0LC8vIyFDdu6VB5eXlwcHBqtEkLU2lVHcC/OD5TGw2e5zzmTSF9puq9QrGqF4gz6Hz8vKiu5DnkEjyHjyILigwJwk4/m5TmayuuTm5oiLkaXpCWdm8hoakvr7q4TtnZmaSi9y6urrnHvn+/fuJiYlkhgMAzJ49Ozk5ubu7eyLf3jCDR5PmzZtXUFCgkcOOYfAEeD8/Pz6fP2Jqq+YzWVhYcLncvr4+bRc2hJ7cVK0nMEb1wr/+9S+9nR05nFRaWVv77qBu06C2tlPd3T9WVq4sLLT/+WdWZeWKrq5rFEUNDHSIRIJ79yKfDlhBUZFzTQ1HIskbLXzFYjEJiMOHD4+/pM7OzuTkZF9fXxJAjo6OiYmJtbW1k/k2m5ubVdfa8fHxmorm5yIT4Mn6JjBsAnxDQ8PmzZvJW+Hh4XTdgYYX9YNhjOqFrVu3AsDRo0fpLkQNT7tNpwuFUFeXWFTk0tT0fwcGOvv728Ti7yWSfIqiqqs3kPQsLLR7+HC7WHxRqXzOndonT55kMBhhYWETuHZWKBQZGRmRkZEkZczMzKKioiY2QerSpUtkNMnV1fX8+fMTOMIk9ff3D5kAn5OTo5rPNPYqJDqAMToYxij9lEolmfqu1mRsPaFUykSiY11dPwiFMHwcv63t1L17L4pExxQKNSYG5ebmTmbSOEVRBQUF8fHxqlVKQ0NDx7+olVQqTUxMJH2Oa9eu1fjMIbVIpdJ//vOfQ9YU37p16zinmmoPxuhgGKP0Ky0tJdMk6S5k4hQKaXHxtAcPXuvs/O/AwGRXUNWUpqamwYtaeXp6JiUliUSiMT5SVla2YMEC0ueovdEkdXV1dX344Yd2dnaOjo4XLlyguxw0FMYo/VJSUkjvG92FTIpMVlNbu7u0dHZBAbO6+iWZ7PmjQ5ROGjVDFrWytraOj48ffiPQkNGkwsJCzZYxSW1tbQDg7OxMdyFoBPgQEfpdu3YNAEZ83oMBsbT09fE5zGZXBwXVUpSitvZtuiv6hZWV1Y4dO0pLS6ZCm3AAAAUMSURBVMmiVjKZ7Pjx42w2e926dZmZmRRFAUBLS0tUVFRCQkJvb298fLxQKFQNl+sJ8rwZjTz3G2kc/lRoRlFUbm4uAERERNBdi2ZYWHi5ur4tlZbQXchQ4eHhmZmZZWVlCQkJNjY2OTk5mzZtWrBgwQcffBAcHHzx4kU3N7cLFy4cO3aMrL6sVzBG9Rn+VGhWWloqEommT58+ZE17wyKTPWpsPNDXV6lU9slk91tbj9jZLaW7qJEFBAQcOXKksbExOTnZx8enrKzs1KlTT548iYyM/Pjjj7OzswsKCuiucQSk1Tx4/VmkP/A59TQjV/Rr1qyhu5BJMTefIpfXVFev7+9vYjKdp0xZ5+39+Tg/u2vXrvfff598LZfLVbfMa1V8fPydO3euXLlSVFRkZ2cnkUji4uLefPNNgUCwcOHCIU+40gfYGtVnGKM0M46OUSbTZcaMryf22YMHD6rWx8zIyCADbtomFoufPHnS3Ny8bds21UYyUbSpqUkHBagLW6P6DGOUTkqlMi8vDww/RifDycnJ29ubfO3s7Kybk5LEfPLkyeCNHh4ewzfqCWyN6jP8qdCppKREJBL5+PiMuEox0h6SmEMantgaRRODMUonckW/du1augsxOSO2RvU5RrE1qs/wop5OXl5eERERv/3tb+kuhDZkaSuV2NjY2NhYHZx3jNYoXtQjdWGM0ikmJoas1ox0zOBao3hRr8/wHzddY7PZJ06cUL08c+YMWb0c6dKIrVE7Ozt7e3upVNrZ2UlTXaPC1qg+w58KMkWjNTz1tkGKrVF9hjGKTJGrq6uFhUV7e7tcLh+8XW+7R7E1qs/wp4JMkZmZmbu7O0VRzc3Ng7ePeLGvD7A1qs8wRmmwa9cu16c4HA7d5ZioERNzrb//gSVLfCUSmooaFbZG9RmO1NOAltsf0RCHly1zsbFxbmsbvJFjZwd37oD+3VSGMarPMEZpQMvtj2iIpTIZ5OdDXd0zWz09AQDwoh6pA/9xQ6ZqxMQkG3GICakDW6PIVHl4AAxLTLJR/1qj7u7uXC6XTCRA+oZBLhYQMjnnzsHmzfDKK3D+/K8bW1vB3R1cXEAkoq8yZGDwGgGZqhEbnq6uYGEB7e0gk9FS1HB425v+wxhFpmrEvlEGA9zdgaLg2fmkCI0BYxSZKk9PYDCguRmG9Gvp62A90lsYo8hUWVkBiwVyObS3P7NdXwfrkd7CGEUmbMTuUf0brMfb3vQcxigyYSM2PPWvNXrw4MGip3g8Ht3loKFw3igyYaO1Rs3NobublopGhLe96TmMUWTCwsJAJoMhc9p37oSEBMD7hdC4YYwiE7Z7N+zePXSjpSUdpSADhncxIQRw+jQcOgSVleDoCK+8Ajwe4LUzGje8ckEmj8+H3bth3z5oboZr16C+HiIj9ecuJqT/sDWKTFtfH3h5weefw1tv/bKltxdmzYKPPwacWoTGB1ujyLQJhdDRAdHRv26xtYVNmyA7m76akIHBGEWmTSQCOztwcHhmo6cnrvCExg9jFJk2V1fo6YEhD19qagJXV5oKQoYHYxSZttBQYLEgPf3XLb29kJEBa9fSVxMyMDhvFJk2Gxv49FPYuxfs7GD9emhqgg8+AHd3eOMNuitDBgNH6hECOHkSPvsMKithyhR45RU4dAhcXOiuCRkMjFGEEJoU7BtFCKFJwRhFCKFJwRhFCKFJwRhFCKFJwRhFCKFJwRhFCKFJ+f/79N5pnzJC6AAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "offmol = OFFMolecule.from_smiles('c1cc(-Cl)ccc1C(=O)CS[C]1=CO[C](F)(F)CC1')\n", "offmol.generate_conformers(n_conformers=1)\n", "offmol.visualize()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Make the OpenMM system and grab the resulting Topology (with VSites added)\n", "\n", "off_top = offmol.to_topology()\n", "off_sys, off_top_w_vsites = ff.create_openmm_system(offmol.to_topology(), return_topology=True)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There are 5 virtual particles\n", "EP None \n", "EP None \n", "EP None \n", "EP None \n", "EP None \n" ] } ], "source": [ "# Retrieve the number of virtualsites that were added\n", "n_virtual_sites = off_top_w_vsites.n_topology_particles - off_top_w_vsites.n_topology_atoms\n", "print(f\"There are {n_virtual_sites} virtual particles\")\n", "\n", "# Get an OpenMM Topology from our openFF Topology\n", "omm_top = off_top.to_openmm()\n", "\n", "# to_openmm doesn't handle virtualsites yet, so we add them manually to the openmm topology here\n", "for top_particle in off_top_w_vsites.topology_particles:\n", " if isinstance(top_particle, TopologyAtom):\n", " continue\n", " print(top_particle.virtual_site.name, app.Element.getByMass(0), [i for i in omm_top.residues()][0])\n", " omm_top.addAtom(top_particle.virtual_site.name, app.Element.getByMass(0), [i for i in omm_top.residues()][0])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from simtk import openmm, unit\n", "\n", "# Propagate the System with Langevin dynamics.\n", "time_step = 2*unit.femtoseconds # simulation timestep\n", "temperature = 300*unit.kelvin # simulation temperature\n", "friction = 1/unit.picosecond # collision rate\n", "integrator = openmm.LangevinIntegrator(temperature, friction, time_step)\n", "\n", "# Length of the simulation.\n", "num_steps = 1000 # number of integration steps to run\n", "\n", "# Logging options.\n", "trj_freq = 1 # number of steps per written trajectory frame\n", "data_freq = 1 # number of steps per written simulation statistics\n", "\n", "# Set up an OpenMM simulation.\n", "simulation = openmm.app.Simulation(omm_top, off_sys, integrator)\n", "\n", "# Set the initial positions.\n", "positions = np.concatenate((offmol.conformers[0], np.zeros((n_virtual_sites, 3))), \n", " axis=0) * unit.angstrom\n", "simulation.context.setPositions(positions)\n", "\n", "# This is necessary since we placed all virtual sites at 0,0,0. This will cause the first simulation\n", "# step to crash if the positions are not calculated beforehand\n", "simulation.context.computeVirtualSites()\n", "\n", "# Randomize the velocities from a Boltzmann distribution at a given temperature.\n", "simulation.context.setVelocitiesToTemperature(temperature)\n", "\n", "# Configure the information in the output files.\n", "pdb_reporter = openmm.app.PDBReporter('trajectory.pdb', trj_freq)\n", "state_data_reporter = openmm.app.StateDataReporter('data.csv', data_freq, step=True,\n", " potentialEnergy=True, temperature=True,\n", " density=True)\n", "simulation.reporters.append(pdb_reporter)\n", "simulation.reporters.append(state_data_reporter)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting simulation\n", "Elapsed time 6.48 seconds\n", "Done!\n" ] } ], "source": [ "import time\n", "\n", "print(\"Starting simulation\")\n", "start = time.process_time()\n", "\n", "# Run the simulation\n", "simulation.step(num_steps)\n", "\n", "end = time.process_time()\n", "print(\"Elapsed time %.2f seconds\" % (end-start))\n", "print(\"Done!\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now go visualize `trajectory.pdb` in VMD or PyMol!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: Numerical comparison of OpenMM's TIP5P and an equivalent SMIRNOFF implementation " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Parameterizes a water box with OpenFF and OpenMM forcefields. Currently set up\n", "to use a TIP5P definition. The code examines the geometry and energy between the\n", "two, and examines cases where minimization is performed. Specifically, the code\n", "compares the four possible combinations:\n", " - oFF and oMM geometry/energy, minimized separately and then compared\n", " - oFF and oMM geometry/energy using no minimization\n", " - Geometry minimized using oFF, then a single point is done with oMM\n", " - Geometry minimized using oMM, then a single point is done with oFF\n", "\n", "The virtual site definitions give differences in geometry and energy mostly due\n", "to how they were defined from their parent atoms. OpenMM uses an OutOfPlaneSite\n", "definition, whereas OpenFF uses the LocalCoordinatesSite definition (both are\n", "OpenMM types). In the OutOfPlaneSite definition, both angle and distance are\n", "variable since the defined out-of-plane angle depends on a weighted vector cross.\n", "The cross is a function of the O-H vectors, so the virtual sites are sensitive\n", "to the molecular geometry. In the OpenFF version, the distance is fixed to a constant\n", "value, and the out-of-plane angle is explicitly required in the OpenFF spec.\n", "\n", "In this example, the OpenFF parameter definition (the \"offxml\") is a string\n", "further below in the `main` function, and can be easily modified to explore\n", "forcefield parameterization. The OpenMM definition is loaded from its internal \n", "default location, and acts as a reference. One can change this this to a different\n", "filename to compare other forcefields.\n", "\n", "This example is somewhat hardcoded to operate on water molecules, but can be easily\n", "modified to examine other cases as well. The only major assumption is that the\n", "system is homogenous, i.e., all of the molecules are same. The reason this is \n", "assumed is mostly due to the difference in how virtual sites are handled between\n", "OpenFF and OpenMM. Whereas OpenMM interleaves the virtual site particles between \n", "the atomic particles, OpenFF instead aggregates all virtual sites and places them\n", "last. The code below does assume that, barring this difference, the virtual sites \n", "are added in the same order.\n", "\n", "The example begins in `run_tests` by defining a grid of water molecules, with\n", "the default being a single water molecule (Nx=1, Ny=1, Nz=1). From this, the\n", "calculations described in the first paragraph above are performed. The energy\n", "difference and distance between the two geometries, per atom, is then reported. \n", "There are commented lines that print the entire set of coordinates, and can be\n", "uncommented if desired." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "The following functions are utilities and other infrastructure to make comparison between OpenFF and OpenMM calculations easier.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from simtk.openmm import app, LangevinIntegrator, XmlSerializer\n", "from simtk import unit\n", "from openforcefield.typing.engines.smirnoff import ForceField as OFFForceField\n", "from openforcefield.topology import Molecule as OFFMolecule\n", "from openforcefield.topology import Topology as OFFTopology\n", "from openforcefield.topology import TopologyAtom\n", "import numpy as np\n", "\n", "# To get a fair comparison , take the OpenMM virtual sites and put them at the\n", "# end so that they follow the same order as the oFF order.\n", "def reorder_openmm_to_openff(xyz, n_atoms_per_mol, n_vptls_per_mol):\n", " \"\"\"\n", " \"\"\"\n", " # constants for simplicity\n", " n_particles_per_mol = (n_atoms_per_mol + n_vptls_per_mol)\n", " n_mols = xyz.shape[0] // n_particles_per_mol\n", "\n", " atom_base_indices = np.arange(n_atoms_per_mol)\n", " atom_indices = np.hstack([atom_base_indices + (i*n_particles_per_mol) for i in range(n_mols)])\n", "\n", " vsite_base_indices = np.arange(n_vptls_per_mol) + n_atoms_per_mol \n", " vsite_indices = np.hstack([vsite_base_indices + (i*n_particles_per_mol) for i in range(n_mols)])\n", "\n", " mask = np.hstack([atom_indices, vsite_indices])\n", "\n", " return xyz[mask]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "def reorder_openff_to_openmm(xyz, n_atoms_per_mol, n_vptls_per_mol):\n", " # constants for simplicity\n", " n_particles_per_mol = (n_atoms_per_mol + n_vptls_per_mol)\n", " n_mols = xyz.shape[0] // n_particles_per_mol\n", "\n", " atom_base_indices = np.arange(n_atoms_per_mol)\n", " vsite_base_indices = np.arange(n_vptls_per_mol) + n_atoms_per_mol * n_mols\n", "\n", " mask = np.hstack([\n", " np.hstack([\n", " atom_base_indices + (i*n_atoms_per_mol), \n", " vsite_base_indices + (i*n_vptls_per_mol)\n", " ]) for i in range(n_mols)])\n", "\n", " return xyz[mask]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "def openmm_evaluate_vsites_and_energy(omm_top,\n", " omm_sys,\n", " atom_xyz_in_nm,\n", " minimize=False):\n", " \"\"\"\n", " Calculate the Virtual Site positions and potential energy\n", " of the given system\n", "\n", " Parameters\n", " ----------\n", " omm_top: OpenMM Topology\n", " omm_sys: OpenMM System\n", " atom_xyz_in_nm: N,3 list of floats interpreted as nanometers\n", " The coordinates of the molecules\n", " minimize: bool \n", " Perform a minimization before calculating energy\n", "\n", " Returns\n", " -------\n", " pos: N,3 List of floats of the positions of all particles\n", " ene: The potential energy of the particle positions\n", " \"\"\"\n", "\n", " integ = LangevinIntegrator(300 * unit.kelvin, 1 / unit.picosecond,\n", " 0.002 * unit.picoseconds)\n", "\n", " sim = app.Simulation(omm_top, omm_sys, integ)\n", " pos = sim.context.getState(getPositions=True).getPositions()\n", "\n", "\n", " sim.context.setPositions(atom_xyz_in_nm)\n", "\n", " # Once the atom positions are known, virtual sites can\n", " # be determined.\n", " sim.context.computeVirtualSites()\n", "\n", "\n", " if minimize:\n", " sim.minimizeEnergy()\n", "\n", " state = sim.context.getState(getEnergy=True, getPositions=True)\n", " ene = state.getPotentialEnergy()\n", " pos = [\n", " list(xyz) for xyz in state.getPositions().value_in_unit(unit.nanometer)\n", " ] * unit.nanometer\n", " return pos, ene" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "\n", "def build_openff_water_lattice(Nx=1, Ny=1, Nz=1, dx=2.0, dy=2.0, dz=2.0):\n", " \"\"\"\n", " Generate a box of water molecules as OpenFF Molecules\n", "\n", " Parameters\n", " ----------\n", " Nx, Ny, Nz: The number of molecules in each directon\n", " dx, dy, dz: The spacing between the molecules in each direction\n", "\n", " Returns\n", " -------\n", " water_box: A list of OpenFF molecules with a 3D conformation\n", " \"\"\"\n", " Lx, Ly, Lz = Nx * dx, Ny * dy, Nz * dz\n", " Z, Y, X = np.mgrid[0:Lz:dz, 0:Ly:dy, 0:Lx:dx]\n", " XYZ = [list(xyz) for xyz in zip(X.flat, Y.flat, Z.flat)]\n", "\n", " water_box = list([None] * len(XYZ))\n", "\n", " water_reference = OFFMolecule.from_smiles('O')\n", " water_reference.atoms[0].name = 'O'\n", " water_reference.atoms[1].name = 'H1'\n", " water_reference.atoms[2].name = 'H2'\n", " water_reference.generate_conformers()\n", "\n", " for i, xyz in enumerate(XYZ):\n", " water_box[i] = OFFMolecule(water_reference)\n", " water_box[i].conformers[0] += xyz * unit.angstrom\n", "\n", " return water_box\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "\n", "def insert_vsite_padding(xyz, n_mol, n_vptl, order=\"OpenFF\"):\n", " \"\"\"\n", " Given a set of atomic coordinates, insert dummy coordinates in indices\n", " expected for either OpenFF or OpenMM. OpenFF places the dummy coordinates\n", " at the end, whereas OpenMM expects virtual sites to be next to their owning\n", " molecule.\n", " Parameters\n", " ----------\n", " xyz: N,3 List, of coordinates of atoms only\n", " n_mol: int, The number of molecules\n", " n_vptl: int, The number of virtual site particle dummy positions to insert\n", " \n", " Returns\n", " -------\n", " xyz: (N+n_vsite,3) list, the coordinates of all particles in the\n", " expected order\n", " \"\"\"\n", "\n", " # heuristic that vsites come at the end; applies to OpenFF\n", " vsite_pad = np.arange(n_vptl*3).reshape(n_vptl, 3)\n", " xyz = np.vstack([xyz, vsite_pad])\n", "\n", " # but does not apply to OpenMM. Need to interleave vsite coordinates\n", " if order == \"OpenMM\":\n", " n_vptls_per_mol = n_vptl // n_mol\n", " n_atoms_per_mol = len(xyz) // n_mol - n_vptls_per_mol\n", " xyz = reorder_openff_to_openmm(xyz, n_atoms_per_mol, n_vptls_per_mol)\n", " \n", " return xyz\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "\n", "def coords_from_off_mols(mols, conformer_id=0, unit=unit.angstrom):\n", " \"\"\"\n", " Small utility to extract the coordinates from the molecule conformers,\n", " with optional unit conversion\n", " \"\"\"\n", " xyz = np.vstack(\n", " [mol.conformers[conformer_id].value_in_unit(unit) for mol in mols])\n", " return xyz * unit\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function to calculate water using the oMM definition\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "\n", "def evaluate_water_omm(water, ff, minimize=False):\n", " \"\"\"\n", " Given a list of molecules and a forcefield definition, calculate the\n", " positions and energy.\n", " \n", " Parameters\n", " ----------\n", " water: List of OFFMolecules, each with a 3D conformation\n", " ff: an OpenMM ForceField object\n", " minimize: Boolean, whether the structure should be minimized\n", "\n", " Returns\n", " -------\n", " xyz: List The coordinates of all particles in the system (OpenMM ordering)\n", " ene: float, The potential energy\n", " \"\"\"\n", "\n", " top = OFFTopology.from_molecules(water).to_openmm()\n", "\n", " # get the oFF topology and the positions\n", " atom_xyz = coords_from_off_mols(water, unit=unit.nanometer)\n", "\n", " # first pass; no virtual sites\n", " omm_modeller = app.Modeller(top, atom_xyz)\n", " omm_modeller.addExtraParticles(ff)\n", "\n", " # second pass: virtual sites now present (according to oMM FF)\n", " top = omm_modeller.getTopology()\n", "\n", "\n", " sys = ff.createSystem(top, nonbondedMethod=app.NoCutoff)\n", " n_vptl = sys.getNumParticles() - len(atom_xyz) \n", " n_mol = len(water)\n", "\n", " # need to insert dummy positions into coordinates, since OpenMM\n", " # systems create virtual site particles in place, i.e. OHHMMOHHMM etc.\n", " atom_xyz = insert_vsite_padding(atom_xyz, n_mol, n_vptl, order=\"OpenMM\")\n", "\n", " # Need an openMM topology, an openMM system, and a set of positions\n", " # to calculate the energy. Also returns the vsite positions based\n", " # on supplied atom coordinates\n", " xyz, ene = openmm_evaluate_vsites_and_energy(top, sys, atom_xyz, minimize=minimize)\n", "\n", " return xyz,ene" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function to calculate water using the oFF definition\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def evaluate_water_off(water, ff, minimize=False):\n", " \"\"\"\n", " Given a list of molecules and a forcefield definition, calculate the\n", " positions and energy.\n", " \n", " Parameters\n", " ----------\n", " water: list of OFFMolecules, each with a 3D conformation\n", " ff: an OpenFF ForceField object\n", " minimize: boolean, whether the structure should be minimized\n", "\n", " Returns\n", " -------\n", " xyz: list The coordinates of all particles in the system (OpenMM ordering)\n", " ene: float, The potential energy\n", " \"\"\"\n", "\n", " # get the oFF topology and the positions\n", " top = OFFTopology.from_molecules(water)\n", " atom_xyz = coords_from_off_mols(water, unit=unit.nanometer)\n", "\n", " # IMPORTANT! The new system has the virtual sites, but the oFF top does not as there is no API yet to do so\n", " sys = ff.create_openmm_system(top, allow_nonintegral_charges=True)\n", "\n", " n_vptl = sys.getNumParticles() - len(atom_xyz) \n", " n_mol = len(water)\n", " atom_xyz = insert_vsite_padding(atom_xyz, n_mol, n_vptl, order=\"OpenFF\")\n", "\n", " # Need an openMM topology, an openMM system, and a set of positions\n", " # to calculate the energy. Also returns the vsite positions based\n", " # on supplied atom coordinates\n", " xyz, ene = openmm_evaluate_vsites_and_energy(top.to_openmm(), sys, atom_xyz, minimize=minimize)\n", "\n", " return xyz,ene\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "\n", "def print_info(xyz, ene, name, crd_units=unit.angstrom):\n", " \"\"\"\n", " \"\"\"\n", " print(f\"Results for: {name}\")\n", " print(f\"Energy: {ene}\")\n", " print(\"Coordinates:\")\n", " print(xyz*crd_units)\n", " print()\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "tip5p_offxml = '''\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "'''\n", "\n", "constraints = \"\"\"\n", " \n", " \n", " \n", " \n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Our custom-crafted oFF force field definition above\n", "\n", "off_ff = OFFForceField(tip5p_offxml)\n", "\n", "# The standard OpenMM definition of tip5p\n", "omm_ff = app.ForceField('tip5p.xml')\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "===> Evaluate OFF and OMM no minimization single points <===\n", "Nxyz = (1, 1, 1) dxyz = (3, 3, 3)\n", "Results for: OFF - OMM norm\n", "Energy: -3.019806626980426e-14 kJ/mol\n", "Coordinates:\n", "[ 0.0000000000 0.0000000000 0.0000000000 0.0256116359 0.0256116359] A\n", "\n" ] } ], "source": [ "\n", "minimize=True\n", "\n", "\n", "Nx, Ny, Nz = (1,1,1) # 2x2x2 = 8 water molecules\n", "dx, dy, dz = (3,3,3) # water spaced 3A apart in each direction\n", "\n", "#Nx, Ny, Nz = (1,1,1) \n", "\n", "gridspec = [Nx,Ny,Nz,dx,dy,dz]\n", "\n", "n_atoms_per_mol = 3\n", "n_vptls_per_mol = 2\n", "np.set_printoptions(formatter={'float_kind': \"{:13.10f}\".format})\n", "\n", "if False:\n", " minimize=True\n", " print(\"===> Evaluate separate OFF and OMM minimizations <===\")\n", " print(f\"Nxyz = {(Nx,Ny,Nz)} dxyz = {(dx,dy,dz)}\")\n", "\n", " # default spacing of water is 2 Angstroms (dx,dy,dz arguments)\n", " water_box = build_openff_water_lattice(*gridspec)\n", "\n", " off_crds, off_ene = evaluate_water_off(water_box, off_ff, minimize=minimize)\n", " off_crds = np.array(off_crds.value_in_unit(unit.angstrom))\n", "\n", " omm_crds, omm_ene = evaluate_water_omm(water_box, omm_ff, minimize=minimize)\n", " omm_crds = np.array(omm_crds.value_in_unit(unit.angstrom))\n", "\n", " # constants for ordering (to compare positions)\n", " # We assume here that the system is homogenous for the reordering\n", " omm_crds = reorder_openmm_to_openff(omm_crds, n_atoms_per_mol, n_vptls_per_mol)\n", "\n", " # print_info(off_crds, off_ene, \"OpenForceField\")\n", " # print_info(omm_crds, omm_ene, \"OpenMM\")\n", " # print_info(off_crds - omm_crds, off_ene - omm_ene, \"OFF - OMM\")\n", " print_info(np.linalg.norm(off_crds - omm_crds, axis=1), off_ene - omm_ene, \"OFF - OMM norm\")\n", "\n", "if True:\n", " minimize=False\n", " print(\"===> Evaluate OFF and OMM no minimization single points <===\")\n", " print(f\"Nxyz = {(Nx,Ny,Nz)} dxyz = {(dx,dy,dz)}\")\n", "\n", " # default spacing of water is 2 Angstroms (dx,dy,dz arguments)\n", " water_box = build_openff_water_lattice(*gridspec)\n", "\n", " off_crds, off_ene = evaluate_water_off(water_box, off_ff, minimize=minimize)\n", " off_crds = np.array(off_crds.value_in_unit(unit.angstrom))\n", "\n", " omm_crds, omm_ene = evaluate_water_omm(water_box, omm_ff, minimize=minimize)\n", " omm_crds = np.array(omm_crds.value_in_unit(unit.angstrom))\n", "\n", " # We assume here that the system is homogenous for the reordering\n", " omm_crds = reorder_openmm_to_openff(omm_crds, n_atoms_per_mol, n_vptls_per_mol)\n", "\n", " # print_info(off_crds, off_ene, \"OpenForceField\")\n", " # print_info(omm_crds, omm_ene, \"OpenMM\")\n", " # print_info(off_crds - omm_crds, off_ene - omm_ene, \"OFF - OMM\")\n", " print_info(np.linalg.norm(off_crds - omm_crds, axis=1), off_ene - omm_ene, \"OFF - OMM norm\")\n", "\n", "if False:\n", " print(\"===> Evaluate OMM minimum with OFF single point <===\")\n", " print(f\"Nxyz = {(Nx,Ny,Nz)} dxyz = {(dx,dy,dz)}\")\n", " \n", " water_box = build_openff_water_lattice(*gridspec)\n", " omm_crds, omm_ene = evaluate_water_omm(water_box, omm_ff, minimize=True)\n", " omm_crds = np.array(omm_crds.value_in_unit(unit.angstrom))\n", "\n", " # Need some machinery here to copy the minimized coordinates into the next single\n", " # point calculation\n", " omm_crds = reorder_openmm_to_openff(omm_crds, n_atoms_per_mol, n_vptls_per_mol)\n", " for i,mol in enumerate(water_box):\n", " pos = omm_crds[3*i:3*i + 3]\n", " mol.conformers[0][:] = pos * unit.angstrom\n", "\n", " off_crds, off_ene = evaluate_water_off(water_box, off_ff, minimize=False)\n", " off_crds = np.array(off_crds.value_in_unit(unit.angstrom))\n", "\n", " # print_info(omm_crds, omm_ene, \"OpenMM\")\n", " # print_info(off_crds, off_ene, \"OpenForceField\")\n", " # print_info(off_crds - omm_crds, off_ene - omm_ene, \"OFF - OMM\")\n", " print_info(np.linalg.norm(off_crds - omm_crds, axis=1), off_ene - omm_ene, \"OFF - OMM norm\")\n", "\n", "if False:\n", " print(\"===> Evaluate OFF minimum with OMM single point <===\")\n", " print(f\"Nxyz = {(Nx,Ny,Nz)} dxyz = {(dx,dy,dz)}\")\n", " \n", " water_box = build_openff_water_lattice(*gridspec)\n", " off_crds, off_ene = evaluate_water_off(water_box, off_ff, minimize=True)\n", " off_crds = np.array(off_crds.value_in_unit(unit.angstrom))\n", "\n", " # copy the minimized coordinates, no shuffling needed\n", " for i,mol in enumerate(water_box):\n", " pos = off_crds[3*i:3*i + 3]\n", " mol.conformers[0][:] = pos * unit.angstrom\n", "\n", " omm_crds, omm_ene = evaluate_water_omm(water_box, omm_ff, minimize=False)\n", " omm_crds = np.array(omm_crds.value_in_unit(unit.angstrom))\n", " omm_crds = reorder_openmm_to_openff(omm_crds, n_atoms_per_mol, n_vptls_per_mol)\n", "\n", " # print_info(omm_crds, omm_ene, \"OpenMM\")\n", " # print_info(off_crds, off_ene, \"OpenForceField\")\n", " # print_info(off_crds - omm_crds, off_ene - omm_ene, \"OFF - OMM\")\n", " print_info(np.linalg.norm(off_crds - omm_crds, axis=1), off_ene - omm_ene, \"OFF - OMM norm\")\n" ] } ], "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.7.8" } }, "nbformat": 4, "nbformat_minor": 4 }