Skip to content

Instantly share code, notes, and snippets.

@iwatobipen

iwatobipen/openmm-psikit.ipynb

Created March 2, 2020 00:37
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save iwatobipen/22e1c2a02045d9d049e66cecc66c80bf to your computer and use it in GitHub Desktop.
Save iwatobipen/22e1c2a02045d9d049e66cecc66c80bf to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
openmm-psikit.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"RDKit WARNING: [09:11:58] Enabling RDKit 2019.09.3 jupyter extensions\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"from psikit import Psikit\n",
"from psiomm import molecule\n",
"from psiomm import psi_omm as po\n",
"\n",
"from simtk.openmm.app import Simulation\n",
"from simtk.openmm.app import ForceField\n",
"from simtk.openmm.app import Modeller\n",
"from simtk.openmm.app import PDBReporter\n",
"\n",
"from simtk.openmm import Platform\n",
"from simtk.openmm import LangevinIntegrator\n",
"\n",
"from simtk.openmm import Vec3\n",
"\n",
"from simtk.unit import picosecond\n",
"from simtk.unit import kelvin\n",
"from simtk.unit import kilocalorie_per_mole\n",
"from simtk.unit import angstrom\n",
"\n",
"from rdkit import Chem\n",
"from rdkit.Chem import AllChem\n",
"from rdkit.Chem.Draw import IPythonConsole"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Work flow \n",
"\n",
"- make molecule from xyz string => solute\n",
"- generate bonds\n",
"- generate atom types\n",
"- generate charges\n",
"\n",
"- make topology\n",
"- define forcefield\n",
"- generate template with FF, topology and solute\n",
"- make modeller\n",
"- perform simulation to get trajectory"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"pk = Psikit()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"pk.read_from_smiles('CCOc1cccnc1')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"pk.rdkit_optimize()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pk.mol.GetNumConformers()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# remove mutiplicity, formal charge\n",
"xyzstring = '\\n'.join(pk.mol2xyz().split('\\n')[2:])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"solute = molecule.Molecule.from_xyz_string(xyzstring)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"solute.generate_bonds()\n",
"solute.generate_atom_types()\n",
"solute.generate_charges()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0 c3 -0.16240415191041357\n",
"1 c3 0.03679384783580364\n",
"2 os -0.263056064493707\n",
"3 ca 0.130662401965024\n",
"4 ca -0.0752094501302798\n",
"5 ca -0.05835160068850964\n",
"6 ca 0.017895212168927088\n",
"7 nb -0.2358970396079183\n",
"8 ca -0.0012324033472115303\n",
"9 hc 0.06626046228444116\n",
"10 hc 0.06408533045243181\n",
"11 hc 0.0662604628830285\n",
"12 h1 0.05665766362774605\n",
"13 h1 0.056657585826414836\n",
"14 ha 0.07811395098040658\n",
"15 ha 0.07313080623699353\n",
"16 h4 0.0746315390071699\n",
"17 h4 0.07500144690962385\n"
]
}
],
"source": [
"# check generated atoms and charges\n",
"for i, at in enumerate(solute.atom_types):\n",
" print(i, at, solute.charges[i])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"topology = po.make_topology(solute)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"forcefield = ForceField('gaff2.xml', 'tip3p.xml')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"po.generateTemplate(forcefield, topology, solute)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 3.24816249 -0.28545284 -0.04993103]\n",
"[1.92071005 0.42919173 0.14531796]\n",
"[ 0.85091939 -0.4907528 -0.05535284]\n",
"[-0.50768098 -0.15878488 0.04887394]\n",
"[-1.46213108 -1.15677838 -0.17258026]\n",
"[-2.8230515 -0.8610512 -0.07698647]\n",
"[-3.21883223 0.43887336 0.24135335]\n",
"[-2.28344822 1.40174908 0.45475069]\n",
"[-0.95180364 1.13851888 0.36746672]\n",
"[ 3.34237371 -1.11698822 0.68025764]\n",
"[4.08472027 0.42799293 0.10557298]\n",
"[ 3.31027413 -0.69564551 -1.08016704]\n",
"[1.87654003 0.83881484 1.17771217]\n",
"[ 1.84434812 1.26134088 -0.58767369]\n",
"[-1.14759338 -2.16296406 -0.41913518]\n",
"[-3.56354787 -1.63168544 -0.24792887]\n",
"[-4.26989979 0.68286051 0.31891397]\n",
"[-0.2500595 1.94076113 0.54667875]\n"
]
}
],
"source": [
"coords = []\n",
"for i in range(len(solute.xyz)):\n",
" print(solute.xyz[i])\n",
" coords.append(Vec3(solute.xyz[i][0], solute.xyz[i][1], solute.xyz[i][2]))\n",
"positions = coords * angstrom"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Unit({BaseUnit(base_dim=BaseDimension(\"length\"), name=\"angstrom\", symbol=\"A\"): 1.0})"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"angstrom"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[Vec3(x=3.2481624923230448, y=-0.2854528425796883, z=-0.04993103005605603),\n",
" Vec3(x=1.920710052723795, y=0.4291917313962072, z=0.1453179551578719),\n",
" Vec3(x=0.8509193877262783, y=-0.49075280476915223, z=-0.055352838733997316),\n",
" Vec3(x=-0.5076809807429097, y=-0.1587848841626658, z=0.04887394206217463),\n",
" Vec3(x=-1.4621310768713758, y=-1.1567783763036417, z=-0.1725802617719825),\n",
" Vec3(x=-2.8230514983120987, y=-0.8610511955441892, z=-0.0769864701930148),\n",
" Vec3(x=-3.2188322259953295, y=0.438873361620278, z=0.24135335156400248),\n",
" Vec3(x=-2.283448224222744, y=1.4017490836032298, z=0.4547506919670112),\n",
" Vec3(x=-0.9518036417405741, y=1.1385188763514051, z=0.36746672198889957),\n",
" Vec3(x=3.3423737135879765, y=-1.1169882182841946, z=0.6802576399828909),\n",
" Vec3(x=4.084720272244058, y=0.4279929337767665, z=0.10557297623138764),\n",
" Vec3(x=3.3102741303712704, y=-0.6956455148747527, z=-1.0801670357578719),\n",
" Vec3(x=1.8765400251897808, y=0.8388148365810248, z=1.1777121658397807),\n",
" Vec3(x=1.8443481168007743, y=1.2613408837702806, z=-0.5876736912515516),\n",
" Vec3(x=-1.1475933776063765, y=-2.1629640631639386, z=-0.4191351810438239),\n",
" Vec3(x=-3.5635478701375596, y=-1.6316854420307834, z=-0.2479288679473922),\n",
" Vec3(x=-4.269899791194574, y=0.6828605089244459, z=0.31891396548747974),\n",
" Vec3(x=-0.2500595041431521, y=1.940761125689509, z=0.5466787524455193)]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coords"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Quantity(value=[Vec3(x=3.2481624923230448, y=-0.2854528425796883, z=-0.04993103005605603), Vec3(x=1.920710052723795, y=0.4291917313962072, z=0.1453179551578719), Vec3(x=0.8509193877262783, y=-0.49075280476915223, z=-0.055352838733997316), Vec3(x=-0.5076809807429097, y=-0.1587848841626658, z=0.04887394206217463), Vec3(x=-1.4621310768713758, y=-1.1567783763036417, z=-0.1725802617719825), Vec3(x=-2.8230514983120987, y=-0.8610511955441892, z=-0.0769864701930148), Vec3(x=-3.2188322259953295, y=0.438873361620278, z=0.24135335156400248), Vec3(x=-2.283448224222744, y=1.4017490836032298, z=0.4547506919670112), Vec3(x=-0.9518036417405741, y=1.1385188763514051, z=0.36746672198889957), Vec3(x=3.3423737135879765, y=-1.1169882182841946, z=0.6802576399828909), Vec3(x=4.084720272244058, y=0.4279929337767665, z=0.10557297623138764), Vec3(x=3.3102741303712704, y=-0.6956455148747527, z=-1.0801670357578719), Vec3(x=1.8765400251897808, y=0.8388148365810248, z=1.1777121658397807), Vec3(x=1.8443481168007743, y=1.2613408837702806, z=-0.5876736912515516), Vec3(x=-1.1475933776063765, y=-2.1629640631639386, z=-0.4191351810438239), Vec3(x=-3.5635478701375596, y=-1.6316854420307834, z=-0.2479288679473922), Vec3(x=-4.269899791194574, y=0.6828605089244459, z=0.31891396548747974), Vec3(x=-0.2500595041431521, y=1.940761125689509, z=0.5466787524455193)], unit=angstrom)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"positions"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"modeller = Modeller(topology, positions)\n",
"modeller.addSolvent(forcefield, numAdded=50)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"omm_sys = forcefield.createSystem(modeller.topology)\n",
"integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picosecond)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"simulation = Simulation(modeller.topology,\n",
" omm_sys,\n",
" integrator,\n",
" platform=Platform.getPlatformByName('Reference'))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"simulation.context.setPositions(modeller.positions)\n",
"simulation.reporters.append(PDBReporter('output.pdb', 10))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"simulation.step(400)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"simulation.minimizeEnergy(tolerance=2*kilocalorie_per_mole)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"new_z_vals, new_xyz = po.get_atom_positions(simulation.topology, simulation)\n",
"po.write_traj('mol-water_opt_geo_mm.xyz', new_z_vals, new_xyz, comment=\"Optimized Geometry\")"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"new_mol = molecule.Molecule(new_z_vals, new_xyz)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"pk.psi4.core.clean()"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<psi4.core.Molecule at 0x7f8a11dd5a40>"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# The string has molecule and 50 water molecules.\n",
"# it seems not good input format for psi4 (IMHO)\n",
"pk.psi4.geometry(new_mol.to_xyz_string())"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" C 2.107668 -1.041729 0.407400\n",
" C 0.939056 -0.080194 0.152568\n",
" O -0.262383 -0.862826 0.045346\n",
" C -1.487680 -0.296092 -0.247091\n",
" C -2.621257 -1.120669 -0.295859\n",
" C -3.859757 -0.541903 -0.609772\n",
" C -3.918938 0.838723 -0.863944\n",
" N -2.834127 1.639124 -0.823689\n",
" C -1.642206 1.081641 -0.522831\n",
" H 1.941580 -1.627038 1.319255\n",
" H 3.045004 -0.484339 0.522443\n",
" H 2.228201 -1.736929 -0.429190\n",
" H 0.856149 0.637291 0.978373\n",
" H 1.119733 0.470401 -0.777223\n",
" H -2.540967 -2.184617 -0.096930\n",
" H -4.758064 -1.149656 -0.659803\n",
" H -4.864366 1.317780 -1.113085\n",
" H -0.796466 1.766584 -0.499850\n",
" O -1.031373 2.089310 -3.782707\n",
" H -1.582826 2.865109 -3.681382\n",
" H -0.257056 2.274190 -3.251215\n",
" O 0.926622 -0.314881 3.574484\n",
" H 0.800566 -0.107839 4.500483\n",
" H 1.875499 -0.388005 3.471930\n",
" O 0.151567 -1.749689 -2.618348\n",
" H 0.203764 -2.696330 -2.750173\n",
" H 0.130621 -1.640574 -1.667619\n",
" O -0.101947 -4.494288 -2.646100\n",
" H 0.739134 -4.617789 -2.206142\n",
" H -0.713550 -5.040532 -2.152347\n",
" O -3.863970 -2.899203 2.642706\n",
" H -4.804076 -2.904981 2.462707\n",
" H -3.666761 -1.988013 2.859669\n",
" O -2.611683 3.210967 1.802606\n",
" H -3.055211 3.064659 0.967077\n",
" H -2.177436 2.379374 1.992646\n",
" O -0.315058 -2.677895 2.988256\n",
" H 0.068271 -1.800860 2.978206\n",
" H 0.045687 -3.106385 2.212053\n",
" O 0.281579 6.397002 0.251003\n",
" H -0.237125 6.377795 -0.553242\n",
" H 1.187958 6.316468 -0.046020\n",
" O 1.356110 2.496508 -2.418858\n",
" H 1.734788 1.655902 -2.676184\n",
" H 1.566838 2.580028 -1.488885\n",
" O 2.111037 0.050263 -3.347123\n",
" H 1.432488 -0.606210 -3.189497\n",
" H 2.253145 0.026088 -4.293407\n",
" O 4.519567 -2.722550 1.913817\n",
" H 4.239710 -3.544484 2.316728\n",
" H 4.685271 -2.948190 0.998470\n",
" O 0.466181 -0.056117 6.301515\n",
" H 1.012128 -0.052706 7.087747\n",
" H 0.211284 -0.972028 6.190312\n",
" O 3.608239 -0.708743 3.536553\n",
" H 3.711290 0.143497 3.113115\n",
" H 4.042019 -1.324317 2.945676\n",
" O -0.089344 4.312531 1.900630\n",
" H -1.015116 4.069268 1.899369\n",
" H -0.059069 5.144077 1.427507\n",
" O -2.210317 -4.645801 1.218156\n",
" H -1.309772 -4.321563 1.207155\n",
" H -2.709286 -3.953711 1.652059\n",
" O -0.337894 -4.608042 -5.439611\n",
" H -0.067044 -4.732460 -4.530000\n",
" H -1.199759 -4.196014 -5.379162\n",
" O -3.635819 -0.250490 3.152433\n",
" H -2.789525 0.054593 2.825413\n",
" H -3.953101 0.465974 3.702198\n",
" O 1.695327 4.962287 -3.446270\n",
" H 1.921505 4.881289 -4.372830\n",
" H 1.571995 4.060607 -3.149630\n",
" O -2.003345 -0.509471 -3.700604\n",
" H -1.209508 -0.883616 -3.318411\n",
" H -1.866827 0.437014 -3.658649\n",
" O 2.163822 -0.050240 -6.129246\n",
" H 2.508662 0.295378 -6.952572\n",
" H 1.230503 0.160965 -6.152485\n",
" O 3.190692 -2.528520 -5.451578\n",
" H 2.868242 -1.682751 -5.762917\n",
" H 2.798498 -3.168365 -6.045729\n",
" O 2.937270 -2.508834 5.555632\n",
" H 1.981153 -2.463914 5.563039\n",
" H 3.204448 -1.833139 4.932508\n",
" O 4.423246 -5.272847 0.161409\n",
" H 5.225031 -5.781739 0.041403\n",
" H 4.229150 -5.347584 1.095739\n",
" O -4.799739 -0.353073 -3.936228\n",
" H -3.851549 -0.390571 -3.810681\n",
" H -5.037123 -1.230829 -4.235268\n",
" O -1.251422 1.045641 2.692736\n",
" H -1.229924 1.459696 3.555481\n",
" H -0.476809 0.483540 2.676940\n",
" O 1.873907 3.238264 0.246014\n",
" H 1.146472 3.512542 0.804442\n",
" H 2.248343 4.056794 -0.079620\n",
" O 0.149416 -2.689716 5.659342\n",
" H -0.216400 -3.491267 6.033411\n",
" H -0.107102 -2.716245 4.737536\n",
" O 2.282617 -4.346465 -1.259057\n",
" H 3.013838 -4.776625 -0.815763\n",
" H 2.693854 -3.677663 -1.806615\n",
" O -0.647058 5.917814 -2.264758\n",
" H -1.354816 5.415137 -2.668027\n",
" H 0.131179 5.685419 -2.771287\n",
" O -6.292920 1.833494 -3.437638\n",
" H -5.713389 1.085784 -3.583597\n",
" H -7.128663 1.567851 -3.821304\n",
" O 0.538599 -4.199064 0.889129\n",
" H 1.109965 -4.157701 0.122277\n",
" H 0.800141 -5.001753 1.340259\n",
" O 2.680066 5.627104 -0.956323\n",
" H 3.507581 6.096614 -1.061253\n",
" H 2.359990 5.502305 -1.849748\n",
" O 1.252813 -6.461808 2.255360\n",
" H 0.977981 -7.343892 2.005105\n",
" H 1.097940 -6.420867 3.199060\n",
" O 3.677273 -5.243285 2.780312\n",
" H 3.731984 -5.256609 3.735854\n",
" H 2.835716 -5.654175 2.582364\n",
" O 1.771999 -4.191975 -7.079010\n",
" H 1.708854 -4.882996 -7.738354\n",
" H 1.008193 -4.325314 -6.517715\n",
" O 3.069183 1.778150 2.402834\n",
" H 2.822422 2.170998 1.565570\n",
" H 2.525705 2.229116 3.048970\n",
" O -0.814743 2.125243 5.155835\n",
" H -0.032313 2.607881 4.889198\n",
" H -0.482748 1.413055 5.702460\n",
" O -2.634206 -4.959639 4.092430\n",
" H -3.086094 -4.154020 3.841420\n",
" H -2.363316 -5.351512 3.262197\n",
" O 3.549161 -5.108062 5.556253\n",
" H 3.365144 -4.169167 5.585358\n",
" H 4.212466 -5.243700 6.232907\n",
" O -3.541981 2.297996 4.392549\n",
" H -2.716646 2.352280 4.874324\n",
" H -3.426419 2.891926 3.650846\n",
" O -2.430234 -3.257457 -4.076099\n",
" H -1.873328 -3.639603 -3.397828\n",
" H -2.315648 -2.312385 -3.976397\n",
" O -2.527491 4.392251 -3.603878\n",
" H -2.846582 4.789933 -4.413978\n",
" H -3.318692 4.218677 -3.093874\n",
" O -1.947302 -5.838373 -1.209011\n",
" H -2.585726 -6.550943 -1.238855\n",
" H -2.075599 -5.437046 -0.349530\n",
" O -4.952766 3.990544 -2.512616\n",
" H -5.534584 4.304471 -1.820396\n",
" H -5.401268 3.226471 -2.874928\n",
" O 3.447180 -2.431228 -2.761160\n",
" H 3.486191 -2.627724 -3.697162\n",
" H 3.336436 -1.481357 -2.719761\n",
" O -1.034389 -5.109979 6.246609\n",
" H -1.609945 -5.068497 5.482903\n",
" H -1.588318 -5.456689 6.946028\n",
" O -0.430305 0.642032 -6.038374\n",
" H -0.621302 1.245366 -5.320224\n",
" H -1.173667 0.039103 -6.049413\n",
" O 1.359154 3.407952 3.977376\n",
" H 0.826869 3.723062 3.246890\n",
" H 1.781078 4.194296 4.323607\n",
" O 1.102281 -6.264524 4.905671\n",
" H 0.414307 -5.888141 5.454542\n",
" H 1.916747 -5.919208 5.271230\n",
" O -5.056353 -2.963325 -4.731185\n",
" H -5.412191 -3.754739 -5.135258\n",
" H -4.149014 -3.186885 -4.523849\n"
]
}
],
"source": [
"print(new_mol.to_xyz_string())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment