ligand_prep.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2aa3a2d1-e64c-42f6-ba14-3be04c8e14cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/numpy/core/getlimits.py:500: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.\n",
      "  setattr(self, word, getattr(machar, word).flat[0])\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.\n",
      "  return self._float_to_str(self.smallest_subnormal)\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/numpy/core/getlimits.py:500: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.\n",
      "  setattr(self, word, getattr(machar, word).flat[0])\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.\n",
      "  return self._float_to_str(self.smallest_subnormal)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "==============================================================\n",
      "Sending anonymous Sire usage statistics to http://siremol.org.\n",
      "For more information, see http://siremol.org/analytics\n",
      "To disable, set the environment variable 'SIRE_DONT_PHONEHOME' to 1\n",
      "To see the information sent, set the environment variable \n",
      "SIRE_VERBOSE_PHONEHOME equal to 1. To silence this message, set\n",
      "the environment variable SIRE_SILENT_PHONEHOME to 1.\n",
      "==============================================================\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "\n",
      "KeyboardInterrupt\n",
      "\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1037: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1170: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1037: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1170: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1037: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_gromacs.py:1170: UserWarning: The process exited with an error!\n",
      "  _warnings.warn(\"The process exited with an error!\")\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n",
      "/home/li/miniconda3/envs/bss/lib/python3.9/site-packages/BioSimSpace/Process/_process.py:611: UserWarning: The system contains a perturbable molecule .We will assume that you intend to simulate the lambda = 0 state. If you want to simulate the lambda = 1 state, then pass {'is_lambda1' : True} in the 'property_map' argument.\n",
      "  _warnings.warn(\"The system contains a perturbable molecule .\"\n"
     ]
    }
   ],
   "source": [
    "import BioSimSpace as BSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "88b27273-efd5-4eb6-a41d-df729cfe4b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "ejm_46 = BSS.IO.readMolecules(\"org_protein_files/docked_ligands/ejm_46.sdf\")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "944168a3-1f54-4aba-aa8a-4d9775fa5579",
   "metadata": {},
   "outputs": [],
   "source": [
    "ejm_50 = BSS.IO.readMolecules(\"org_protein_files/docked_ligands/ejm_50.sdf\")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a8e33871-cd50-46d9-aee7-63fb12d1b3fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#paramaterising protein - not working. done in tleap\n",
    "AF_protein = BSS.IO.readMolecules([\"org_protein_files/af_dry.rst7\", \"org_protein_files/af_dry.prm7\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "26b6e3f4-4128-4a08-a304-6c8bab807a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#paramaterising ligands\n",
    "ejm_46 = BSS.Parameters.gaff2(ejm_46).getMolecule()\n",
    "ejm_50 = BSS.Parameters.gaff2(ejm_50).getMolecule()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6fa13370-7d06-4725-955a-b320d7bdca0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping\n",
    "mapping = BSS.Align.matchAtoms(ejm_46, ejm_50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0ecdf365-9d08-4640-b818-4b9e91ff7229",
   "metadata": {},
   "outputs": [],
   "source": [
    "# inverse mapping\n",
    "inv_mapping = {v: k for k, v in mapping.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ccd5e532-74ed-4e6a-a188-9a869366910f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Align ejm46 to ejm50 based on the inverse mapping.\n",
    "ejm_50_aligned = BSS.Align.rmsdAlign(ejm_50, ejm_46, inv_mapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d6095e64-fbc0-43c0-a0d1-21a64838082a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge the ejm46 and ejm50 based on the mapping.\n",
    "merged = BSS.Align.merge(ejm_46, ejm_50_aligned, mapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "88aec626-dcd3-4612-b104-16b691945941",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = merged + AF_protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9fa67f88-d615-4edc-9284-fa5258f5abc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check ligand dimensions\n",
    "box_min, box_max = merged.getAxisAlignedBoundingBox()\n",
    "box_size = [y - x for x, y in zip(box_min,box_max)]\n",
    "box_sizes = [x + 20 * BSS.Units.Length.angstrom for x in box_size]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4be131e5-649b-4c5e-8242-cbc3ffdc6563",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating box for free ligands\n",
    "box, angles = BSS.Box.cubic(max(box_sizes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b377d446-8421-4474-b433-3cdb1546092a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.7184e-16 |e|"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solvation for unbounde state\n",
    "solvated = BSS.Solvent.solvate('tip3p', molecule=merged, box=box, angles=angles, is_neutral=True)\n",
    "\n",
    "# Check charge\n",
    "solvated.charge()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0cf01590-d6cd-4946-9731-86af8c9732a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.5915e-07 |e|"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check system dimensions\n",
    "box_min, box_max = system.getAxisAlignedBoundingBox()\n",
    "sys_box_size = [y - x for x, y in zip(box_min,box_max)]\n",
    "sys_box_sizes = [x + 20 * BSS.Units.Length.angstrom for x in sys_box_size]\n",
    "\n",
    "# Creating box for bounded system\n",
    "box, angles = BSS.Box.cubic(max(sys_box_sizes))\n",
    "\n",
    "# Solvation for bound state\n",
    "system_solvated = BSS.Solvent.solvate('tip3p', molecule=system, box=box, angles=angles, is_neutral=True)\n",
    "\n",
    "# Check charge\n",
    "system_solvated.charge()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5d22ea79-d985-4788-a7c3-9ede29535016",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save perturbable system to file\n",
    "BSS.IO.savePerturbableSystem(\"pert_unbonded\", solvated)\n",
    "BSS.IO.savePerturbableSystem(\"pert_bonded\", system_solvated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dfc40c51-616c-4246-bfd7-93071d072b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define minimisation protocol\n",
    "min_protocol = BSS.Protocol.Minimisation(steps=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "32dc1b4a-28b3-42bd-8f93-90370096efa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Minimisation - merged ligands\n",
    "process = BSS.Process.Gromacs(system=solvated, protocol=min_protocol, work_dir=\"Minimisation/free_lig\")\n",
    "process.start()\n",
    "process.wait()\n",
    "minimised_ligands = process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8b541fb8-5bae-49b5-91e9-ce8ab49baf2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Minimisation - system\n",
    "process = BSS.Process.Gromacs(system=system_solvated, protocol=min_protocol, work_dir=\"Minimisation/bound_sys\").start()\n",
    "process.wait()\n",
    "minimised_system = process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9fb034e1-bf7b-486e-ae2b-7553ba636c0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - free ligands - restrained NVT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 5*BSS.Units.Time.picosecond, temperature_start=0*BSS.Units.Temperature.kelvin, temperature_end=298*BSS.Units.Temperature.kelvin, restraint=\"all\")\n",
    "eq_process = BSS.Process.Gromacs(system=minimised_ligands, protocol=eq_protocol, work_dir=\"Equilibration/free_lig/res_NVT\").start()\n",
    "eq_process.wait()\n",
    "eq_resNVT_unbound = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f0990fe-8a10-4013-8008-e52ba83cb014",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibration - free ligands - unrestrained NVT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 50*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin)\n",
    "eq_process = BSS.Process.Gromacs(system=eq_resNVT_unbound, protocol=eq_protocol, work_dir=\"Equilibration/free_lig/unres_NVT\").start()\n",
    "eq_process.wait()\n",
    "eq_unresNVT_unbound = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "96f08f5c-c2b6-418b-80b3-3b1ac0e3aede",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibration - free ligands - restrained NPT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 200*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin, pressure=1*BSS.Units.Pressure.atm, restraint=\"heavy\")\n",
    "eq_process = BSS.Process.Gromacs(system=eq_unresNVT_unbound, protocol=eq_protocol, work_dir=\"Equilibration/free_lig/res_NPT\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8e94d8c3-31c2-44ec-b81d-a78f07c11067",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibration - free ligands - restrained NPT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 200*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin, pressure=1*BSS.Units.Pressure.atm, restraint=\"heavy\")\n",
    "eq_process = BSS.Process.Gromacs(system=eq_unresNVT_unbound, protocol=eq_protocol, work_dir=\"Equilibration/free_lig/res_NPT\").start()\n",
    "eq_process.wait()\n",
    "eq_resNPT_unbound = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f59fee6e-b91f-4c7c-b7b8-0a88923a48be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibration - free ligands - restrained NPT\n",
    "# gmx grompp -f gromacs.mdp -c ../unres_NVT/gromacs.gro -r ../unres_NVT/gromacs.gro -p gromacs.top -t ../unres_NVT/gromacs.cpt -o gromacs.tpr\n",
    "# nohup gmx mdrun -deffnm gromacs -v > nohup_res_NPT.out &"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "7d25cfc7-e439-46ac-ad47-5aca47f4bd8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibration - free ligands - unrestrained NPT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 200*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin, pressure=1*BSS.Units.Pressure.atm)\n",
    "eq_process = BSS.Process.Gromacs(system=eq_resNPT_unbound, protocol=eq_protocol, work_dir=\"Equilibration/free_lig/unres_NPT\").start()\n",
    "eq_process.wait()\n",
    "eq_unbound = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d917d6e-62f3-461d-ac93-2bcf315e93fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "61d68216-812e-4fd4-ac10-e4e71785ae2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - system - restrained NVT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 5*BSS.Units.Time.picosecond, temperature_start=0*BSS.Units.Temperature.kelvin, temperature_end=298*BSS.Units.Temperature.kelvin, restraint=\"all\")\n",
    "eq_process = BSS.Process.Gromacs(system=minimised_system, protocol=eq_protocol, work_dir=\"Equilibration/bound_sys/res_NVT\").start()\n",
    "eq_process.wait()\n",
    "eq_resNVT_system = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "d78af786-30f2-45a4-b68c-3f191c19cfa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - system - backbone restrained NVT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 50*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin, restraint=\"backbone\")\n",
    "eq_process = BSS.Process.Gromacs(system=eq_resNVT_system, protocol=eq_protocol, work_dir=\"Equilibration/bound_sys/b_res_NVT\").start()\n",
    "eq_process.wait()\n",
    "eq_bResNVT_system = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "8ac8c0a1-eed8-4dd4-9c2e-dde3f7796727",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - system - unrestrained NVT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 50*BSS.Units.Time.picosecond, temperature=298*BSS.Units.Temperature.kelvin)\n",
    "eq_process = BSS.Process.Gromacs(system=eq_bResNVT_system, protocol=eq_protocol, work_dir=\"Equilibration/bound_sys/unres_NVT\").start()\n",
    "eq_process.wait()\n",
    "eq_unresNVT_system = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "eccd0567-d8e4-415e-b789-4c6ee995f4f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - system - restrained NPT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 200*BSS.Units.Time.picosecond, pressure=1*BSS.Units.Pressure.atm, temperature=298*BSS.Units.Temperature.kelvin, restraint=\"heavy\")\n",
    "eq_process = BSS.Process.Gromacs(system=eq_unresNVT_system, protocol=eq_protocol, work_dir=\"Equilibration/bound_sys/res_NPT\").start()\n",
    "eq_process.wait()\n",
    "eq_resNPT_system = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "3e23c546-fc79-41bb-9430-db5b47ce242b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibiration - system - unrestrained NPT\n",
    "eq_protocol = BSS.Protocol.Equilibration(runtime= 200*BSS.Units.Time.picosecond, pressure=1*BSS.Units.Pressure.atm, temperature=298*BSS.Units.Temperature.kelvin)\n",
    "eq_process = BSS.Process.Gromacs(system=eq_resNPT_system, protocol=eq_protocol, work_dir=\"Equilibration/bound_sys/unres_NPT\").start()\n",
    "eq_process.wait()\n",
    "eq_system = eq_process.getSystem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75b1224-a9d6-4cc3-9b44-35698afc7528",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "da6bc4f2-74ab-4e2d-8045-adfe45f1b317",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Free energy protocol\n",
    "FE_protocol = BSS.Protocol.FreeEnergy(timestep=2*BSS.Units.Time.femtosecond, runtime=4*BSS.Units.Time.nanosecond, num_lam=9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "235f341e-4bce-46a4-952b-832c692cbd16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Config 2 legs\n",
    "fep_unbound = BSS.FreeEnergy.Relative(eq_unbound, FE_protocol, work_dir=\"ejm46_ejm50/unbound\")\n",
    "fep_bound  = BSS.FreeEnergy.Relative(eq_system, FE_protocol, work_dir=\"ejm46_ejm50/bound\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "7d23c496-9ed9-42ea-8a25-63ed95e40a40",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ProcessRunner already started!\n"
     ]
    }
   ],
   "source": [
    "# Run FEP\n",
    "#fep_unbound.run()\n",
    "#fep_bound.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af055df5-386a-49a7-b059-34610654f4b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analysis\n",
    "pmf_unbound, overlap_unbound = fep_unbound.analyse()\n",
    "pmf_bound,  overlap_bound  = fep_bound.analyse()\n",
    "\n",
    "FE_binding = BSS.FreeEnergy.Relative.difference(pmf_bound, pmf_unbound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d2ebf88-fa58-43dc-bfff-a35e5f485777",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}