{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation of the $\\mu$ opioid receptor with ligand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:45:59,431 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but \"NUMEXPR_MAX_THREADS\" not set, so enforcing safe limit of 8.\n",
      "2024-12-05 12:45:59,431 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.\n",
      "/home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).\n",
      "  from pandas.core.computation.check import NUMEXPR_INSTALLED\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049\n",
      "HTMD Documentation at: https://software.acellera.com/htmd/\n",
      "\n",
      "You are on the latest HTMD version (2.4.2+20.g74f7135a7.dirty).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.home import home\n",
    "#get the files\n",
    "\n",
    "shutil.copytree(home(dataDir=\"mor\"),'/tmp/testmor/pdb')\n",
    "os.chdir('/tmp/testmor')\n",
    "path='./01_prepare/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4dkl.pdb  ff.rtf                 MOL.cif     QM-min.pdb\n",
      "ff.prm    membrane80by80C36.pdb  MOL.frcmod  sod.pdb\n"
     ]
    }
   ],
   "source": [
    "%ls /tmp/testmor/pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:02,233 - moleculekit.molecule - INFO - Removed 5120 atoms. 2262 atoms remaining in the molecule.\n",
      "2024-12-05 12:46:02,268 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 65 and 263.\n",
      "2024-12-05 12:46:02,268 - moleculekit.tools.autosegment - INFO - Created segment P1 between resid 270 and 352.\n",
      "2024-12-05 12:46:02,308 - moleculekit.tools.preparation - WARNING - Both chains and segments are defined in Molecule.chain / Molecule.segid, however they are inconsistent. Protein preparation will use the chain information.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "---- Molecule chain report ----\n",
      "Chain A:\n",
      "    First residue: HOH   717  \n",
      "    Final residue: HOH   717  \n",
      "Chain B:\n",
      "    First residue: MET    65  \n",
      "    Final residue: HOH   735  \n",
      "---- End of chain report ----\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:04,219 - moleculekit.tools.preparation - INFO - Modified residue ASP   114 B to ASH\n",
      "2024-12-05 12:46:04,220 - moleculekit.tools.preparation - INFO - Modified residue CYS   140 B to CYX\n",
      "2024-12-05 12:46:04,220 - moleculekit.tools.preparation - INFO - Modified residue HIS   171 B to HID\n",
      "2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue CYS   217 B to CYX\n",
      "2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue HIS   223 B to HID\n",
      "2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue HIS   297 B to HID\n",
      "2024-12-05 12:46:04,222 - moleculekit.tools.preparation - INFO - Modified residue HIS   319 B to HIE\n",
      "2024-12-05 12:46:04,223 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 1 residues is within 1.0 units of pH 7.4.\n",
      "2024-12-05 12:46:04,224 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    ASP   114 B (pKa= 7.85)\n"
     ]
    }
   ],
   "source": [
    "#Protein 4dkl is taken from opm\n",
    "prot = Molecule('pdb/4dkl.pdb')\n",
    "prot.filter('(protein and noh and chain B) or (water and within 5 of (chain B and protein))')\n",
    "pcenter = np.mean(prot.get('coords','protein'), axis=0)\n",
    "prot = autoSegment(prot, sel='protein') \n",
    "\n",
    "prot = systemPrepare(prot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/vmd/vmd_LINUXAMD64: /lib/x86_64-linux-gnu/libGL.so.1: no version information available (required by /usr/local/lib/vmd/vmd_LINUXAMD64)\n"
     ]
    }
   ],
   "source": [
    "prot.view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:07,287 - moleculekit.molecule - INFO - Removed 327 residues from appended Molecule due to collisions.\n",
      "2024-12-05 12:46:07,493 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb\n",
      "2024-12-05 12:46:07,945 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n",
      "\n",
      "2024-12-05 12:46:16,364 - htmd.builder.solvate - INFO - 9644 water molecules were added to the system.\n",
      "2024-12-05 12:46:18,858 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n",
      "2024-12-05 12:46:18,875 - htmd.builder.builder - INFO - One disulfide bond was added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'B', resid: 78, insertion: '', segid: 'P0'>\n",
      "                   and: UniqueResidueID<resname: 'CYX', chain: 'B', resid: 155, insertion: '', segid: 'P0'>\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:19,910 - htmd.builder.amber - INFO - Starting the build.\n",
      "2024-12-05 12:46:22,046 - htmd.builder.amber - INFO - Finished building.\n",
      "2024-12-05 12:46:23,017 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n",
      "2024-12-05 12:46:28,148 - htmd.builder.builder - WARNING - Found cis peptide bond in 1 frames: [0] in the omega diheral \"Angle of (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) \" with indexes [2528, 2541, 2543, 2553]\n",
      "2024-12-05 12:46:28,452 - htmd.builder.ionize - INFO - Adding 15 anions + 0 cations for neutralizing and 68 ions for the given salt concentration 0.15 M.\n",
      "2024-12-05 12:46:37,889 - htmd.builder.amber - INFO - Starting the build.\n",
      "2024-12-05 12:46:39,968 - htmd.builder.amber - INFO - Finished building.\n",
      "2024-12-05 12:46:40,942 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n",
      "2024-12-05 12:46:45,901 - htmd.builder.builder - WARNING - Found cis peptide bond in 1 frames: [0] in the omega diheral \"Angle of (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) \" with indexes [2528, 2541, 2543, 2553]\n",
      "2024-12-05 12:46:45,936 - py.warnings - WARNING - /home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/Bio/pairwise2.py:278: BiopythonDeprecationWarning: Bio.pairwise2 has been deprecated, and we intend to remove it in a future release of Biopython. As an alternative, please consider using Bio.Align.PairwiseAligner as a replacement, and contact the Biopython developers if you still need the Bio.pairwise2 module.\n",
      "  warnings.warn(\n",
      "\n",
      "2024-12-05 12:46:48,396 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 282 residues: 2-285\n"
     ]
    }
   ],
   "source": [
    "#Add sodium in the receptor\n",
    "sod = Molecule('pdb/sod.pdb')\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "\n",
    "#Use a POPC membrane created with vmd and C36\n",
    "memb = Molecule('pdb/membrane80by80C36.pdb')\n",
    "mcenter = np.mean(memb.get('coords'),axis=0)\n",
    "memb.moveBy(pcenter-mcenter)\n",
    "mol = prot.copy()\n",
    "mol.append(memb, collisions=True)  # Append membrane and remove colliding atoms\n",
    "\n",
    "#Add ligand, previously parametrized using gaamp\n",
    "lig = Molecule('pdb/MOL.cif') \n",
    "lig.set('segid','L')\n",
    "lcenter = np.mean(lig.get('coords'),axis=0)\n",
    "newlcenter = [np.random.uniform(-10, 10), np.random.uniform(-10, 10),  43]\n",
    "lig.rotateBy(uniformRandomRotation(), lcenter)\n",
    "lig.moveBy(newlcenter - lcenter)\n",
    "mol.append(lig) \n",
    "\n",
    "#Add water\n",
    "coo = mol.get('coords','lipids or protein')\n",
    "m = np.min(coo,axis=0) + [0,0,-5]\n",
    "M = np.max(coo,axis=0) + [0,0,20]\n",
    "mol = solvate(mol, minmax=np.vstack((m,M)))\n",
    "\n",
    "#Build\n",
    "topos  = amber.defaultTopo() + ['pdb/MOL.cif']\n",
    "params = amber.defaultParam() + ['pdb/MOL.frcmod']\n",
    "mol = amber.build(mol, topo=topos, param=params, outdir=os.path.join(path,'build'), saltconc=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:48,911 - rdkit - INFO - Enabling RDKit 2024.03.5 jupyter extensions\n",
      "2024-12-05 12:46:48,993 - acemd - INFO - # You are on the latest ACEMD version (4.0.1).\n",
      "2024-12-05 12:46:48,994 - acemd - INFO - Copying ./01_prepare/build/structure.pdb to ./01_prepare/equil/structure.pdb\n",
      "2024-12-05 12:46:48,997 - acemd - INFO - Copying ./01_prepare/build/structure.prmtop to ./01_prepare/equil/structure.prmtop\n",
      "2024-12-05 12:46:50,432 - acemd - WARNING - Found cis peptide bond with dihedral angle -3.32 deg in the omega diheral (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) with indexes [2528 2541 2543 2553]\n"
     ]
    }
   ],
   "source": [
    "from acemd.protocols import setup_equilibration\n",
    "\n",
    "# Use a 10A flat bottom potential on the ligand coordinates to prevent the ligand from diffusing \n",
    "# from its original position during equilibration.\n",
    "# You can refer to https://software.acellera.com/acemd/manual.html#extforces-options \n",
    "# for more information on the restraint options\n",
    "restr = {\n",
    "    \"type\": \"positionalRestraint\",\n",
    "    \"sel\": \"resname MOL and noh\",\n",
    "    \"fbwidth\": [10, 10, 10],\n",
    "    \"setpoints\": [\"5@0ns\"],\n",
    "}\n",
    "\n",
    "setup_equilibration(\n",
    "    os.path.join(path,'build'), \n",
    "    os.path.join(path,'equil'),\n",
    "    run=\"40ns\",\n",
    "    temperature=300,\n",
    "    barostatconstratio=True,\n",
    "    extforces=[restr]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:46:51,491 - jobqueues.util - INFO - Trying to determine all GPU devices\n",
      "2024-12-05 12:46:51,547 - jobqueues.localqueue - INFO - Using GPU devices 0\n",
      "2024-12-05 12:46:51,548 - jobqueues.util - INFO - Trying to determine all GPU devices\n",
      "2024-12-05 12:46:51,585 - jobqueues.localqueue - INFO - Queueing /tmp/testmor/01_prepare/equil\n",
      "2024-12-05 12:46:51,586 - jobqueues.localqueue - INFO - Running /tmp/testmor/01_prepare/equil on device 0\n",
      "2024-12-05 12:48:19,471 - jobqueues.localqueue - INFO - Completed /tmp/testmor/01_prepare/equil\n"
     ]
    }
   ],
   "source": [
    "mdx = LocalGPUQueue()\n",
    "mdx.submit(os.path.join(path, 'equil'))\n",
    "mdx.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:50:16,763 - acemd - INFO - Copying ./01_prepare/equil/output.coor to ./01_prepare/prod/input.coor\n",
      "2024-12-05 12:50:16,766 - acemd - INFO - Copying ./01_prepare/equil/output.xsc to ./01_prepare/prod/input.xsc\n",
      "2024-12-05 12:50:16,767 - acemd - INFO - Copying ./01_prepare/equil/structure.prmtop to ./01_prepare/prod/structure.prmtop\n",
      "2024-12-05 12:50:17,710 - acemd - WARNING - Found cis peptide bond with dihedral angle -10.44 deg in the omega diheral (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) with indexes [2528 2541 2543 2553]\n"
     ]
    }
   ],
   "source": [
    "from acemd.protocols import setup_production\n",
    "\n",
    "# Apply a flat bottom potential to prevent the ligand from changing from the \n",
    "# extra- to the intra-cellular side of the membrane from the periodic image.\n",
    "# Link the box to the center of mass of the membrane but offset the box\n",
    "# center in the positive z direction by 30 Angstrom so that it's situated\n",
    "# above the membrane. The box will have 60A width in the z direction.\n",
    "\n",
    "restr = {\n",
    "    \"type\": \"positionalRestraint\",\n",
    "    \"sel\": \"resname MOL and noh\",\n",
    "    \"fbwidth\": [70, 70, 60],\n",
    "    \"fbcenter\": \"lipids\",\n",
    "    \"fbcenteroffset\": [0, 0, 30],\n",
    "    \"setpoints\": [\"5@0ns\"],\n",
    "}\n",
    "\n",
    "setup_production(\n",
    "    os.path.join(path,'equil'), \n",
    "    os.path.join(path,'prod'),\n",
    "    run=\"50ns\",\n",
    "    temperature=300,\n",
    "    extforces=[restr]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-05 12:55:59,602 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n"
     ]
    }
   ],
   "source": [
    "# Read in the last frame of the equilibration\n",
    "mol = Molecule(os.path.join(path,'equil','structure.prmtop'))\n",
    "mol.read(os.path.join(path,'equil','output.coor'))\n",
    "mol.read(os.path.join(path,'equil','output.xsc'))\n",
    "mol.view('not water')\n",
    "\n",
    "# Visualize the FB box\n",
    "width = np.array([70, 70, 60])\n",
    "fbcenter = mol.getCenter(\"lipids\", com=True)\n",
    "fbcenteroffset = np.array([0, 0, 30])\n",
    "b = VMDBox(np.vstack((fbcentre + fbcenteroffset - width/2, fbcentre + fbcenteroffset + width/2)).T.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}