rtp_to_top/PDB from RTP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extracting GRO/PDB file from RTP\n",
    "\n",
    "GROMACS natively have force-field parameters for small molecules included in rtp files. However, there are no coordinates of these molecules (PDB/GRO) included. Here is very simple algorithm to create these coordinates.\n",
    "Sketch:\n",
    "1. Extract atom names\n",
    "2. Put atoms on a surface of sphere\n",
    "3. Minimaze energy\n",
    "4. Hope that there are no entanglements (like long alkane crossing in the middle of a ring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Firstly, let's see what are avaible molecules in RTP file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "File = open(\"charmm36-mar2019.ff/merged.rtp\")\n",
    "text = File.read()\n",
    "File.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_residues = []\n",
    "for line in text.splitlines():\n",
    "    if len(line)>0:\n",
    "        if line[0] == \"[\":\n",
    "            all_residues.append(line.split()[1])\n",
    "all_residues=all_residues[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Number of residues:', 1360)\n",
      "['11BPO', '123TBB', '123TCB', '123TIB', '124TBB', '124TCB', '124TIB', '12A', '12DBB', '12DCB', '12DIB', '12MU', '135TBB', '135TCB', '135TIB', '13BPO', '13DB', '13DBB', '13DCB', '13DIB']\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of residues:\", len(all_residues))\n",
    "print(all_residues[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In CHARMM36 there are 1360 parametrized molecules. Let's select one of the residues, for example here we chose residue with name BENZ, which describes benzene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "selected = \"BENZ\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's extract atom names from RTP file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def atom_names(residue):\n",
    "    inside = 0\n",
    "    atom_names = [] \n",
    "    for line in text.splitlines():\n",
    "        if inside>0:\n",
    "            if \"[\" in line:\n",
    "                inside+=1\n",
    "            if inside==2 and \"[\" not in line:\n",
    "                atom_names.append(line.split()[0])\n",
    "        if residue in line:\n",
    "            inside = 1\n",
    "    return atom_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what are the atoms in BENZ:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CG', 'HG', 'CD1', 'HD1', 'CD2', 'HD2', 'CE1', 'HE1', 'CE2', 'HE2', 'CZ', 'HZ']\n"
     ]
    }
   ],
   "source": [
    "print(atom_names(selected))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's radomly sample points from a sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Some things are faster to google then to implement:\n",
    "# http://mathworld.wolfram.com/SpherePointPicking.html\n",
    "# and the implementation:\n",
    "# https://stackoverflow.com/questions/33976911/generate-a-random-sample-of-points-distributed-on-the-surface-of-a-unit-sphere\n",
    "\n",
    "def sample_spherical():\n",
    "    vec = np.random.randn(3, 1)\n",
    "    vec /= np.linalg.norm(vec, axis=0)\n",
    "    return 2*np.array(vec.T[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And add to every atom coordinante of radomly chosen place on a sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def save_sphere_pdb(residue, atom_names):\n",
    "    File = open(residue+\".pdb\",\"w\")\n",
    "    File.write(\"COMPND    {}\\n\".format(residue))\n",
    "    File.write(\"AUTHOR    GENERATED\\n\")\n",
    "    for atom in range(len(atom_names)):\n",
    "        rand = sample_spherical()\n",
    "        File.write(\"HETATM   {:2d}  {:3s} {}    0      {: 5.3f}  {: 5.3f}  {: 5.3f}  1.00  0.00           {}\\n\".format(atom, atom_names[atom],residue,rand[0],rand[1],rand[2], atom_names[atom][0]))\n",
    "    File.write(\"END\\n\")\n",
    "    File.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "save_sphere_pdb(selected,atom_names(selected))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's visualze the results. To do that I use [nglview](https://github.com/arose/nglview) library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import nglview\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nglview.show_structure_file(selected+\".pdb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is definitly questionable structure of benzene. But that's just start. Now we can perform energy minimalization and hope we obtain perfect structure. To do so, we will use [GROMACS](http://www.gromacs.org/) and python [GromacsWrapper](https://github.com/Becksteinlab/GromacsWrapper):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import gromacs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def generate_gro(residue, gro=None,v=True):\n",
    "    if gro==None:\n",
    "        gro=residue\n",
    "    pdb_name = residue+\".pdb\"\n",
    "    \n",
    "    save_sphere_pdb(residue,atom_names(residue))\n",
    "    # Generate topology for guessed PDB:\n",
    "    gromacs.pdb2gmx(f=pdb_name, ff=\"charmm36-mar2019\", water=\"none\", o =residue+\".gro\", p =residue+\".top\")\n",
    "    gromacs.editconf(f=residue+\".gro\",o=residue+\".gro\", box =[5,5,5])\n",
    "    # Minimaze energy:\n",
    "    gromacs.grompp(f=\"em.mdp\", p=residue+\".top\",c=residue+\".gro\",o=gro+\".tpr\")\n",
    "    gromacs.mdrun(deffnm=gro)\n",
    "    # As return, let's visualize (if v true, otherwise let's don't do anything):\n",
    "    if v:\n",
    "        return nglview.show_structure_file(gro+\".gro\")\n",
    "    else:\n",
    "        return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "generate_gro(selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's how you can obtain perfect GRO file. You have to remember that this method works for molecules without any rings and it MIGHT work for other molecules (but no promises) \n",
    "\n",
    "You can also easly convert it to PDB file by command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, None, None)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gromacs.editconf(f=selected+\".gro\",o=selected+\".pdb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nglview.show_structure_file(selected+\".pdb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create many coordinates files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "selected = [\"MEOH\",\"ETOH\",\"PRO2\",\"BENZ\", \"PHEN\", \"OCOH\",\"BALD\"]\n",
    "for select in selected:\n",
    "    generate_gro(select, v=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "widgets": {
   "state": {
    "2ddaf301a8be4456ae97d3cccbda791b": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "adf7920015464c2687fc010a88f1c262": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    },
    "e082b20f584840278f7077dae87cce76": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}