Water box example for topology PR121

topology-water-box-test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "from foyer import Forcefield\n",
    "import unyt as u\n",
    "\n",
    "import topology\n",
    "from topology import ForceField\n",
    "from topology.external import convert_mbuild\n",
    "from topology.external import convert_parmed\n",
    "from topology.formats.gro import write_gro\n",
    "from topology.formats.top import write_top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write out the SPC/E water mol2 file\n",
    "\n",
    "with open(\"spce.mol2\", \"w\") as f:\n",
    "    f.write(\n",
    "\"\"\"@<TRIPOS>MOLECULE\n",
    "RES\n",
    "3 0 1 0 1\n",
    "SMALL\n",
    "NO_CHARGES\n",
    "@<TRIPOS>CRYSIN\n",
    "    3.0130     3.0130     3.0130    90.0000    90.0000    90.0000  1  1\n",
    "@<TRIPOS>ATOM\n",
    "       1 O            0.0000     0.0000     0.0000 O             1 RES     \n",
    "       2 H           -0.6126    -0.7357     0.0000 H             1 RES     \n",
    "       3 H           -0.5469     0.7762     0.0000 H             1 RES     \n",
    "@<TRIPOS>BOND\n",
    "     1     1     2    1\n",
    "     2     1     3    1\n",
    "@<TRIPOS>SUBSTRUCTURE\n",
    "       1 RES             1 RESIDUE    0 **** ROOT      0\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in single water (SPC/E) structure\n",
    "\n",
    "water = mb.load(\"spce.mol2\")\n",
    "water = water.children[0]\n",
    "water.name = \"water\"\n",
    "water[0].name = \"opls_116\"\n",
    "water[1].name = water[2].name = \"opls_117\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill a box with 100 water molecule\n",
    "\n",
    "water_box = mb.fill_box(\n",
    "    compound=water,     \n",
    "    n_compounds=100,\n",
    "    density=1000,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write out the SPC/E water xml file\n",
    "\n",
    "with open(\"spce.xml\", \"w\") as f:\n",
    "    f.write(\n",
    "\"\"\"<?xml version='1.0' encoding='UTF-8'?>\n",
    "<ForceField name=\"spce\" version=\"0.0.1\">\n",
    "  <FFMetaData>\n",
    "    <Units energy=\"kcal/mol\" mass=\"amu\" charge=\"coulomb\" distance=\"nm\"/>\n",
    "  </FFMetaData>\n",
    "  <AtomTypes expression=\"4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)\">\n",
    "    <ParametersUnitDef parameter=\"epsilon\" unit=\"kJ/mol\"/>\n",
    "    <ParametersUnitDef parameter=\"sigma\" unit=\"nm\"/>\n",
    "    <AtomType name=\"opls_116\" atomclass=\"OW\" element=\"O\" charge=\"-0.8476\" definition=\"O\">\n",
    "      <Parameters>\n",
    "        <Parameter name=\"epsilon\" value=\"0.650194\"/>\n",
    "        <Parameter name=\"sigma\" value=\"0.316557\"/>\n",
    "      </Parameters>\n",
    "    </AtomType>\n",
    "    <AtomType name=\"opls_117\" atomclass=\"HW\" element=\"H\" charge=\"0.4238\" definition=\"H\">\n",
    "      <Parameters>\n",
    "        <Parameter name=\"epsilon\" value=\"0.0\"/>\n",
    "        <Parameter name=\"sigma\" value=\"0.1\"/>\n",
    "      </Parameters>\n",
    "    </AtomType>\n",
    "  </AtomTypes>\n",
    "  <BondTypes expression=\"0.5 * k * (r-r_eq)**2\">\n",
    "    <ParametersUnitDef parameter=\"k\" unit=\"kJ/mol/nm**2\"/>\n",
    "    <ParametersUnitDef parameter=\"r_eq\" unit=\"nm\"/>\n",
    "    <BondType name=\"BondType-Harmonic-1\" type1=\"opls_116\" type2=\"opls_117\">\n",
    "      <Parameters>\n",
    "        <Parameter name=\"k\" value=\"345000.0\"/>\n",
    "        <Parameter name=\"r_eq\" value=\"0.1\"/>\n",
    "      </Parameters>\n",
    "    </BondType>\n",
    "  </BondTypes>\n",
    "  <AngleTypes expression=\"0.5 * k * (theta - theta_eq)**2\">\n",
    "    <ParametersUnitDef parameter=\"k\" unit=\"kJ/mol/degree**2\"/>\n",
    "    <ParametersUnitDef parameter=\"theta_eq\" unit=\"degree\"/>\n",
    "    <AngleType name=\"AngleType-Harmonic-1\" type1=\"opls_116\" type2=\"opls_117\" type3=\"opls_117\">\n",
    "      <Parameters>\n",
    "        <Parameter name=\"k\" value=\"383.0\"/>\n",
    "        <Parameter name=\"theta_eq\" value=\"1.91061193\"/>\n",
    "      </Parameters>\n",
    "    </AngleType>\n",
    "  </AngleTypes>\n",
    "</ForceField>\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in topology forcefield\n",
    "\n",
    "forcefield = ForceField(\"spce.xml\")\n",
    "\n",
    "# Charges aren't being parsed correctly?\n",
    "forcefield.atom_types[\"opls_116\"].charge = -0.8476\n",
    "forcefield.atom_types[\"opls_117\"].charge = 0.4238"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a topology from the mbuild compound\n",
    "\n",
    "top = convert_mbuild.from_mbuild(water_box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign atom types\n",
    "for atom in top.sites:\n",
    "    atom.atom_type = forcefield.atom_types[atom.name]\n",
    "\n",
    "top.update_sites()\n",
    "top.update_atom_types()\n",
    "    \n",
    "# Assign bond types\n",
    "for bond in top.bonds:\n",
    "    bond.bond_type = bond.connection_type = forcefield.bond_types[\"opls_116~opls_117\"]\n",
    "\n",
    "top.update_bonds()\n",
    "top.update_bond_types()\n",
    "\n",
    "# Create angles with correct atom type and add to top\n",
    "for subtop in top.subtops:\n",
    "    angle = topology.core.angle.Angle(\n",
    "        connection_members=[site for site in subtop.sites],\n",
    "        name=\"opls_116~opls_117~opls_117\",\n",
    "        connection_type=forcefield.angle_types[\"opls_116~opls_117~opls_117\"]\n",
    "    )\n",
    "    top.add_connection(angle)\n",
    "\n",
    "top.update_angles()\n",
    "top.update_angle_types()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write out gro and top files\n",
    "\n",
    "write_gro(top, \"system.gro\")\n",
    "write_top(top, \"system.top\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat system.top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write out an MDP file for an NPT simulation\n",
    "\n",
    "with open(\"npt.mdp\", \"w\") as f:\n",
    "    f.write(\n",
    "\"\"\"\n",
    "constraints         = h-bonds\n",
    "integrator          = md\n",
    "nsteps              = 1000000\n",
    "dt                  = 0.001\n",
    "\n",
    "nstxtcout           = 1000\n",
    "nstenergy           = 1000\n",
    "nstlog              = 1000\n",
    "\n",
    "cutoff-scheme       = Verlet\n",
    "ns_type             = grid\n",
    "nstlist             = 10\n",
    "rcoulomb            = 0.9\n",
    "rvdw                = 0.9 \n",
    "\n",
    "coulombtype         = PME\n",
    "\n",
    "gen_vel             = no\n",
    "\n",
    "tcoupl              = nose-hoover\n",
    "tc-grps             = System\n",
    "tau_t               = 1\n",
    "ref_t               = 300\n",
    "\n",
    "pcoupl              = Parrinello-Rahman\n",
    "pcoupltype          = isotropic\n",
    "tau-p               = 10\n",
    "ref-p               = 1\n",
    "compressibility     = 1e-5\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gmx grompp -f npt.mdp -c system.gro -p system.top -o npt.tpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gmx mdrun -v -deffnm npt"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}