{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical Molecular Interaction Potentials tutorial using BioExcel Building Blocks (biobb)\n",
    "\n",
    "***\n",
    "This tutorial aims to illustrate the process of computing **classical molecular interaction potentials** from **protein structures**, step by step, using the **BioExcel Building Blocks library (biobb)**. Examples shown are **Molecular Interaction Potentials (MIPs) grids, protein-protein/ligand interaction potentials, and protein titration**. The particular structures used are the **Lysozyme** protein (PDB code [1AKI](https://www.rcsb.org/structure/1aki), [https://doi.org/10.2210/pdb1AKI/pdb](https://doi.org/10.2210/pdb1AKI/pdb)), the **Epidermal Growth Factor Receptor** kinase domain (PDB code [4HJO](https://www.rcsb.org/structure/4HJO), [https://doi.org/10.2210/pdb4HJO/pdb](https://doi.org/10.2210/pdb4HJO/pdb)) complexed with the **Erlotinib** inhibitor (PDB code [AQ4](https://www.rcsb.org/ligand/AQ4), DrugBank Ligand Code [DB00530](https://go.drugbank.com/drugs/DB00530)),  and a MD simulation of the complex formed by the **SARS-CoV-2 Receptor Binding Domain and the human Angiotensin Converting Enzyme 2** (PDB code [6VW1](https://www.rcsb.org/structure/6vw1), [https://doi.org/10.2210/pdb6VW1/pdb](https://doi.org/10.2210/pdb6VW1/pdb)). \n",
    "\n",
    "The code wrapped is the ***Classical Molecular Interaction Potentials ([CMIP](https://anaconda.org/bioconda/cmip))*** code:\n",
    "\n",
    "**Classical molecular interaction potentials: Improved setup procedure in molecular dynamics simulations of proteins.**\n",
    "*Gelpí, J.L., Kalko, S.G., Barril, X., Cirera, J., de la Cruz, X., Luque, F.J. and Orozco, M. (2001)*\n",
    "*Proteins, 45: 428-437. https://doi.org/10.1002/prot.1159*\n",
    "***\n",
    "\n",
    "## Settings\n",
    "\n",
    "### Biobb modules used\n",
    "\n",
    " - [biobb_io](https://github.com/bioexcel/biobb_io): Tools to fetch biomolecular data from public databases.\n",
    " - [biobb_cmip](https://github.com/bioexcel/biobb_cmip): Tools to compute classical molecular interaction potentials from protein structures.\n",
    " - [biobb_structure_utils](https://github.com/bioexcel/biobb_structure_utils): Tools to modify or extract information from a PDB structure.\n",
    " - [biobb_chemistry](https://github.com/bioexcel/biobb_chemistry): Tools to perform chemoinformatics on molecular structures.\n",
    " - [biobb_amber](https://github.com/bioexcel/biobb_amber): Tools to setup and simulate atomistic MD simulations using AMBER MD package.\n",
    "  \n",
    "### Auxiliary libraries used\n",
    "\n",
    "* [jupyter](https://jupyter.org/): Free software, open standards, and web services for interactive computing across all programming languages.\n",
    "* [nglview](http://nglviewer.org/#nglview): Jupyter/IPython widget to interactively view molecular structures and trajectories in notebooks.\n",
    "* [plotly](https://plot.ly/python/offline/): Python interactive graphing library integrated in Jupyter notebooks.\n",
    "* [simpletraj](https://github.com/arose/simpletraj): Lightweight coordinate-only trajectory reader based on code from GROMACS, MDAnalysis and VMD.\n",
    "\n",
    "### Conda Installation and Launch\n",
    "\n",
    "```console\n",
    "git clone https://github.com/bioexcel/biobb_wf_cmip.git\n",
    "cd biobb_wf_cmip\n",
    "conda env create -f conda_env/environment.yml\n",
    "conda activate biobb_wf_cmip\n",
    "jupyter-notebook biobb_wf_cmip/notebooks/biobb_wf_cmip.ipynb\n",
    "```\n",
    "\n",
    "***\n",
    "## Pipeline steps\n",
    " 1. [Input Parameters](#input)\n",
    " 2. [Fetching PDB structure](#fetch)\n",
    " 3. [CMIP PDB preparation (from PDB structure)](#preparePDB)\n",
    " 4. [Structural water molecules & ions](#titration)\n",
    " 5. [Molecular Interaction Potentials](#mips)\n",
    "    1. [Positive MIP (MIP+)](#mip_pos) \n",
    "    2. [Negative MIP (MIP-)](#mip_neg) \n",
    "    3. [Neutral MIP (MIPn)](#mip_neutral) \n",
    " 6. [Interaction Potential Energies](#interaction)\n",
    "     1. [Protein-Ligand Interaction Energies](#prot-lig) \n",
    "        1. [Fetching PDB structure](#fetch2) \n",
    "        2. [Removing water molecules](#delWater)\n",
    "        3. [Creating ligand topology](#ligandTop)\n",
    "        4. [Generating system topology](#systemTop)\n",
    "        5. [Minimizing the energy of the system](#minSystem)\n",
    "        6. [Preparing the structures for CMIP](#prepProtLigCMIP)\n",
    "        7. [Computing the Protein-Ligand interaction energies](#cmip_protein_ligand)\n",
    "     2. [Protein-Protein Interaction Energies](#prot-prot)\n",
    "        1. [CMIP PDB preparation (from MD)](#preparePDB_MD)\n",
    "        2. [CMIP Boxes](#cmip_boxes)\n",
    "        3. [RDB Interaction Potential Energies](#RBDinteraction)\n",
    "        4. [hACE2 Interaction Potential Energies](#hACE2interaction)\n",
    " 8. [Questions & Comments](#questions)\n",
    " \n",
    "***\n",
    "<img src=\"https://bioexcel.eu/wp-content/uploads/2019/04/Bioexcell_logo_1080px_transp.png\" alt=\"Bioexcel2 logo\"\n",
    "\ttitle=\"Bioexcel2 logo\" width=\"400\" />\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Initializing colab\n",
    "The cell below is used only in case this notebook is executed via **Google Colab**. This process can take a **few minutes** since **miniforge** is **downloaded** and **installed** first, and then the **conda environment** must be created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only executed when using google colab\n",
    "import sys\n",
    "import os\n",
    "if 'google.colab' in sys.modules:\n",
    "  !git clone https://github.com/bioexcel/biobb_wf_cmip.git\n",
    "  !wget https://github.com/conda-forge/miniforge/releases/download/26.1.0-0/Miniforge3-26.1.0-0-Linux-x86_64.sh\n",
    "  !bash Miniforge3-26.1.0-0-Linux-x86_64.sh -b -p /usr/local/miniforge\n",
    "  !/usr/local/miniforge/condabin/mamba create -f biobb_wf_cmip/conda_env/environment.yml -y\n",
    "  %env PATH=/usr/local/bin:/opt/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/tools/node/bin:/tools/google-cloud-sdk/bin:/usr/local/miniforge/envs/biobb_wf_cmip/bin\n",
    "  _ = (sys.path.append(\"/usr/local/miniforge/envs/biobb_wf_cmip/lib/python3.12/site-packages/\"))\n",
    "  from google.colab import output\n",
    "  output.enable_custom_widget_manager()\n",
    "  os.chdir('/content/biobb_wf_cmip/biobb_wf_cmip/notebooks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"input\"></a>\n",
    "## Input parameters\n",
    "**Input parameters** needed:\n",
    " - **pdbCode**: PDB code of the protein structure (e.g. 1AKI, [https://doi.org/10.2210/pdb1AKI/pdb](https://doi.org/10.2210/pdb1AKI/pdb))\n",
    " - **complexCode**: PDB code of the protein-ligand complex (e.g. 4HJO [*EGFR + Erlotinib*], [https://doi.org/10.2210/pdb4HJO/pdb](https://doi.org/10.2210/pdb4HJO/pdb)) \n",
    " - **ligandCode**: PDB code of the ligand (e.g. AQ4 [*Erlotinib*], DrugBank Ligand Code [DB00530](https://go.drugbank.com/drugs/DB00530))\n",
    " - **mol_charge**: Charge of the small molecule\n",
    " - **MDCode**: Code of the Molecular Dynamics trajectory (e.g. RBD-hACE2)\n",
    "     - **inputPDB_MD**: MD reference structure (PDB format)\n",
    "     - **inputTOP_MD**: MD topology (Amber Parmtop7 format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview\n",
    "import ipywidgets\n",
    "import sys\n",
    "\n",
    "pdbCode = \"1aki\"       # Structure of the orthorhombic form of Hen Egg-white Lysozyme\n",
    "complexCode = \"4hjo\"   # Crystal structure of the inactive EGFR tyrosine kinase domain with erlotinib\n",
    "ligandCode = \"AQ4\"     # Erlotinib epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor\n",
    "mol_charge = 0         # Molecular charge for Erlotinib\n",
    "\n",
    "MDCode = \"RBD-hACE2-ZN\" \n",
    "inputPDB_MD = \"Files/\" + MDCode + \".pdb\" # Structure of the SARS-CoV-2 RBD-hACE2 complex\n",
    "inputTOP_MD = \"Files/\" + MDCode + \".top\" # MD Topology of the SARS-CoV-2 RBD-hACE2 complex\n",
    "\n",
    "vdw_params_path='/usr/local/miniforge/envs/biobb_wf_cmip/share/cmip/dat/vdwprm' if 'google.colab' in sys.modules else None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fetch\"></a>\n",
    "***\n",
    "## Fetching PDB structure\n",
    "Downloading **PDB structure** with the **protein molecule** from the RCSB PDB database.<br>\n",
    "Alternatively, a **PDB file** can be used as starting structure. <br>\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [Pdb](https://biobb-io.readthedocs.io/en/latest/api.html#module-api.pdb) from **biobb_io.api.pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downloading desired PDB file \n",
    "# Import module\n",
    "from biobb_io.api.pdb import pdb\n",
    "\n",
    "# Create properties dict and inputs/outputs\n",
    "downloaded_pdb = pdbCode+'.pdb'\n",
    "prop = {\n",
    "    'pdb_code': pdbCode,\n",
    "    'api_id' : 'pdbe'\n",
    "}\n",
    "\n",
    "#Create and launch bb\n",
    "pdb(output_pdb_path=downloaded_pdb,\n",
    "    properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the downloaded/given **PDB structure** using **NGL**: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(downloaded_pdb)\n",
    "view.add_representation(repr_type='ball+stick', selection='all')\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"preparePDB\"></a>\n",
    "***\n",
    "## CMIP PDB Preparation (from PDB structure)\n",
    "**CMIP** tool needs additional information (e.g. charges, elements) to be included in the **structure PDB file** to properly run. A specific **BioBB building block** (prepare_pdb) is used in the next cell to prepare the **input PDB file**, adding this extra information. **Charges and elements** are taken from an internal **CMIP library** based on the **AMBER force fields**. \n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_prepare_pdb](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#cmip-prepare-pdb-module) from **biobb_cmip.cmip.cmip_prepare_pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_prepare_pdb import cmip_prepare_pdb\n",
    "\n",
    "cmipPDB = pdbCode + \".cmip.pdb\"\n",
    "\n",
    "cmip_prepare_pdb(input_pdb_path=downloaded_pdb,\n",
    "            output_cmip_pdb_path=cmipPDB\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"titration\"></a>\n",
    "***\n",
    "## Structural water molecules & ions\n",
    "One of the many steps involved in the **MD structure setup process** is the addition of **solvent and counterions** (when working with explicit solvent). **Solvent molecules** and **counterions** are usually integrated on the **structure surface** in two steps:\n",
    "- **Structural waters/ions**: A **first shell** of **water molecules and ions** is commonly added in the **energetically most favorable positions** on the surface of the structure. It is a computationally expensive process and is usually reduced to just tens of **water molecules** and **ions** (depending on the structure size).\n",
    "- **Solvent box/ionic concentration**: A box of **solvent molecules** is created surrounding the original structure, and an additional number of **ions** are added until reaching a desired **ionic concentration**.\n",
    "\n",
    "Whereas the second step is integrated in all the **MD packages**, the first one is rarely available. **CMIP** and the **biobb_titration building block** is helping in this task.   \n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_titration](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#cmip-titration-module) from **biobb_cmip.cmip.cmip_titration**\n",
    " - [cat_pdb](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.cat_pdb) from **biobb_structure_utils.utils.cat_pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"run_titration\"></a>\n",
    "### Computing structural water molecules & ions positions\n",
    "Computing the positions of **20 structural water molecules**, **5 positive** and **5 negative ions** in the most **energetically favourable** regions of the **structure surface**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_titration import cmip_titration\n",
    "\n",
    "wat_ions_pdb = pdbCode + \".wat_ions.pdb\"\n",
    "wat_ions_log = pdbCode + \".wat_ions.log\"\n",
    "\n",
    "prop = { \n",
    "#    'neutral' : True, # Can be also used to neutralize the system\n",
    "    'num_positive_ions' : 5,\n",
    "    'num_negative_ions' : 5,    \n",
    "    'num_wats' : 20\n",
    "}\n",
    "\n",
    "cmip_titration(input_pdb_path=cmipPDB,\n",
    "          input_vdw_params_path=vdw_params_path,\n",
    "          output_pdb_path=wat_ions_pdb,\n",
    "          output_log_path=wat_ions_log,\n",
    "          properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"catPDB_tit\"></a>\n",
    "### Adding structural water molecules & ions\n",
    "Adding the 20 + 10 computed **structural water molecules** and **ions** to the original **PDB file**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_structure_utils.utils.cat_pdb import cat_pdb\n",
    "\n",
    "titPDB = pdbCode + \".tit.pdb\"\n",
    "\n",
    "cat_pdb(input_structure1=cmipPDB,\n",
    "       input_structure2=wat_ions_pdb,\n",
    "       output_structure_path=titPDB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visTIT\"></a>\n",
    "### Visualizing structural water molecules & ions\n",
    "Visualizing the recently added **structural water molecules and ions**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(titPDB)\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view.add_representation(repr_type='spacefill', selection='water')\n",
    "view.add_representation(repr_type='spacefill', selection='.Na', color='element')\n",
    "view.add_representation(repr_type='spacefill', selection='.Cl', color='element')\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"mips\"></a>\n",
    "***\n",
    "## Molecular Interaction Potentials (MIPs)\n",
    "**Molecular interaction potentials (MIP)** are field properties arising from the interaction of a **probe** (e.g., methyl, proton or water) with a molecule. These are calculated in the surface of the molecule, with a grid defined around the structure.\n",
    "\n",
    "**MIPs** are one of the most important molecular properties in the relationship between **molecular and binding data** (e.g. *3D Quantitative Structure-Activity Relationships, 3D-QSAR*), and is extensively applied in **drug discovery** processes.  \n",
    "\n",
    "In this example, three different **MIPs** are used, with a **Water Oxygen atom** as a probe:\n",
    " - **Positive** MIP - highlighting the protein regions with **higher affinity** to **negatively charged groups**.\n",
    " - **Negative** MIP - highlighting the protein regions with **higher affinity** to **positively charged groups**.\n",
    " - **Neutral** MIP - highlighting the protein regions with **lower affinity** to **electrocharged groups**.\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.cmip) from **biobb_cmip.cmip.cmip**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"mip_pos\"></a>\n",
    "### Positive MIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "mip_pos_log = pdbCode + \".mip_pos.log\"\n",
    "mip_pos_cube = pdbCode + \".mip_pos.cube\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'mip_pos'\n",
    "}\n",
    "\n",
    "cmip_run( input_pdb_path=cmipPDB,\n",
    "          input_vdw_params_path=vdw_params_path,\n",
    "          output_log_path=mip_pos_log,\n",
    "          output_cube_path=mip_pos_cube,\n",
    "          properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visMIP1\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the **positive MIP grid**, with protein regions with **higher affinity** to **negatively charged groups** highlighted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(mip_pos_cube)\n",
    "view.add_component(nglview.FileStructure(cmipPDB))\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view.add_surface(isolevelType=\"value\", isolevel=-5, color=\"blue\")\n",
    "view.component_1.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"mip_neg\"></a>\n",
    "### Negative MIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "mip_neg_log = pdbCode + \".mip_neg.log\"\n",
    "mip_neg_cube = pdbCode + \".mip_neg.cube\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'mip_neg'\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_log_path=mip_neg_log,\n",
    "         output_cube_path=mip_neg_cube,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visMIP2\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the **negative MIP grid**, with protein regions with **higher affinity** to **positively charged groups** highlighted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(mip_neg_cube)\n",
    "view.add_component(nglview.FileStructure(cmipPDB))\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view.add_surface(isolevelType=\"value\", isolevel=-10, color=\"red\")\n",
    "view.component_1.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"mip_neutral\"></a>\n",
    "### Neutral MIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "mip_neutral_log = pdbCode + \".mip_neutral.log\"\n",
    "mip_neutral_cube = pdbCode + \".mip_neutral.cube\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'mip_neu'\n",
    "}\n",
    "\n",
    "cmip_run( input_pdb_path=cmipPDB,\n",
    "          input_vdw_params_path=vdw_params_path,\n",
    "          output_log_path=mip_neutral_log,\n",
    "          output_cube_path=mip_neutral_cube,\n",
    "          properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visMIP3\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the **neutral MIP grid**, with protein regions with **lower affinity** to **electrocharged groups** highlighted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(mip_neutral_cube)\n",
    "view.add_component(nglview.FileStructure(cmipPDB))\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view.add_surface(isolevelType=\"value\", isolevel=-1, color=\"grey\")\n",
    "view.component_1.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visMIP4\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing all **MIP grids**, for comparison purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Show different structures generated (for comparison)\n",
    "view1 = nglview.show_structure_file(cmipPDB)\n",
    "view1.add_component(nglview.FileStructure(mip_pos_cube))\n",
    "view1.component_0.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view1.component_1.add_surface(isolevelType=\"value\", isolevel=-5, color=\"blue\")\n",
    "view1.component_0.center()\n",
    "view1._remote_call('setSize', target='Widget', args=['350px','400px'])\n",
    "view1.camera='orthographic'\n",
    "view1\n",
    "view2 = nglview.show_structure_file(cmipPDB)\n",
    "view2.add_component(nglview.FileStructure(mip_neg_cube))\n",
    "view2.component_0.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view2.component_1.add_surface(isolevelType=\"value\", isolevel=-10, color=\"red\")\n",
    "view2.component_0.center()\n",
    "view2._remote_call('setSize', target='Widget', args=['350px','400px'])\n",
    "view2.camera='orthographic'\n",
    "view2\n",
    "view3 = nglview.show_structure_file(cmipPDB)\n",
    "view3.add_component(nglview.FileStructure(mip_neutral_cube))\n",
    "view3.component_0.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view3.component_1.add_surface(isolevelType=\"value\", isolevel=-1, color=\"grey\")\n",
    "view3.component_0.center()\n",
    "view3._remote_call('setSize', target='Widget', args=['350px','400px'])\n",
    "view3.camera='orthographic'\n",
    "view3\n",
    "ipywidgets.HBox([view1, view2, view3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"interaction\"></a>\n",
    "***\n",
    "## Interaction Potential Energies \n",
    "\n",
    "Closely related to the previous study of **Molecular Interaction Potentials**, the **Interaction Potential Energies** calculation computes the contributions to the **total energy** of the system from the different **interactions between the subunits of the molecule** considered. These **interaction energies** usually depend on the relation between the **charge** and **positions** of the units studied (e.g. *electrostatic, van der Waals*) and the **solvation energy** (energy released when a compound is dissolved in a solvent).\n",
    "\n",
    "**Interaction Potential Energies** give useful insights in the **macromolecular interaction** process, with the possibility to identify **key residues** involved in the interaction, and thus being another key component of the **drug discovery** process.  \n",
    "\n",
    "To illustrate the calculation of the **interaction potentials** between two subunits of a **structure complex** (e.g. protein-protein, protein-ligand), two different examples are used:\n",
    "\n",
    "- **Epidermal Growth Factor Receptor** kinase domain with **Erlotinib** inhibitor.\n",
    "- **SARS-CoV-2 Receptor Binding Domain** with **human Angiotensin Converting Enzyme 2**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"prot-lig\"></a>\n",
    "***\n",
    "## Protein-Ligand Interaction Energies\n",
    "\n",
    "This example illustrates the steps needed to compute the **protein-ligand interaction energies** from a **crystal structure** of the **complex** taken directly from the **PDB data bank**. The particular example used is the **Epidermal Growth Factor Receptor** kinase domain (PDB code [4HJO](https://www.rcsb.org/structure/4HJO, [https://doi.org/10.2210/pdb4HJO/pdb](https://doi.org/10.2210/pdb4HJO/pdb)) with **Erlotinib** inhibitor (PDB code [AQ4](https://www.rcsb.org/ligand/AQ4), DrugBank Ligand Code [DB00530](https://go.drugbank.com/drugs/DB00530)). \n",
    "\n",
    "To properly compute the **interaction energies**, the **protein-ligand** complex needs to be pre-processed, adding the (typically) missing **hydrogen atoms**, finding out the **atomic charges** and **elements**, and **relaxing** the structure to avoid **steric issues** due to **crystallographic packing**. All these steps are performed in the following cells using **BioBB building blocks**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fetch2\"></a>\n",
    "***\n",
    "## Fetching PDB structure\n",
    "Downloading **PDB structure** with the **protein molecule** from the RCSB PDB database.<br>\n",
    "Alternatively, a **PDB file** can be used as starting structure. <br>\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [Pdb](https://biobb-io.readthedocs.io/en/latest/api.html#module-api.pdb) from **biobb_io.api.pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downloading desired PDB file \n",
    "# Import module\n",
    "from biobb_io.api.pdb import pdb\n",
    "\n",
    "# Create properties dict and inputs/outputs\n",
    "downloaded_pdb = complexCode+'.pdb'\n",
    "prop = {\n",
    "    'pdb_code': complexCode,\n",
    "    'api_id' : 'pdbe',\n",
    "    'filter': ['ATOM', 'HETATM'],\n",
    "}\n",
    "\n",
    "#Create and launch bb\n",
    "pdb(output_pdb_path=downloaded_pdb,\n",
    "    properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the downloaded/given **PDB structure** using **NGL**: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(downloaded_pdb)\n",
    "view.add_representation(repr_type='ball+stick', selection='hetero')\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"delWater\"></a>\n",
    "***\n",
    "## Removing Water Molecules\n",
    "Removing the **water molecules** from the **downloaded/given structure**.<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [remove_pdb_water](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.remove_pdb_water) from **biobb_structure_utils.utils**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Import module\n",
    "from biobb_structure_utils.utils.remove_pdb_water import remove_pdb_water\n",
    "\n",
    "# Removing Waters:\n",
    "\n",
    "# Create properties dict and inputs/outputs\n",
    "nohoh_pdb = complexCode+'.noHOH.pdb'\n",
    "\n",
    "#Create and launch bb\n",
    "remove_pdb_water(\n",
    "    input_pdb_path=downloaded_pdb,\n",
    "    output_pdb_path=nohoh_pdb\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the **structure** without **water molecules** using **NGL**: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(nohoh_pdb)\n",
    "view.add_representation(repr_type='ball+stick', selection='hetero')\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"ligandTop\"></a>\n",
    "***\n",
    "## Creating Ligand Topology\n",
    "Obtaining parameters for the **small molecule** to be used in the **potential energy** calculations.<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [extract_heteroatoms](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.extract_heteroatoms) from **biobb_structure_utils.utils**\n",
    " - [reduce_add_hydrogens](https://biobb-chemistry.readthedocs.io/en/latest/ambertools.html#module-ambertools.reduce_add_hydrogens) from **biobb_chemistry.ambertools**\n",
    " - [acpype_params_ac](https://biobb-chemistry.readthedocs.io/en/latest/acpype.html#module-acpype.acpype_params_ac) from **biobb_chemistry.acpype**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Extracting Ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Ligand system topology, STEP 1\n",
    "# Extracting Ligand \n",
    "# Import module\n",
    "from biobb_structure_utils.utils.extract_heteroatoms import extract_heteroatoms\n",
    "\n",
    "# Create properties dict and inputs/outputs\n",
    "ligandFile = ligandCode+'.pdb'\n",
    "\n",
    "prop = {\n",
    "     'heteroatoms' : [{\"name\": ligandCode}]\n",
    "}\n",
    "\n",
    "extract_heteroatoms(\n",
    "    input_structure_path=nohoh_pdb,\n",
    "    output_heteroatom_path=ligandFile,\n",
    "    properties=prop\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Add hydrogen atoms to the ligand\n",
    "\n",
    "**NOTE**: Check the protonation state of the small molecule after this step. Tools to automatically add hydrogen atoms may not give the most appropriate protonation state.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Ligand system topology, STEP 2\n",
    "# Reduce_add_hydrogens: add Hydrogen atoms to a small molecule (using Reduce tool from Ambertools package)\n",
    "# Import module\n",
    "from biobb_chemistry.ambertools.reduce_add_hydrogens import reduce_add_hydrogens\n",
    "\n",
    "# Create prop dict and inputs/outputs\n",
    "output_reduce_h = ligandCode+'.reduce.H.pdb' \n",
    "\n",
    "# Create and launch bb\n",
    "reduce_add_hydrogens(\n",
    "    input_path=ligandFile,\n",
    "    output_path=output_reduce_h\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Generating ligand topology "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create Ligand system topology, STEP 3\n",
    "# Acpype_params_gmx: Generation of topologies for AMBER with ACPype\n",
    "# Import module\n",
    "from biobb_chemistry.acpype.acpype_params_ac import acpype_params_ac\n",
    "\n",
    "# Create prop dict and inputs/outputs\n",
    "output_acpype_inpcrd = ligandCode+'params.inpcrd'\n",
    "output_acpype_frcmod = ligandCode+'params.frcmod'\n",
    "output_acpype_lib = ligandCode+'params.lib'\n",
    "output_acpype_prmtop = ligandCode+'params.prmtop'\n",
    "output_acpype = ligandCode+'params'\n",
    "\n",
    "prop = {\n",
    "    'basename' : output_acpype,\n",
    "    'charge' : mol_charge\n",
    "}\n",
    "\n",
    "# Create and launch bb\n",
    "acpype_params_ac(input_path=output_reduce_h, \n",
    "                output_path_inpcrd=output_acpype_inpcrd,\n",
    "                output_path_frcmod=output_acpype_frcmod,\n",
    "                output_path_lib=output_acpype_lib,\n",
    "                output_path_prmtop=output_acpype_prmtop,\n",
    "                 properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing 3D structures\n",
    "Visualizing the **ligand** structures:\n",
    "\n",
    "**Original ligand** PDB structure (left)\n",
    "**Ligand** PDB structure with **hydrogen atoms** (right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show different structures generated (for comparison)\n",
    "view1 = nglview.show_structure_file(ligandFile)\n",
    "view1.add_representation(repr_type='ball+stick')\n",
    "view1._remote_call('setSize', target='Widget', args=['400px','400px'])\n",
    "view1.camera='orthographic'\n",
    "view1\n",
    "view2 = nglview.show_structure_file(output_reduce_h)\n",
    "view2.add_representation(repr_type='ball+stick')\n",
    "view2._remote_call('setSize', target='Widget', args=['400px','400px'])\n",
    "view2.camera='orthographic'\n",
    "view2\n",
    "ipywidgets.HBox([view1, view2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"systemTop\"></a>\n",
    "***\n",
    "## Generating system topology\n",
    "Generating the **system topology** for the **protein-ligand complex** using **tleap** from the **AmberTools package**. The **topology** will then be used to extract the needed **parameters** (e.g. atom charges) for the **potential energy** calculation.<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [leap_gen_top](https://biobb-amber.readthedocs.io/en/latest/leap.html#module-leap.leap_gen_top) from **biobb_amber.leap**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment ALWAYS when running in Google Colab:\n",
    "# import os\n",
    "# os.environ['AMBERHOME'] = \"/usr/local/miniforge/envs/biobb_wf_cmip\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment in case of experiencing issues with undefined AMBERHOME variable in the cell below (when running in Jupyter Notebook / Lab):\n",
    "# import os\n",
    "# os.environ['AMBERHOME'] = \"/path/to/anaconda3/envs/biobb_wf_cmip\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import module\n",
    "from biobb_amber.leap.leap_gen_top import leap_gen_top\n",
    "\n",
    "# Create prop dict and inputs/outputs\n",
    "output_pdb_path = 'structure.leap.pdb'\n",
    "output_top_path = 'structure.leap.prmtop'\n",
    "output_crd_path = 'structure.leap.crd'\n",
    "\n",
    "prop = {\n",
    "    \"forcefield\" : [\"protein.ff14SB\",\"gaff2\"]\n",
    "}\n",
    "\n",
    "# Create and launch bb\n",
    "leap_gen_top(input_pdb_path=nohoh_pdb,\n",
    "           input_lib_path=output_acpype_lib,\n",
    "           input_frcmod_path=output_acpype_frcmod,\n",
    "           output_pdb_path=output_pdb_path,\n",
    "           output_top_path=output_top_path,\n",
    "           output_crd_path=output_crd_path,\n",
    "           properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"minSystem\"></a>\n",
    "***\n",
    "## Minimizing the energy of the system\n",
    "Running a short **energetic minimization** to relax the system and remove **possible clashes** introduced by **crystallographic** constraints.<br>\n",
    "                                                                                           \n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [sander_mdrun](https://biobb-amber.readthedocs.io/en/latest/sander.html#module-sander.sander_mdrun) from **biobb_amber.sander**\n",
    " - [amber_to_pdb](https://biobb-amber.readthedocs.io/en/latest/ambpdb.html#module-ambpdb.amber_to_pdb) from **biobb_amber.ambpdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_amber.sander.sander_mdrun import sander_mdrun\n",
    "\n",
    "trj_amber = \"structure.amber.crd\"\n",
    "rst_amber = \"structure.amber.rst\"\n",
    "log_amber = \"structure.amber.log\"\n",
    "\n",
    "prop = {\n",
    "    'simulation_type' : 'minimization',\n",
    "    'mdin' : {\n",
    "        'ntb' : 0,        # Periodic Boundary. No periodicity is applied and PME (Particle Mesh Ewald) is off.\n",
    "        'cut' : 12,       # Nonbonded cutoff, in Angstroms.\n",
    "        'maxcyc' : 500,   # Maximum number of cycles of minimization.\n",
    "        'ncyc' : 50,      # Minimization will be switched from steepest descent to conjugate gradient after ncyc cycles.\n",
    "        'ntwx' : 50       # Trajectory output frequency\n",
    "    }\n",
    "}\n",
    "\n",
    "sander_mdrun(\n",
    "     input_top_path=output_top_path,\n",
    "     input_crd_path=output_crd_path,\n",
    "     output_traj_path=trj_amber,\n",
    "     output_rst_path=rst_amber,\n",
    "     output_log_path=log_amber,\n",
    "     properties=prop\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting a PDB file from the minimized system "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_amber.ambpdb.amber_to_pdb import amber_to_pdb\n",
    "\n",
    "pdb_amber_min = \"structure.amber-min.pdb\"\n",
    "\n",
    "amber_to_pdb(\n",
    "      input_top_path=output_top_path,\n",
    "      input_crd_path=rst_amber,\n",
    "      output_pdb_path=pdb_amber_min\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"prepProtLigCMIP\"></a>\n",
    "***\n",
    "## Preparing the structures for CMIP\n",
    "Preparing the structures for the **CMIP calculations**, using the energetically relaxed **PDB structure** coordinates and the information from the generated **topology**. Files needed are:\n",
    "\n",
    "* **Protein-ligand Complex**: Protein-Ligand complex, with information needed by **CMIP** (charges, elements)\n",
    "* **Protein isolated**: Structure of the protein in **apo state** (without the ligand)\n",
    "* **Protein-ligand Complex** with protein structure in **dielectric mode** (charges neutralized); Atoms in **dielectric mode** are occupying a **volume** but charges are replaced by a **continuum dielectric**, important for the **interaction energy** calculation. These atoms are identified by **CMIP** with an \"X\" letter just before the element symbol in the input PDB file (automatically inserted by the *ignore_residues* building block)\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_prepare_structure](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.prepare_structure) from **biobb_cmip.cmip.cmip_prepare_structure**\n",
    " - [remove_ligand](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.remove_ligand) from **biobb_structure_utils.utils.remove_ligand**\n",
    " - [cmip_ignore_residues](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.ignore_residues) from **biobb_cmip.cmip.cmip_ignore_residues**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CMIP PDB structure Preparation 1\n",
    "\n",
    "Creating the **Protein-Ligand Complex** PDB file from **PDB structure** and **AMBER topology**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_prepare_structure import cmip_prepare_structure\n",
    "\n",
    "cmipPDB = complexCode + \".cmip.pdb\"\n",
    "\n",
    "cmip_prepare_structure(input_pdb_path=pdb_amber_min,\n",
    "            input_topology_path=output_top_path,            \n",
    "            output_cmip_pdb_path=cmipPDB\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CMIP PDB structure Preparation 2\n",
    "\n",
    "Creating the **Protein isolated** PDB file from the previously generated **CMIP PDB structure**, removing the **ligand**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import module\n",
    "from biobb_structure_utils.utils.remove_ligand import remove_ligand\n",
    "\n",
    "# Create properties dict and inputs/outputs\n",
    "cmipPDB_EGFR_Prot = complexCode+'.noLIG.pdb'\n",
    "\n",
    "prop = {\n",
    "    'ligand' : ligandCode\n",
    "}\n",
    "\n",
    "#Create and launch bb\n",
    "remove_ligand(input_structure_path=cmipPDB,\n",
    "    output_structure_path=cmipPDB_EGFR_Prot,\n",
    "    properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CMIP PDB structure Preparation 3\n",
    "\n",
    "Creating the **Protein-ligand Complex** with protein structure in **dielectric mode**. For this, all residues belonging to the **protein structure** (in this case, residues with id from 1 to 279) are marked with an 'X' symbol in the PDB. This is done by the *ignore_residues* building block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_ignore_residues import cmip_ignore_residues\n",
    "\n",
    "cmipPDB_EGFR_Prot_ignored = complexCode + \".Prot_ignored.cmip.pdb\"\n",
    "\n",
    "protein_residues = list(range(1,279))\n",
    "    \n",
    "prop = {\n",
    "    'residue_list': protein_residues # WARNING: Change residue number accordingly\n",
    "}\n",
    "\n",
    "cmip_ignore_residues(input_cmip_pdb_path = cmipPDB,\n",
    "                     output_cmip_pdb_path = cmipPDB_EGFR_Prot_ignored,\n",
    "                     properties = prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_ligand\"></a>\n",
    "***\n",
    "## Computing the Protein-Ligand interaction energies\n",
    "\n",
    "Computing the **Protein-Ligand interaction energies** using the complex with the protein in **dielectric mode** as a **CMIP host structure** and the **protein structure** as a **CMIP ligand** (probe). \n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.cmip) from **biobb_cmip.cmip.cmip**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "EGFR_energies_log = complexCode + \".EGFR.energies.log\"\n",
    "EGFR_byat_out = complexCode + \".EGFR.ener.byat.out\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'pb_interaction_energy'\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_EGFR_Prot_ignored,\n",
    "         input_probe_pdb_path=cmipPDB_EGFR_Prot,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_log_path=EGFR_energies_log,\n",
    "         output_byat_path=EGFR_byat_out,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_ligand_plot\"></a>\n",
    "***\n",
    "### Plotting the Protein-Ligand interaction energies by atom\n",
    "\n",
    "Plotting the **Protein-Ligand interaction energies** with the **protein atom id** in the x-axis and the **potential energies** in KCal/mol in the y-axis. An auxiliary function included in the **biobb_cmip** module is used to extract the energies by atom (*get_energies_byat*).\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "from biobb_cmip.utils.representation import get_energies_byat\n",
    "\n",
    "atom_list, energy_dict = get_energies_byat(EGFR_byat_out, cutoff=50)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=atom_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Atom Number\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_ligand_plot2\"></a>\n",
    "***\n",
    "### Plotting the Protein-Ligand interaction energies by residue\n",
    "\n",
    "Plotting the **Protein-Ligand interaction energies** with the **protein residue sequence** in the x-axis and the **potential energies** in KCal/mol in the y-axis. An auxiliary function included in the **biobb_cmip** module is used to extract the energies by residues (*get_energies_byres*).\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "import re\n",
    "from biobb_cmip.utils.representation import get_energies_byres\n",
    "\n",
    "res_list, energy_dict = get_energies_byres(EGFR_byat_out, cutoff=55)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=res_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Residue ID\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()\n",
    "\n",
    "energy_cutoff = -1.5\n",
    "identified_residues = [res_list[energy_dict['ES&VDW'].index(item)] \n",
    "                       for item in energy_dict['ES&VDW'] if item < energy_cutoff]\n",
    "identified_residues_ngl = str([int(re.sub(r'[A-Z]+ ', '', item)) for item in identified_residues]).replace(',','')\n",
    "\n",
    "print (\"Residues involved in the interaction (energy < \" + str(energy_cutoff) + \" KCal/mol):\")\n",
    "print (identified_residues_ngl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_ligand_ngl\"></a>\n",
    "***\n",
    "### Visualizing the Protein-Ligand interaction residues\n",
    "\n",
    "Viewing the residues involved in the interaction with NGL. **Ligand** is represented in CPK, whereas **protein residues** are represented in Licorice. \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(cmipPDB)\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='all',color='sstruc',opacity=0.3)\n",
    "view.add_representation(repr_type='ball+stick', selection=ligandCode)\n",
    "view.add_representation(repr_type='licorice', selection=identified_residues_ngl)\n",
    "view.center(ligandCode)\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"prot-prot\"></a>\n",
    "***\n",
    "## Protein-Protein Interaction Energies\n",
    "\n",
    "This example illustrates the steps needed to compute the **protein-protein interaction energies** from a **structure conformation** of the **complex** taken directly from a **Molecular Dynamics simulation**. The particular example used is the **SARS-CoV-2 Receptor Binding Domain (RBD)** complexed with the human **Angiotensin Converting Enzyme 2 (hACE2)** (PDB code [6VW1](https://www.rcsb.org/structure/6vw1), MD code [MCV1900410](https://bioexcel-cv19.bsc.es/#/id/MCV1900410/overview), [https://doi.org/10.2210/pdb6VW1/pdb](https://doi.org/10.2210/pdb6VW1/pdb)). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visRBD-hACE2\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the original **SARS-CoV-2 Receptor Binding Domain** (blue chain) and the **human Angiotensin Converting Enzyme 2** (red chain) structure complex, taken from a **MD simulation trajectory**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(inputPDB_MD)\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='protein', color='chainname')\n",
    "view._remote_call('setSize', target='Widget', args=['','400px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"preparePDB_MD\"></a>\n",
    "***\n",
    "## CMIP PDB Preparation (from Molecular Dynamics topology)\n",
    "When working with **structure conformations** taken from **MD simulations**, it is recommended to use the **charges, atom types and elements** considered in the simulation. Usually this information is stored in the so-called ***topology*** files. A specific **building block** (***prepare_structure***) is available to extract these information from an **MD topology file** and use it for the **CMIP calculations**.\n",
    " \n",
    "The next cells are taking one frame of the **MD simulation**, splitting the subunits (chains) in two different **PDB files**, and preparing them to be used in **CMIP**, taking the **charges and elements** used in the simulations from the **MD topology file**. \n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [prepare_structure](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.prepare_structure) from **biobb_cmip.cmip.prepare_structure** \n",
    " - [extract_chain](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.extract_chain) from **biobb_structure_utils.utils.extract_chain**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing structure\n",
    "Preparing the structure for the **CMIP calculations**, using the original **MD topology file** information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_prepare_structure import cmip_prepare_structure\n",
    "\n",
    "cmipPDB_MD = MDCode + \".cmip.pdb\"\n",
    "\n",
    "cmip_prepare_structure(input_pdb_path=inputPDB_MD,\n",
    "                       input_topology_path=inputTOP_MD,\n",
    "                       output_cmip_pdb_path=cmipPDB_MD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting chains \n",
    "Saving **hACE2** (chain A) and **RBD** (chain B) in two different **PDB files**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_structure_utils.utils.extract_chain import extract_chain\n",
    "\n",
    "cmipPDB_MD_hACE2 = MDCode + \".hACE2.cmip.pdb\"\n",
    "cmipPDB_MD_RBD = MDCode + \".RBD.cmip.pdb\"\n",
    "\n",
    "prop = {\n",
    "    'chains': [ 'A' ],\n",
    "    'permissive' : True    \n",
    "}\n",
    "extract_chain(input_structure_path=cmipPDB_MD,\n",
    "            output_structure_path=cmipPDB_MD_hACE2,\n",
    "            properties=prop)\n",
    "\n",
    "prop = {\n",
    "    'chains': [ 'B' ],\n",
    "    'permissive' : True\n",
    "}\n",
    "extract_chain(input_structure_path=cmipPDB_MD,\n",
    "            output_structure_path=cmipPDB_MD_RBD,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing structures\n",
    "Visualizing the generated structures: **hACE2** and **RBD**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show different structures generated \n",
    "view1 = nglview.show_structure_file(cmipPDB_MD_hACE2)\n",
    "view1.component_0.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view1.component_0.center()\n",
    "view1._remote_call('setSize', target='Widget', args=['400px','400px'])\n",
    "view1.camera='orthographic'\n",
    "view1\n",
    "view2 = nglview.show_structure_file(cmipPDB_MD_RBD)\n",
    "view2.component_0.add_representation(repr_type='cartoon', selection='protein', color='sstruc')\n",
    "view2.component_0.center()\n",
    "view2._remote_call('setSize', target='Widget', args=['400px','400px'])\n",
    "view2.camera='orthographic'\n",
    "view2\n",
    "ipywidgets.HBox([view1, view2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_boxes\"></a>\n",
    "***\n",
    "## CMIP Boxes\n",
    "\n",
    "**CMIP** computes **interaction energies** using **3D grids** (see image). **CMIP** tries to automatically built a **grid** enclosing the  structure under study, which is fine in most of the cases. However, when working with **protein-protein interactions**, the **grid** needs to be adjusted to the system. The next cells are computing and showing the **3D grids** needed to enclose the two **monomers** and the **complex**, and are useful to inspect them and re-define them, increasing the size, if needed.\n",
    "\n",
    "A specific property of the **cmip building block** (*execution_type = check_only*) allows the calculation of the **3D grids**, and accepts params such as *perfill*, which specifies the **portion of the grid** filled by the input (host) structure or *box_size_factor*, which increases the **size of the grid**. The computed **grids** can be then used as an input for the final **interaction energies** calculations.\n",
    "\n",
    "***\n",
    "<img src=\"Files/grid2.png\" alt=\"3D Grid\" title=\"3D Grid\" width=\"200\" /></img>\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_run](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.cmip) from **biobb_cmip.cmip.cmip_run**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box (3D Grid) for the RBD\n",
    "Computing and visualizing the **box** enclosing the **RBD monomer**. Applying a *perfill* factor of 0.8. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "cmip_RBD_box_log = \"RBD.box.log\"\n",
    "cmip_RBD_box_out = \"RBD.box.byat.out\"\n",
    "cmip_RBD_box_json = \"RBD.box.json\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'check_only',\n",
    "    'params' : {\n",
    "         'perfill' : 0.8\n",
    "    }\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_MD_RBD,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_log_path=cmip_RBD_box_log,\n",
    "         output_json_box_path=cmip_RBD_box_json,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the RBD Box\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "from biobb_cmip.utils.representation import create_box_representation\n",
    "\n",
    "boxedFilename, atomPair = create_box_representation(cmip_RBD_box_log, inputPDB_MD)\n",
    "# Represent the new file in ngl\n",
    "view = nv.show_structure_file(boxedFilename, default=False)\n",
    "# Structure\n",
    "view.add_representation(repr_type='cartoon', selection='not het', color='#cccccc', opacity=.2)\n",
    "# vertices box\n",
    "view.add_representation(repr_type='ball+stick', selection='9999', aspectRatio = 10)\n",
    "# lines box\n",
    "view.add_representation(repr_type='distance', atomPair= atomPair, labelColor= 'transparent', color= 'black')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box (3D Grid) for the hACE2\n",
    "Computing and visualizing the **box** enclosing the **hACE2 monomer**. Applying a *perfill* factor of 0.8. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "cmip_hACE2_box_log = \"hACE2.box.log\"\n",
    "cmip_hACE2_box_out = \"hACE2.box.byat.out\"\n",
    "cmip_hACE2_box_json = \"hACE2.box.json\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'check_only',\n",
    "    'params' : {\n",
    "         'perfill' : 0.8,\n",
    "    }\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_MD_hACE2,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_log_path=cmip_hACE2_box_log,\n",
    "         output_json_box_path=cmip_hACE2_box_json,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the hACE2 box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "from biobb_cmip.utils.representation import create_box_representation\n",
    "\n",
    "boxedFilename, atomPair = create_box_representation(cmip_hACE2_box_log, inputPDB_MD)\n",
    "# Represent the new file in ngl\n",
    "view = nv.show_structure_file(boxedFilename, default=False)\n",
    "# Structure\n",
    "view.add_representation(repr_type='cartoon', selection='not het', color='#cccccc', opacity=.2)\n",
    "# vertices box\n",
    "view.add_representation(repr_type='ball+stick', selection='9999', aspectRatio = 10)\n",
    "# lines box\n",
    "view.add_representation(repr_type='distance', atomPair= atomPair, labelColor= 'transparent', color= 'black')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box (3D Grid) for the RBD-hACE2 complex\n",
    "Computing and visualizing the **box** enclosing the **RBD-hACE2 complex**. Applying a *perfill* factor of 0.6. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "cmip_COMPLEX_box_log = \"COMPLEX.box.log\"\n",
    "cmip_COMPLEX_box_out = \"COMPLEX.box.byat.out\"\n",
    "cmip_COMPLEX_box_json = \"COMPLEX.box.json\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'check_only',\n",
    "    'params' : {\n",
    "         'perfill' : 0.6,\n",
    "    }\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_MD,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_log_path=cmip_COMPLEX_box_log,\n",
    "         output_json_box_path=cmip_COMPLEX_box_json,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the RBD-hACE2 complex box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "from biobb_cmip.utils.representation import create_box_representation\n",
    "\n",
    "boxedFilename, atomPair = create_box_representation(cmip_COMPLEX_box_log, inputPDB_MD)\n",
    "# Represent the new file in ngl\n",
    "view = nv.show_structure_file(boxedFilename, default=False)\n",
    "# Structure\n",
    "view.add_representation(repr_type='cartoon', selection='not het', color='#cccccc', opacity=.2)\n",
    "# vertices box\n",
    "view.add_representation(repr_type='ball+stick', selection='9999', aspectRatio = 10)\n",
    "# lines box\n",
    "view.add_representation(repr_type='distance', atomPair= atomPair, labelColor= 'transparent', color= 'black')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"RBDinteraction\"></a>\n",
    "***\n",
    "## RDB Interaction Potential Energies \n",
    "The first analysis computes the **interaction potential energies** for the **RBD atoms** with respect to the **hACE2 enzyme**. \n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_ignore_residues](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.ignore_residues) from **biobb_cmip.cmip.cmip_ignore_residues**\n",
    " - [cmip_run](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.cmip) from **biobb_cmip.cmip.cmip_run**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CMIP PDB structure Preparation\n",
    "\n",
    "Creating the **RBD-hACE2 Complex** with **RBD** in **dielectric mode**. For this, all residues belonging to the **RBD** (in this case, the chain B) are marked with an 'X' symbol in the PDB. This is done by the *ignore_residues* building block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_ignore_residues import cmip_ignore_residues\n",
    "\n",
    "cmipPDB_MD_RBD_ignored = MDCode + \".RBD_ignored.cmip.pdb\"\n",
    "\n",
    "prop = {\n",
    "    'residue_list': \"B:\"\n",
    "}\n",
    "\n",
    "cmip_ignore_residues(input_cmip_pdb_path = cmipPDB_MD,\n",
    "               output_cmip_pdb_path = cmipPDB_MD_RBD_ignored,\n",
    "               properties = prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_protein_RBD\"></a>\n",
    "***\n",
    "### Computing the Protein-Protein interaction energies\n",
    "\n",
    "Computing the **Protein-Protein interaction energies** using the **RBD-hACE2** complex with the **RBD** in **dielectric mode** as a **CMIP host structure** and the **RBD** as a **CMIP ligand** (probe). \n",
    "\n",
    "**Grids sizes** previously computed for the **RBD** (*input_json_box_path*) and for the **RBD-hACE2 complex** (*input_json_external_box_path*) are used as input information for the **CMIP** calculation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "RBD_energies_log = \"RBD.energies.log\"\n",
    "RBD_byat_out = \"RBD.energies.byat.out\"\n",
    "output_RBD_box_json = \"RBD.energies.box.output.json\"\n",
    "output_COMPLEX_box_json = \"COMPLEX.energies.box.output.json\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'pb_interaction_energy'\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_MD_RBD_ignored,\n",
    "         input_probe_pdb_path=cmipPDB_MD_RBD,\n",
    "         input_json_box_path=cmip_RBD_box_json,\n",
    "         input_json_external_box_path=cmip_COMPLEX_box_json,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_json_box_path=output_RBD_box_json,\n",
    "         output_json_external_box_path=output_COMPLEX_box_json,\n",
    "         output_log_path=RBD_energies_log,\n",
    "         output_byat_path=RBD_byat_out,\n",
    "         properties=prop)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visBOX1\"></a>\n",
    "### Visualizing CMIP Grid Box\n",
    "Visualizing the **grid box** used by **CMIP** to compute the **Interaction Potential Energies** (taken from the log file). It is important to check that the **grid box** includes the whole **interaction region**, which is the region of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "from biobb_cmip.utils.representation import create_box_representation\n",
    "\n",
    "boxedFilename, atomPair = create_box_representation(output_COMPLEX_box_json, inputPDB_MD)\n",
    "# Represent the new file in ngl\n",
    "view = nv.show_structure_file(boxedFilename, default=False)\n",
    "# Structure\n",
    "view.add_representation(repr_type='cartoon', selection='not het', color='#cccccc', opacity=.2)\n",
    "# vertices box\n",
    "view.add_representation(repr_type='ball+stick', selection='9999', aspectRatio = 10)\n",
    "# lines box\n",
    "view.add_representation(repr_type='distance', atomPair= atomPair, labelColor= 'transparent', color= 'black')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"plotRBD_atoms\"></a>\n",
    "### Interaction energies by atom\n",
    "Visualizing the **interaction potential energies** computed by **CMIP**. The plot shows **interactions energies** (in kcal/mol, Y axis) for **each of the atoms** of the **RBD protein** (X axis). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "from biobb_cmip.utils.representation import get_energies_byat\n",
    "\n",
    "atom_list, energy_dict = get_energies_byat(RBD_byat_out, cutoff=55)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=atom_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Atom Number\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"plotRBD_residues\"></a>\n",
    "### Interaction energies by residue\n",
    "Visualizing the **interaction potential energies** computed by **CMIP**. The plot shows **interactions energies** (in kcal/mol, Y axis) for **each of the residues** (computed summing the contributions of all atoms included in the residue) of the **RBD protein** (X axis). \n",
    "\n",
    "NOTE: This plot reproduces the analysis (for one single snapshot) presented here: [https://bioexcel-cv19.bsc.es/#/id/MCV1900002/energies](https://bioexcel-cv19.bsc.es/#/id/MCV1900002/energies) with **Host: RBD and Guest: hACE2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "from biobb_cmip.utils.representation import get_energies_byres\n",
    "\n",
    "res_list, energy_dict = get_energies_byres(RBD_byat_out, cutoff=55)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=res_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Residue ID\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()\n",
    "\n",
    "energy_cutoff = -2\n",
    "identified_residues = [res_list[energy_dict['ES&VDW'].index(item)] \n",
    "                       for item in energy_dict['ES&VDW'] if item < energy_cutoff]\n",
    "identified_residues_ngl = str([int(re.sub(r'[A-Z]+', '', item)) for item in identified_residues]).replace(',','')\n",
    "\n",
    "print (\"Residues involved in the interaction (energy < \" + str(energy_cutoff) + \" KCal/mol):\")\n",
    "print (identified_residues_ngl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_protein_RBD_ngl\"></a>\n",
    "***\n",
    "### Visualizing the RBD residues involved in the interaction with hACE2\n",
    " \n",
    "Viewing the **RBD residues** involved in the **interaction** with NGL. Residues in the interaction are represented in CPK (Ball & Sticks). Residues involved in the interaction have been selected in the previous cell as the ones having an **interaction energy lower than a particular cutoff (-2 Kcal/mol)**. The **cutoff** can be changed to adapt it to the particular system.  \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(cmipPDB_MD)\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='all',color='sstruc',opacity=0.3)\n",
    "view.add_representation(repr_type='licorice', radius='0.5', selection=identified_residues_ngl)\n",
    "view.add_representation(repr_type='surface', selection=identified_residues_ngl, opacity=0.2, color='element')\n",
    "view.center(identified_residues_ngl)\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"hACE2interaction\"></a>\n",
    "***\n",
    "## hACE2 Interaction Potential Energies \n",
    "The second analysis computes the **interaction potential energies** for the **hACE2 atoms** with respect to the **RBD receptor**. \n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cmip_ignore_residues](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.ignore_residues) from **biobb_cmip.cmip.cmip_ignore_residues**\n",
    " - [cmip_run](https://biobb-cmip.readthedocs.io/en/latest/cmip.html#module-cmip.cmip) from **biobb_cmip.cmip.cmip_run**\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CMIP PDB structure Preparation\n",
    "\n",
    "Creating the **RBD-hACE2 Complex** with **hACE2** in **dielectric mode**. For this, all residues belonging to the **hACE2** (in this case, the chain A) are marked with an 'X' symbol in the PDB. This is done by the *ignore_residues* building block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_ignore_residues import cmip_ignore_residues\n",
    "\n",
    "cmipPDB_MD_hACE2_ignored = MDCode + \".hACE2_ign.cmip.pdb\"\n",
    "\n",
    "prop = {\n",
    "    'residue_list': \"A:\"\n",
    "}\n",
    "\n",
    "cmip_ignore_residues(input_cmip_pdb_path = cmipPDB_MD,\n",
    "                     output_cmip_pdb_path = cmipPDB_MD_hACE2_ignored,\n",
    "                     properties = prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_protein_hACE2\"></a>\n",
    "***\n",
    "### Computing the Protein-Protein interaction energies\n",
    "\n",
    "Computing the **Protein-Protein interaction energies** using the complex with the **hACE2** in **dielectric mode** as a **CMIP host structure** and the **hACE2** as a **CMIP ligand** (probe). \n",
    "\n",
    "**Grids sizes** previously computed for the **hACE2** (*input_json_box_path*) and for the **RBD-hACE2 complex** (*input_json_external_box_path*) are used as input information for the **CMIP** calculation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_cmip.cmip.cmip_run import cmip_run\n",
    "\n",
    "hACE2_energies_log = \"hACE2.energies.log\"\n",
    "hACE2_byat_out = \"hACE2.energies.byat.out\"\n",
    "output_hACE2_box_json = \"hACE2.energies.box.output.json\"\n",
    "output_COMPLEX_2_box_json = \"COMPLEX_2.energies.box.output.json\"\n",
    "\n",
    "prop = { \n",
    "    'execution_type' : 'pb_interaction_energy'\n",
    "}\n",
    "\n",
    "cmip_run(input_pdb_path=cmipPDB_MD_hACE2_ignored,\n",
    "         input_probe_pdb_path=cmipPDB_MD_hACE2,\n",
    "         input_json_box_path=cmip_hACE2_box_json,\n",
    "         input_json_box_external_path=cmip_COMPLEX_box_json,\n",
    "         input_vdw_params_path=vdw_params_path,\n",
    "         output_json_box_path=output_hACE2_box_json,\n",
    "         output_json_external_box_path=output_COMPLEX_2_box_json,\n",
    "         output_log_path=hACE2_energies_log,\n",
    "         output_byat_path=hACE2_byat_out,\n",
    "         properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"visBOX2\"></a>\n",
    "### Visualizing CMIP Box\n",
    "Visualizing the **box** used by **CMIP** to compute the **Interaction Potential Energies** (taken from the log file). It is important to check that the box includes the whole **interaction region**, which is the region of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "from biobb_cmip.utils.representation import create_box_representation\n",
    "\n",
    "boxedFilename, atomPair = create_box_representation(hACE2_energies_log, inputPDB_MD)\n",
    "# Represent the new file in ngl\n",
    "view = nv.show_structure_file(boxedFilename, default=False)\n",
    "# Structure\n",
    "view.add_representation(repr_type='cartoon', selection='not het', color='#cccccc', opacity=.2)\n",
    "# vertices box\n",
    "view.add_representation(repr_type='ball+stick', selection='9999', aspectRatio = 10)\n",
    "# lines box\n",
    "view.add_representation(repr_type='distance', atomPair= atomPair, labelColor= 'transparent', color= 'black')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"plotRBD_atoms\"></a>\n",
    "### Interaction energies by atom\n",
    "Visualizing the **interaction potential energies** computed by **CMIP**. The plot shows **interactions energies** (in kcal/mol, Y axis) for **each of the atoms** of the **hACE2 protein** (X axis). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "from biobb_cmip.utils.representation import get_energies_byat\n",
    "\n",
    "atom_list, energy_dict = get_energies_byat(hACE2_byat_out, cutoff=55)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=atom_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Atom Number\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"plotRBD_residues\"></a>\n",
    "### Interaction energies by residue\n",
    "Visualizing the **interaction potential energies** computed by **CMIP**. The plot shows **interactions energies** (in kcal/mol, Y axis) for **each of the residues** (computed summing the contributions of all atoms included in the residue) of the **hACE2 protein** (X axis). \n",
    "\n",
    "NOTE: This plot reproduces the analysis (for one single snapshot) presented here: [https://bioexcel-cv19.bsc.es/#/id/MCV1900002/energies](https://bioexcel-cv19.bsc.es/#/id/MCV1900002/energies) with **Host: hACE2 and Guest: RBD**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "from biobb_cmip.utils.representation import get_energies_byres\n",
    "\n",
    "res_list, energy_dict = get_energies_byres(hACE2_byat_out, cutoff=55)\n",
    "\n",
    "# Create a scatter plot\n",
    "fig = go.Figure(data=go.Scatter(x=res_list, y=energy_dict['ES&VDW'], mode='lines'))\n",
    "\n",
    "# Update layout\n",
    "fig.update_layout(title=\"CMIP Interaction Potential\",\n",
    "                  xaxis_title=\"Residue ID\",\n",
    "                  yaxis_title=\"Potential Energy Kcal/mol\",\n",
    "                  height=600)\n",
    "\n",
    "# Show the figure\n",
    "fig.show()\n",
    "\n",
    "energy_cutoff = -1.5\n",
    "identified_residues2 = [res_list[energy_dict['ES&VDW'].index(item)] \n",
    "                       for item in energy_dict['ES&VDW'] if item < energy_cutoff]\n",
    "identified_residues2_ngl = str([int(re.sub(r'[A-Z]+', '', item)) for item in identified_residues2]).replace(',','')\n",
    "\n",
    "print (\"Residues involved in the interaction (energy < \" + str(energy_cutoff) + \" KCal/mol):\")\n",
    "print (identified_residues_ngl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cmip_protein_protein_hACE2_ngl\"></a>\n",
    "***\n",
    "### Visualizing the hACE2 residues involved in the interaction with RBD\n",
    "\n",
    "Viewing the **hACE2 residues** involved in the **interaction** with NGL. Residues in the interaction are represented in CPK (Ball & Sticks). Residues involved in the interaction have been selected in the previous cell as the ones having an **interaction energy lower than a particular cutoff (-1.5 Kcal/mol)**. The **cutoff** can be changed to adapt it to the particular system.  \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show protein\n",
    "view = nglview.show_structure_file(cmipPDB_MD)\n",
    "view.clear_representations()\n",
    "view.add_representation(repr_type='cartoon', selection='all',color='sstruc',opacity=0.3)\n",
    "view.add_representation(repr_type='ball+stick', selection=ligandCode)\n",
    "view.add_representation(repr_type='licorice', radius='0.5', selection=identified_residues2_ngl)\n",
    "view.add_representation(repr_type='surface', selection=identified_residues2_ngl, opacity=0.3, color='element')\n",
    "view.center(identified_residues_ngl)\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "<a id=\"questions\"></a>\n",
    "\n",
    "## Questions & Comments\n",
    "\n",
    "Questions, issues, suggestions and comments are really welcome!\n",
    "\n",
    "* GitHub issues:\n",
    "    * [https://github.com/bioexcel/biobb](https://github.com/bioexcel/biobb)\n",
    "\n",
    "* BioExcel forum:\n",
    "    * [https://ask.bioexcel.eu/c/BioExcel-Building-Blocks-library](https://ask.bioexcel.eu/c/BioExcel-Building-Blocks-library)\n"
   ]
  }
 ],
 "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
