Glycam.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "name": "Glycam.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/pablo-arantes/making-it-rain/blob/main/Glycam.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pj2BWZxUDbEE"
      },
      "source": [
        "# **Hello there!**\n",
        "\n",
        "This is a Jupyter notebook for running Molecular Dynamics (MD) simulations using OpenMM engine and AMBER/GLYCAM force field for **Protein** and **Glycosylated** systems. This notebook is a supplementary material of the paper \"***Making it rain: Cloud-based molecular simulations for everyone***\" ([link here](https://doi.org/10.1021/acs.jcim.1c00998)) and we encourage you to read it before using this pipeline.\n",
        "\n",
        "The main goal of this notebook is to demonstrate how to harness the power of cloud-computing to run microsecond-long MD simulations in a cheap and yet feasible fashion.\n",
        "\n",
        "---\n",
        "\n",
        " **This notebook is NOT a standard protocol for MD simulations!** It is just a simple MD pipeline illustrating each step of a simulation protocol.\n",
        "\n",
        "--- \n",
        "**Bugs**\n",
        "- If you encounter any bugs, please report the issue to https://github.com/pablo-arantes/making-it-rain/issues\n",
        "\n",
        "**Acknowledgments**\n",
        "- We would like to thank the OpenMM team for developing an excellent and open source engine. \n",
        "\n",
        "- Thank you to [GLYCAM-web](https://glycam.org/#) server team, Professor Robert J. Woods, for their amazing web-tool and force field.\n",
        "\n",
        "- Making-it-rain by **Pablo R. Arantes** ([@pablitoarantes](https://twitter.com/pablitoarantes)), **Marcelo D. Polêto** ([@mdpoleto](https://twitter.com/mdpoleto)), **Conrado Pedebos** ([@ConradoPedebos](https://twitter.com/ConradoPedebos)) and **Rodrigo Ligabue-Braun** ([@ligabue_braun](https://twitter.com/ligabue_braun)).\n",
        "\n",
        "\n",
        "- Also, credit to [David Koes](https://github.com/dkoes) for his awesome [py3Dmol](https://3dmol.csb.pitt.edu/) plugin.\n",
        "\n",
        "- For related notebooks see: [Making-it-rain](https://github.com/pablo-arantes/making-it-rain)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hoyY6XonD1UX"
      },
      "source": [
        "# **Introduction**\n",
        "\n",
        "In general, MD simulations rely on 1) a set of atomic coordinates of all atoms on a simulation box and 2) a set of force field parameters that describes the interaction energies between atoms.\n",
        "\n",
        "In terms of inputs, we wil need:\n",
        "*  A PDB ID or a .pdb file containing a set of atomic coordinates of your glycosylated system from [GLYCAM](https://glycam.org/#) server.\n",
        "\n",
        "In this notebook, we will simulate a glycosylated protein, Ferrulate esterase (PDB ID 2HL6). To build our simulation box, we will use LEaP program (https://ambermd.org/tutorials/pengfei/index.php). The LEaP program is a portal between many chemical structure file types (.pdb and .mol2, primarily), and the Amber model parameter file types such as .lib, .prepi, parm.dat, and .frcmod. Each of the parameter files contains pieces of information needed for constructing a simulation, whether for energy minimization or molecular dynamics. LEaP functions within a larger workflow described in Section 1.1 of the [Amber Manual](https://ambermd.org/doc12/Amber20.pdf);.\n",
        "\n",
        "You can download the input files examples from [here](https://github.com/pablo-arantes/making-it-rain/tree/main/AMBER_GLYCAM);\n",
        "## ---\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lh96y6mGFY1D"
      },
      "source": [
        "---\n",
        "---\n",
        "# **Setting the environment for MD calculation**\n",
        "\n",
        "Firstly, we need to install all necessary libraries and packages for our simulation. The main packages we will be installing are:\n",
        "\n",
        "1.    Anaconda (https://docs.conda.io/en/latest/miniconda.html)\n",
        "2.    OpenMM (https://openmm.org/)\n",
        "3.    PyTraj (https://amber-md.github.io/pytraj/latest/index.html)\n",
        "4.    py3Dmol (https://pypi.org/project/py3Dmol/)\n",
        "5.    Numpy (https://numpy.org/)\n",
        "6.    Matplotlib (https://matplotlib.org/)\n",
        "7.    AmberTools (https://ambermd.org/AmberTools.php)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wH1oMiVUlxO_",
        "cellView": "form"
      },
      "source": [
        "#@title **Install dependencies**\n",
        "#@markdown It will take a few minutes, please, drink a coffee and wait. ;-)\n",
        "# install dependencies\n",
        "!pip -q install py3Dmol 2>&1 1>/dev/null\n",
        "!pip install --upgrade MDAnalysis 2>&1 1>/dev/null\n",
        "!pip install git+https://github.com/pablo-arantes/biopandas 2>&1 1>/dev/null\n",
        "# install conda\n",
        "!wget -qnc https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh \n",
        "!bash Miniconda3-latest-Linux-x86_64.sh -bfp /usr/local 2>&1 1>/dev/null\n",
        "!rm Miniconda3-latest-Linux-x86_64.sh\n",
        "!conda install -y -q -c conda-forge openmm=7.6 python=3.7 pdbfixer 2>&1 1>/dev/null\n",
        "!conda install -c conda-forge ambertools --yes 2>&1 1>/dev/null\n",
        "!conda install -c ambermd pytraj  --yes 2>&1 1>/dev/null\n",
        "#load dependencies\n",
        "import sys\n",
        "sys.path.append('/usr/local/lib/python3.7/site-packages/')\n",
        "from biopandas.pdb import PandasPdb\n",
        "import openmm as mm\n",
        "from openmm import *\n",
        "from openmm.app import *\n",
        "from openmm.unit import *\n",
        "import os\n",
        "import urllib.request  \n",
        "import numpy as np\n",
        "import MDAnalysis as mda\n",
        "import py3Dmol\n",
        "from __future__ import print_function\n",
        "import pytraj as pt\n",
        "import platform\n",
        "import scipy.cluster.hierarchy\n",
        "from scipy.spatial.distance import squareform\n",
        "import scipy.stats as stats\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "from scipy.interpolate import griddata\n",
        "import seaborn as sb\n",
        "from statistics import mean, stdev\n",
        "from pytraj import matrix\n",
        "from matplotlib import colors\n",
        "from IPython.display import set_matplotlib_formats\n",
        "#%matplotlib inline\n",
        "#set_matplotlib_formats('png')\n",
        "#plt.figure(figsize=(5,7))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fDQnAKJLFxtt"
      },
      "source": [
        "## Using Google Drive to store simulation data\n",
        "\n",
        "Google Colab does not allow users to keep data on their computing nodes. However, we can use Google Drive to read, write, and store our simulations files. Therefore, we suggest to you to:\n",
        "\n",
        "1.   Create a folder in your own Google Drive and copy the necessary input files there.\n",
        "2.   Copy the path of your created directory. We will use it below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Lm7Akepv_vl-",
        "cellView": "form"
      },
      "source": [
        "#@title ### **Import Google Drive**\n",
        "#@markdown Click in the \"Run\" buttom to make your Google Drive accessible.\n",
        "from google.colab import drive\n",
        "\n",
        "drive.flush_and_unmount()\n",
        "drive.mount('/content/drive', force_remount=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lOKg9eH_ueRn",
        "cellView": "form"
      },
      "source": [
        "#@title **Check if you correctly allocated GPU nodes**\n",
        "\n",
        "gpu_info = !nvidia-smi\n",
        "gpu_info = '\\n'.join(gpu_info)\n",
        "if gpu_info.find('failed') >= 0:\n",
        "  print('Select the Runtime > \"Change runtime type\" menu to enable a GPU accelerator, ')\n",
        "  print('and then re-execute this cell.')\n",
        "else:\n",
        "  print(gpu_info)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BZwl66HTGI7v"
      },
      "source": [
        "---\n",
        "---\n",
        "# **Loading the necessary input files**\n",
        "\n",
        "At this point, we should have all libraries and dependencies installed and all necessary input files already at your Google Drive folder.\n",
        "\n",
        "**Important**: Make sure the PDB file points to the correct pathway. If necessary, correct the pathway and re-upload the files.\n",
        "\n",
        "Below, you should provide the names of all input files and the pathway of your Google Drive folder containing them.\n",
        "\n",
        "**Please, don't use spaces in the files and folders names, use _ instead i.e., protein_glycan.pdb, MyDrive/protein_glycan and so on.**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wZfALbPc6t5H",
        "cellView": "form"
      },
      "source": [
        "#@title **Please, provide the necessary input files below**: \n",
        "#@markdown **Important:** This notebook uses as input, the pdb format from AMBER/GLYCAM package. You **MUST** use the pdb from GLYCAM server (https://glycam.org/). Probably you will have an error in topology generation, if you use a different PDB format.\n",
        "PDB_file_name = 'structure.pdb' #@param {type:\"string\"}\n",
        "file_name = PDB_file_name\n",
        "Google_Drive_Path = '/content/drive/MyDrive/protein_glycan' #@param {type:\"string\"}\n",
        "workDir = Google_Drive_Path\n",
        "initial_pdb = os.path.join(workDir, str(file_name))\n",
        "# starting = os.path.join(workDir, \"starting.pdb\")\n",
        "# starting2 = os.path.join(workDir, str(file_name))\n",
        "starting_end = os.path.join(workDir, str(file_name))\n",
        "tleap = os.path.join(workDir, \"tleap.in\")\n",
        "top_nw = os.path.join(workDir, \"SYS_nw.prmtop\")\n",
        "crd_nw = os.path.join(workDir, \"SYS_nw.crd\")\n",
        "pdb_nw = os.path.join(workDir, \"SYS_nw.pdb\")\n",
        "top = os.path.join(workDir, \"SYS.prmtop\")\n",
        "crd = os.path.join(workDir, \"SYS.crd\")\n",
        "pdb = os.path.join(workDir, \"SYS.pdb\")\n",
        "\n",
        "# pdb4amber_cmd = \"pdb4amber -i \" + str(starting2) + \" -o \" + str(starting_end) + \" -p\"\n",
        "\n",
        "# original_stdout = sys.stdout # Save a reference to the original standard output\n",
        "\n",
        "# with open('pdb4amber.sh', 'w') as f:\n",
        "#     sys.stdout = f # Change the standard output to the file we created.\n",
        "#     print(pdb4amber_cmd)\n",
        "#     sys.stdout = original_stdout # Reset the standard output to its original value\n",
        "\n",
        "# !chmod 700 pdb4amber.sh 2>&1 1>/dev/null\n",
        "# !bash pdb4amber.sh\n",
        "\n",
        "#@markdown ---"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CU7wej8xdjrW",
        "cellView": "form"
      },
      "source": [
        "#@title **Parameters to generate the Amber topology:**\n",
        "\n",
        "Force_field = \"ff19SB\" #@param [\"ff19SB\", \"ff14SB\"]\n",
        "if Force_field == \"ff19SB\":\n",
        "  ff = \"leaprc.protein.ff19SB\"\n",
        "else:\n",
        "  ff = \"leaprc.protein.ff14SB\"\n",
        "\n",
        "Water_type = \"TIP3P\" #@param [\"TIP3P\", \"OPC\"]\n",
        "if Water_type == \"TIP3P\":\n",
        "  water = \"leaprc.water.tip3p\"\n",
        "  water_box = \"TIP3PBOX\"\n",
        "else:\n",
        "  water = \"leaprc.water.opc\"\n",
        "  water_box = \"OPCBOX\"\n",
        "\n",
        "#@markdown Size Box (Angstrons):\n",
        "\n",
        "Size_box = 12 #@param {type:\"slider\", min:10, max:20, step:1}\n",
        "size_box = Size_box\n",
        "\n",
        "#@markdown **ATTENTION**: Give the concentration in Molar units, AMBER tleap will neutralize your system automatically:\n",
        "\n",
        "Ions = \"NaCl\" #@param [\"NaCl\", \"KCl\" ]\n",
        "\n",
        "Concentration = \"0.15\" #@param {type:\"string\"}\n",
        "\n",
        "#@markdown ---\n",
        "\n",
        "f = open(tleap, \"w\")\n",
        "f.write(\"\"\"source \"\"\" + str(ff) + \"\\n\"\n",
        "\"\"\"source leaprc.DNA.OL15\n",
        "source leaprc.RNA.OL3\n",
        "source leaprc.GLYCAM_06j-1 \n",
        "source leaprc.gaff2\n",
        "source \"\"\"  + str(water) + \"\\n\"\n",
        "\"\"\"SYS = loadpdb \"\"\"  + str(starting_end) + \"\\n\"\n",
        "\"\"\"alignaxes SYS\n",
        "savepdb SYS \"\"\" + str(pdb_nw) + \"\\n\"\n",
        "\"\"\"saveamberparm SYS \"\"\" + str(top_nw) + \" \" + str(crd_nw) + \"\\n\"\n",
        "\"\"\"solvatebox SYS \"\"\" + str(water_box) + \" \" + str(size_box) +  \"\"\" 0.7\n",
        "saveamberparm SYS \"\"\" + str(top) + \" \" + str(crd) + \"\\n\"\n",
        "\"\"\"savepdb SYS \"\"\" + str(pdb) + \"\\n\"\n",
        "\"\"\"quit\"\"\")\n",
        "f.close()\n",
        "\n",
        "tleap_command = \"tleap -f \" + str(tleap)\n",
        "\n",
        "original_stdout = sys.stdout # Save a reference to the original standard output\n",
        "\n",
        "with open('run_tleap.sh', 'w') as f:\n",
        "    sys.stdout = f # Change the standard output to the file we created.\n",
        "    print(tleap_command)\n",
        "    sys.stdout = original_stdout # Reset the standard output to its original value\n",
        "\n",
        "SYS = os.path.join(workDir, \"SYS*\")\n",
        "rm_sys = \"rm \" + SYS\n",
        "\n",
        "original_stdout = sys.stdout # Save a reference to the original standard output\n",
        "\n",
        "with open('rm_sys.sh', 'w') as f:\n",
        "    sys.stdout = f # Change the standard output to the file we created.\n",
        "    print(rm_sys)\n",
        "    sys.stdout = original_stdout # Reset the standard output to its original value\n",
        "\n",
        "!chmod 700 rm_sys.sh 2>&1 1>/dev/null\n",
        "!bash rm_sys.sh 2> /dev/null\n",
        "\n",
        "!chmod 700 run_tleap.sh 2>&1 1>/dev/null\n",
        "!bash run_tleap.sh 2>&1 1>/dev/null\n",
        "\n",
        "!grep \"Volume:\" leap.log > temp.txt\n",
        "with open(\"temp.txt\", 'r') as f:\n",
        "  for line in f:\n",
        "        vol = float(line.split()[1])\n",
        "\n",
        "vol_lit  = vol * pow(10, -27)\n",
        "atom_lit = 9.03 * pow(10, 22)\n",
        "conc = float(Concentration)\n",
        "num_ion = int(vol_lit * (conc/0.15) * atom_lit)\n",
        "\n",
        "if Ions == \"NaCl\":\n",
        "  pos_neut = \"Na+ 0\"\n",
        "  pos_num = \"Na+ \" + str(num_ion)\n",
        "  Cl_num = num_ion\n",
        "else:\n",
        "  pos_neut = \"K+ 0\"\n",
        "  pos_num = \"K+ \" + str(num_ion)\n",
        "  Cl_num = num_ion\n",
        "\n",
        "f = open(tleap, \"w\")\n",
        "f.write(\"\"\"source \"\"\" + str(ff) + \"\\n\"\n",
        "\"\"\"source leaprc.DNA.OL15\n",
        "source leaprc.RNA.OL3\n",
        "source leaprc.GLYCAM_06j-1 \n",
        "source leaprc.gaff2\n",
        "source \"\"\"  + str(water) + \"\\n\"\n",
        "\"\"\"SYS = loadpdb \"\"\"  + str(starting_end) + \"\\n\"\n",
        "\"\"\"alignaxes SYS\n",
        "check SYS \n",
        "charge SYS\n",
        "addions SYS \"\"\" + str(pos_neut) + \"\\n\"\n",
        "\"\"\"addions SYS Cl- 0\n",
        "check SYS\n",
        "charge SYS\n",
        "savepdb SYS \"\"\" + str(pdb_nw) + \"\\n\"\n",
        "\"\"\"saveamberparm SYS \"\"\" + str(top_nw) + \" \" + str(crd_nw) + \"\\n\"\n",
        "\"\"\"solvatebox SYS \"\"\" + str(water_box) + \" \" + str(size_box) +  \"\"\" 0.7 \"\"\" + \"\\n\"\n",
        "\"\"\"addIonsRand SYS \"\"\" + str(pos_num) + \"\"\" Cl- \"\"\" + str(Cl_num) + \"\\n\"\n",
        "\"\"\"saveamberparm SYS \"\"\" + str(top) + \" \" + str(crd) + \"\\n\"\n",
        "\"\"\"savepdb SYS \"\"\" + str(pdb) + \"\\n\"\n",
        "\"\"\"quit\"\"\")\n",
        "f.close()\n",
        "\n",
        "!chmod 700 run_tleap.sh 2>&1 1>/dev/null\n",
        "# !bash run_tleap.sh 2>&1 1>/dev/null\n",
        "!bash run_tleap.sh 2>&1 1>/dev/null\n",
        "\n",
        "\n",
        "pdb_amber = os.path.exists(pdb)\n",
        "top_amber = os.path.exists(top)\n",
        "crd_amber = os.path.exists(crd)\n",
        "\n",
        "if pdb_amber == True and top_amber == True and crd_amber == True:\n",
        "  print(\"Successfully generated topology! :-)\")\n",
        "else:\n",
        "  print(\"ERROR: Check if your input file is in AMBER/GLYCAM format! \")\n",
        "!rm *.sh temp.txt >/dev/null 2>&1  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C8kKR7bpI86W"
      },
      "source": [
        "## Let's take a look on our simulation box:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vmQ27nZLssjv",
        "cellView": "form"
      },
      "source": [
        "#@title **Show 3D structure**\n",
        "import ipywidgets\n",
        "from ipywidgets import interact, fixed\n",
        "import warnings\n",
        "warnings.filterwarnings('ignore')\n",
        "\n",
        "def show_pdb(show_glycan=True,\n",
        "             show_sidechains=False,\n",
        "             show_mainchain=False,\n",
        "             show_box=True,\n",
        "             color=\"None\"):\n",
        "\n",
        "  def mainchain(p, color=\"white\", model=0):\n",
        "    BB = ['C','O','N','CA']\n",
        "    p.addStyle({\"model\":model,'atom':BB},\n",
        "                       {'stick':{'colorscheme':f\"{color}Carbon\",'radius':0.3}})\n",
        "    p.setViewStyle({'style':'outline','color':'black','width':0.1})\n",
        "\n",
        "  def glycan(p, model=0):\n",
        "    HP = ['0AA', '0AB', '0AD', '0AE', '0AF', '0AU', '0BA', '0BB', '0BC', '0BD', '0BU', '0CA', \n",
        "          '0CB', '0CD', '0CU', '0DA', '0DB', '0DD', '0DU', '0EA', '0EB', '0FA', '0FB', '0GA', \n",
        "          '0GB', '0GL', '0HA', '0HB', '0JA', '0JB', '0JD', '0JU', '0KA', '0KB', '0LA', '0LB', \n",
        "          '0MA', '0MB', '0NA', '0NB', '0OA', '0OB', '0PA', '0PB', '0PD', '0PU', '0QA', '0QB', \n",
        "          '0RA', '0RB', '0RD', '0RU', '0SA', '0SB', '0TA', '0TB', '0TV', '0Tv', '0UA', '0UB', \n",
        "          '0VA', '0VB', '0WA', '0WB', '0XA', '0XB', '0XD', '0XU', '0YA', '0YB', '0ZA', '0ZB', \n",
        "          '0aA', '0aB', '0aD', '0aU', '0bA', '0bB', '0bC', '0bD', '0bU', '0cA', '0cB', '0cD', \n",
        "          '0cU', '0dA', '0dB', '0dD', '0dR', '0dU', '0eA', '0eB', '0fA', '0fB', '0gA', '0gB', \n",
        "          '0gL', '0hA', '0hB', '0jA', '0jB', '0jD', '0jU', '0kA', '0kB', '0lA', '0lB', '0mA', \n",
        "          '0mB', '0nA', '0nB', '0oA', '0oB', '0pA', '0pB', '0pD', '0pU', '0qA', '0qB', '0rA', \n",
        "          '0rB', '0rD', '0rU', '0sA', '0sB', '0tA', '0tB', '0tV', '0tv', '0uA', '0uB', '0vA', \n",
        "          '0vB', '0wA', '0wB', '0xA', '0xB', '0xD', '0xU', '0yA', '0yB', '0zA', '0zB', '1AA', \n",
        "          '1AB', '1AD', '1AU', '1BA', '1BB', '1BD', '1BU', '1CA', '1CB', '1CD', '1CU', '1DA', \n",
        "          '1DB', '1DD', '1DU', '1EA', '1EB', '1FA', '1FB', '1GA', '1GB', '1HA', '1HB', '1JA', \n",
        "          '1JB', '1JD', '1JU', '1KA', '1KB', '1LA', '1LB', '1MA', '1MB', '1NA', '1NB', '1OA', \n",
        "          '1OB', '1PA', '1PB', '1PD', '1PU', '1QA', '1QB', '1RA', '1RB', '1RD', '1RU', '1TA', \n",
        "          '1TB', '1TV', '1Tv', '1UA', '1UB', '1VA', '1VB', '1WA', '1WB', '1XA', '1XB', '1XD', \n",
        "          '1XU', '1YA', '1YB', '1ZA', '1ZB', '1aA', '1aB', '1aD', '1aU', '1bA', '1bB', '1bD', \n",
        "          '1bU', '1cA', '1cB', '1cD', '1cU', '1dA', '1dB', '1dD', '1dU', '1eA', '1eB', '1fA', \n",
        "          '1fB', '1gA', '1gB', '1hA', '1hB', '1jA', '1jB', '1jD', '1jU', '1kA', '1kB', '1lA', \n",
        "          '1lB', '1mA', '1mB', '1nA', '1nB', '1oA', '1oB', '1pA', '1pB', '1pD', '1pU', '1qA', \n",
        "          '1qB', '1rA', '1rB', '1rD', '1rU', '1tA', '1tB', '1tV', '1tv', '1uA', '1uB', '1vA', \n",
        "          '1vB', '1wA', '1wB', '1xA', '1xB', '1xD', '1xU', '1yA', '1yB', '1zA', '1zB', '2AA', \n",
        "          '2AB', '2AD', '2AE', '2AF', '2AU', '2BA', '2BB', '2BD', '2BU', '2CA', '2CB', '2CD', \n",
        "          '2CU', '2DA', '2DB', '2DD', '2DU', '2EA', '2EB', '2FA', '2FB', '2GA', '2GB', '2HA', \n",
        "          '2HB', '2JA', '2JB', '2JD', '2JU', '2KA', '2KB', '2LA', '2LB', '2MA', '2MB', '2NA', \n",
        "          '2NB', '2OA', '2OB', '2PA', '2PB', '2PD', '2PU', '2QA', '2QB', '2RA', '2RB', '2RD', \n",
        "          '2RU', '2TA', '2TB', '2TV', '2Tv', '2UA', '2UB', '2XA', '2XB', '2XD', '2XU', '2ZA', \n",
        "          '2ZB', '2aA', '2aB', '2aD', '2aU', '2bA', '2bB', '2bD', '2bU', '2cA', '2cB', '2cD', \n",
        "          '2cU', '2dA', '2dB', '2dD', '2dU', '2eA', '2eB', '2fA', '2fB', '2gA', '2gB', '2hA', \n",
        "          '2hB', '2jA', '2jB', '2jD', '2jU', '2kA', '2kB', '2lA', '2lB', '2mA', '2mB', '2nA', \n",
        "          '2nB', '2oA', '2oB', '2pA', '2pB', '2pD', '2pU', '2qA', '2qB', '2rA', '2rB', '2rD', \n",
        "          '2rU', '2tA', '2tB', '2tV', '2tv', '2uA', '2uB', '2xA', '2xB', '2xD', '2xU', '2zA', \n",
        "          '2zB', '3AA', '3AB', '3AD', '3AU', '3BA', '3BB', '3BC', '3BD', '3BU', '3CA', '3CB', \n",
        "          '3CD', '3CU', '3DA', '3DB', '3DD', '3DU', '3EA', '3EB', '3FA', '3FB', '3GA', '3GB', \n",
        "          '3HA', '3HB', '3JA', '3JB', '3JD', '3JU', '3KA', '3KB', '3LA', '3LB', '3MA', '3MB', \n",
        "          '3NA', '3NB', '3OA', '3OB', '3PA', '3PB', '3PD', '3PU', '3QA', '3QB', '3RA', '3RB', \n",
        "          '3RD', '3RU', '3TA', '3TB', '3UA', '3UB', '3VA', '3VB', '3WA', '3WB', '3XA', '3XB', \n",
        "          '3XD', '3XU', '3YA', '3YB', '3ZA', '3ZB', '3aA', '3aB', '3aD', '3aU', '3bA', '3bB', \n",
        "          '3bC', '3bD', '3bU', '3cA', '3cB', '3cD', '3cU', '3dA', '3dB', '3dD', '3dR', '3dU', \n",
        "          '3eA', '3eB', '3fA', '3fB', '3gA', '3gB', '3hA', '3hB', '3jA', '3jB', '3jD', '3jU', \n",
        "          '3kA', '3kB', '3lA', '3lB', '3mA', '3mB', '3nA', '3nB', '3oA', '3oB', '3pA', '3pB', \n",
        "          '3pD', '3pU', '3qA', '3qB', '3rA', '3rB', '3rD', '3rU', '3tA', '3tB', '3uA', '3uB', \n",
        "          '3vA', '3vB', '3wA', '3wB', '3xA', '3xB', '3xD', '3xU', '3yA', '3yB', '3zA', '3zB', \n",
        "          '4AA', '4AB', '4AE', '4AF', '4BA', '4BB', '4BD', '4BU', '4CA', '4CB', '4CD', '4CU', \n",
        "          '4DA', '4DB', '4EA', '4EB', '4FA', '4FB', '4GA', '4GB', '4GL', '4HA', '4HB', '4JA', \n",
        "          '4JB', '4JD', '4JU', '4KA', '4KB', '4LA', '4LB', '4MA', '4MB', '4NA', '4NB', '4OA', \n",
        "          '4OB', '4PA', '4PB', '4PD', '4PU', '4QA', '4QB', '4RA', '4RB', '4SA', '4SB', '4TA', \n",
        "          '4TB', '4TV', '4Tv', '4UA', '4UB', '4VA', '4VB', '4WA', '4WB', '4XA', '4XB', '4YA', \n",
        "          '4YB', '4ZA', '4ZB', '4aA', '4aB', '4bA', '4bB', '4bD', '4bU', '4cA', '4cB', '4cD', \n",
        "          '4cU', '4dA', '4dB', '4dR', '4eA', '4eB', '4fA', '4fB', '4gA', '4gB', '4gL', '4hA', \n",
        "          '4hB', '4jA', '4jB', '4jD', '4jU', '4kA', '4kB', '4lA', '4lB', '4mA', '4mB', '4nA', \n",
        "          '4nB', '4oA', '4oB', '4pA', '4pB', '4pD', '4pU', '4qA', '4qB', '4rA', '4rB', '4sA', \n",
        "          '4sB', '4tA', '4tB', '4tV', '4tv', '4uA', '4uB', '4vA', '4vB', '4wA', '4wB', '4xA', \n",
        "          '4xB', '4yA', '4yB', '4zA', '4zB', '5AD', '5AU', '5BA', '5BB', '5CA', '5CB', '5DD', \n",
        "          '5DU', '5JA', '5JB', '5PA', '5PB', '5RD', '5RU', '5XD', '5XU', '5aD', '5aU', '5bA', \n",
        "          '5bB', '5cA', '5cB', '5dD', '5dU', '5jA', '5jB', '5pA', '5pB', '5rD', '5rU', '5xD', \n",
        "          '5xU', '6BD', '6BU', '6CD', '6CU', '6EA', '6EB', '6GA', '6GB', '6JD', '6JU', '6KA', \n",
        "          '6KB', '6LA', '6LB', '6MA', '6MB', '6NA', '6NB', '6PD', '6PU', '6TA', '6TB', '6VA', \n",
        "          '6VB', '6WA', '6WB', '6YA', '6YB', '6bD', '6bU', '6cD', '6cU', '6eA', '6eB', '6gA', \n",
        "          '6gB', '6jD', '6jU', '6kA', '6kB', '6lA', '6lB', '6mA', '6mB', '6nA', '6nB', '6pD', \n",
        "          '6pU', '6tA', '6tB', '6vA', '6vB', '6wA', '6wB', '6yA', '6yB', '7GL', '7SA', '7SB', \n",
        "          '7gL', '7sA', '7sB', '8GL', '8SA', '8SB', '8gL', '8sA', '8sB', '9GL', '9SA', '9SB', \n",
        "          '9gL', '9sA', '9sB', 'ACX', 'AGL', 'ASA', 'ASB', 'AgL', 'AsA', 'AsB', 'BGL', 'BSA', \n",
        "          'BSB', 'BgL', 'BsA', 'BsB', 'CA2', 'CGL', 'CSA', 'CSB', 'CgL', 'CsA', 'CsB', 'DGL', \n",
        "          'DSA', 'DSB', 'DgL', 'DsA', 'DsB', 'EGL', 'ESA', 'ESB', 'EgL', 'EsA', 'EsB', 'FGL', \n",
        "          'FSA', 'FSB', 'FgL', 'FsA', 'FsB', 'GGL', 'GSA', 'GSB', 'GgL', 'GsA', 'GsB', 'HGL', \n",
        "          'HSA', 'HSB', 'HgL', 'HsA', 'HsB', 'IGL', 'ISA', 'ISB', 'IgL', 'IsA', 'IsB', 'JGL', \n",
        "          'JSA', 'JSB', 'JgL', 'JsA', 'JsB', 'KGL', 'KSA', 'KSB', 'KgL', 'KsA', 'KsB', 'MEX', \n",
        "          'NLN', 'OLS', 'OLT', 'OME', 'PEA', 'PEB', 'PGA', 'PGB', 'PKA', 'PKB', 'PLA', 'PLB', \n",
        "          'PMA', 'PMB', 'PNA', 'PNB', 'PTA', 'PTB', 'PeA', 'PeB', 'PgA', 'PgB', 'PkA', 'PkB', \n",
        "          'PlA', 'PlB', 'PmA', 'PmB', 'PnA', 'PnB', 'PtA', 'PtB', 'QBD', 'QBU', 'QCD', 'QCU', \n",
        "          'QEA', 'QEB', 'QGA', 'QGB', 'QJD', 'QJU', 'QKA', 'QKB', 'QLA', 'QLB', 'QMA', 'QMB', \n",
        "          'QNA', 'QNB', 'QPD', 'QPU', 'QTA', 'QTB', 'QVA', 'QVB', 'QWA', 'QWB', 'QYA', 'QYB', \n",
        "          'QbD', 'QbU', 'QcD', 'QcU', 'QeA', 'QeB', 'QgA', 'QgB', 'QjD', 'QjU', 'QkA', 'QkB', \n",
        "          'QlA', 'QlB', 'QmA', 'QmB', 'QnA', 'QnB', 'QpD', 'QpU', 'QtA', 'QtB', 'QvA', 'QvB', \n",
        "          'QwA', 'QwB', 'QyA', 'QyB', 'REA', 'REB', 'RGA', 'RGB', 'RKA', 'RKB', 'RLA', 'RLB', \n",
        "          'RMA', 'RMB', 'RNA', 'RNB', 'ROH', 'RTA', 'RTB', 'ReA', 'ReB', 'RgA', 'RgB', 'RkA', \n",
        "          'RkB', 'RlA', 'RlB', 'RmA', 'RmB', 'RnA', 'RnB', 'RtA', 'RtB', 'SEA', 'SEB', 'SGA', \n",
        "          'SGB', 'SKA', 'SKB', 'SLA', 'SLB', 'SMA', 'SMB', 'SNA', 'SNB', 'SO3', 'STA', 'STB', \n",
        "          'SeA', 'SeB', 'SgA', 'SgB', 'SkA', 'SkB', 'SlA', 'SlB', 'SmA', 'SmB', 'SnA', 'SnB', \n",
        "          'StA', 'StB', 'TAA', 'TAB', 'TBT', 'TDA', 'TDB', 'TEA', 'TEB', 'TFA', 'TFB', 'TGA', \n",
        "          'TGB', 'THA', 'THB', 'TKA', 'TKB', 'TLA', 'TLB', 'TMA', 'TMB', 'TNA', 'TNB', 'TOA', \n",
        "          'TOB', 'TQA', 'TQB', 'TRA', 'TRB', 'TTA', 'TTB', 'TUA', 'TUB', 'TXA', 'TXB', 'TZA', \n",
        "          'TZB', 'TaA', 'TaB', 'TdA', 'TdB', 'TeA', 'TeB', 'TfA', 'TfB', 'TgA', 'TgB', 'ThA', \n",
        "          'ThB', 'TkA', 'TkB', 'TlA', 'TlB', 'TmA', 'TmB', 'TnA', 'TnB', 'ToA', 'ToB', 'TqA', \n",
        "          'TqB', 'TrA', 'TrB', 'TtA', 'TtB', 'TuA', 'TuB', 'TxA', 'TxB', 'TzA', 'TzB', 'UBD', \n",
        "          'UBU', 'UCD', 'UCU', 'UEA', 'UEB', 'UGA', 'UGB', 'UJD', 'UJU', 'UKA', 'UKB', 'ULA', \n",
        "          'ULB', 'UMA', 'UMB', 'UNA', 'UNB', 'UPD', 'UPU', 'UTA', 'UTB', 'UVA', 'UVB', 'UWA', \n",
        "          'UWB', 'UYA', 'UYB', 'UbD', 'UbU', 'UcD', 'UcU', 'UeA', 'UeB', 'UgA', 'UgB', 'UjD', \n",
        "          'UjU', 'UkA', 'UkB', 'UlA', 'UlB', 'UmA', 'UmB', 'UnA', 'UnB', 'UpD', 'UpU', 'UtA', \n",
        "          'UtB', 'UvA', 'UvB', 'UwA', 'UwB', 'UyA', 'UyB', 'VBD', 'VBU', 'VCD', 'VCU', 'VEA', \n",
        "          'VEB', 'VGA', 'VGB', 'VJD', 'VJU', 'VKA', 'VKB', 'VLA', 'VLB', 'VMA', 'VMB', 'VNA', \n",
        "          'VNB', 'VPD', 'VPU', 'VTA', 'VTB', 'VVA', 'VVB', 'VWA', 'VWB', 'VYA', 'VYB', 'VbD', \n",
        "          'VbU', 'VcD', 'VcU', 'VeA', 'VeB', 'VgA', 'VgB', 'VjD', 'VjU', 'VkA', 'VkB', 'VlA', \n",
        "          'VlB', 'VmA', 'VmB', 'VnA', 'VnB', 'VpD', 'VpU', 'VtA', 'VtB', 'VvA', 'VvB', 'VwA', \n",
        "          'VwB', 'VyA', 'VyB', 'WAA', 'WAB', 'WBA', 'WBB', 'WBD', 'WBU', 'WCA', 'WCB', 'WCD', \n",
        "          'WCU', 'WDA', 'WDB', 'WEA', 'WEB', 'WFA', 'WFB', 'WGA', 'WGB', 'WHA', 'WHB', 'WJA', \n",
        "          'WJB', 'WJD', 'WJU', 'WKA', 'WKB', 'WLA', 'WLB', 'WMA', 'WMB', 'WNA', 'WNB', 'WOA', \n",
        "          'WOB', 'WPA', 'WPB', 'WPD', 'WPU', 'WQA', 'WQB', 'WRA', 'WRB', 'WTA', 'WTB', 'WUA', \n",
        "          'WUB', 'WVA', 'WVB', 'WWA', 'WWB', 'WXA', 'WXB', 'WYA', 'WYB', 'WZA', 'WZB', 'WaA', \n",
        "          'WaB', 'WbA', 'WbB', 'WbD', 'WbU', 'WcA', 'WcB', 'WcD', 'WcU', 'WdA', 'WdB', 'WdR', \n",
        "          'WeA', 'WeB', 'WfA', 'WfB', 'WgA', 'WgB', 'WhA', 'WhB', 'WjA', 'WjB', 'WjD', 'WjU', \n",
        "          'WkA', 'WkB', 'WlA', 'WlB', 'WmA', 'WmB', 'WnA', 'WnB', 'WoA', 'WoB', 'WpA', 'WpB', \n",
        "          'WpD', 'WpU', 'WqA', 'WqB', 'WrA', 'WrB', 'WtA', 'WtB', 'WuA', 'WuB', 'WvA', 'WvB', \n",
        "          'WwA', 'WwB', 'WxA', 'WxB', 'WyA', 'WyB', 'WzA', 'WzB', 'XEA', 'XEB', 'XGA', 'XGB', \n",
        "          'XKA', 'XKB', 'XLA', 'XLB', 'XMA', 'XMB', 'XNA', 'XNB', 'XTA', 'XTB', 'XeA', 'XeB', \n",
        "          'XgA', 'XgB', 'XkA', 'XkB', 'XlA', 'XlB', 'XmA', 'XmB', 'XnA', 'XnB', 'XtA', 'XtB', \n",
        "          'YAA', 'YAB', 'YAF', 'YDA', 'YDB', 'YEA', 'YEB', 'YFA', 'YFB', 'YGA', 'YGB', 'YGa', \n",
        "          'YHA', 'YHB', 'YKA', 'YKB', 'YLA', 'YLB', 'YMA', 'YMB', 'YNA', 'YNB', 'YOA', 'YOB', \n",
        "          'YQA', 'YQB', 'YRA', 'YRB', 'YTA', 'YTB', 'YTV', 'YTv', 'YUA', 'YUB', 'YXA', 'YXB', \n",
        "          'YZA', 'YZB', 'YaA', 'YaB', 'YdA', 'YdB', 'YeA', 'YeB', 'YfA', 'YfB', 'YgA', 'YgB',\n",
        "          'YhA', 'YhB', 'YkA', 'YkB', 'YlA', 'YlB', 'YmA', 'YmB', 'YnA', 'YnB', 'YoA', 'YoB', \n",
        "          'YqA', 'YqB', 'YrA', 'YrB', 'YtA', 'YtB', 'YtV', 'Ytv', 'YuA', 'YuB', 'YxA', 'YxB', \n",
        "          'YzA', 'YzB', 'ZAA', 'ZAB', 'ZAD', 'ZAU', 'ZDA', 'ZDB', 'ZDD', 'ZDU', 'ZEA', 'ZEB', \n",
        "          'ZFA', 'ZFB', 'ZGA', 'ZGB', 'ZHA', 'ZHB', 'ZKA', 'ZKB', 'ZLA', 'ZLB', 'ZMA', 'ZMB', \n",
        "          'ZNA', 'ZNB', 'ZOA', 'ZOB', 'ZOL', 'ZOL', 'ZQA', 'ZQB', 'ZRA', 'ZRB', 'ZRD', 'ZRU', \n",
        "          'ZTA', 'ZTB', 'ZUA', 'ZUB', 'ZXA', 'ZXB', 'ZXD', 'ZXU', 'ZZA', 'ZZB', 'ZaA', 'ZaB', \n",
        "          'ZaD', 'ZaU', 'ZdA', 'ZdB', 'ZdD', 'ZdU', 'ZeA', 'ZeB', 'ZfA', 'ZfB', 'ZgA', 'ZgB', \n",
        "          'ZhA', 'ZhB', 'ZkA', 'ZkB', 'ZlA', 'ZlB', 'ZmA', 'ZmB', 'ZnA', 'ZnB', 'ZoA', 'ZoB', \n",
        "          'ZqA', 'ZqB', 'ZrA', 'ZrB', 'ZrD', 'ZrU', 'ZtA', 'ZtB', 'ZuA', 'ZuB', 'ZxA', 'ZxB', \n",
        "          'ZxD', 'ZxU', 'ZzA', 'ZzB', 'HYP', 'NLN', 'OLP', 'OLS', 'OLT', 'CHYP', 'CNLN', 'COLP', \n",
        "          'COLS', 'COLT', 'NHYP', 'NNLN', 'NOLP', 'NOLS', 'NOLT']\n",
        "    p.addStyle({\"model\":model,'and':[{'resn':HP}]},\n",
        "                       {'stick':{'colorscheme':'greenCarbon','radius':0.3}})\n",
        "    p.setViewStyle({'style':'outline','color':'black','width':0.1})\n",
        "\n",
        "\n",
        "  def box(p, model=0):\n",
        "      p.addModelsAsFrames(pdb)\n",
        "      p.addSurface(py3Dmol.SAS, {'opacity': 0.6, 'color':'white'}) #comment this line if you dont want to see the water box\n",
        "      p.setViewStyle({'style':'outline','color':'black','width':0.1})\n",
        "\n",
        "\n",
        "  def sidechain(p, model=0):\n",
        "    HP = [\"ALA\",\"GLY\",\"VAL\",\"ILE\",\"LEU\",\"PHE\",\"MET\",\"PRO\",\"TRP\",\"CYS\",\"TYR\"]\n",
        "    BB = ['C','O','N']\n",
        "    p.addStyle({\"model\":model,'and':[{'resn':HP},{'atom':BB,'invert':True}]},\n",
        "              {'stick':{'colorscheme':\"whiteCarbon\",'radius':0.3}})\n",
        "    p.addStyle({\"model\":model,'and':[{'resn':\"GLY\"},{'atom':'CA'}]},\n",
        "              {'sphere':{'colorscheme':\"whiteCarbon\",'radius':0.3}})\n",
        "    p.addStyle({\"model\":model,'and':[{'resn':\"PRO\"},{'atom':['C','O'],'invert':True}]},\n",
        "              {'stick':{'colorscheme':\"whiteCarbon\",'radius':0.3}})  \n",
        "    p.addStyle({\"model\":model,'and':[{'resn':HP,'invert':True},{'atom':BB,'invert':True}]},\n",
        "              {'stick':{'colorscheme':\"whiteCarbon\",'radius':0.3}})\n",
        "\n",
        "  p = py3Dmol.view(js='https://3dmol.org/build/3Dmol.js')\n",
        "  p.addModel(open(pdb,'r').read(),'pdb')\n",
        "  if color == \"rainbow\":\n",
        "    p.setStyle({'cartoon': {'color':'spectrum'}})\n",
        "  else:\n",
        "    p.setStyle({'cartoon':{}})\n",
        "\n",
        "  if show_sidechains: sidechain(p)\n",
        "  if show_mainchain: mainchain(p)\n",
        "  if show_glycan: glycan(p)\n",
        "  if show_box: box(p)\n",
        "  p.zoomTo()\n",
        "  return p.show()\n",
        "\n",
        "interact(show_pdb,\n",
        "         show_box=ipywidgets.Checkbox(value=True),\n",
        "         show_glycan=ipywidgets.Checkbox(value=True),\n",
        "         show_sidechains=ipywidgets.Checkbox(value=False),\n",
        "         show_mainchain=ipywidgets.Checkbox(value=False),\n",
        "         color=ipywidgets.Dropdown(options=['None', 'rainbow'], value='None'))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n85MrAO7M7uQ"
      },
      "source": [
        "---\n",
        "---\n",
        "# **Equilibrating the simulation box**\n",
        "\n",
        "Proper MD equilibration protocol is designed to equilibrate both temperature and pressure throughout the simulation box while preserving the protein experimental conformation. In addition, we also allow the solvent to accomodate around the protein, creating proper solvation layers.\n",
        "\n",
        "Below, we will set up the MD equilibration parameters, such as temperature, pressure and the desired simulation time. We will define the force constant used to restraint protein heavy-atoms in place and the frequency at which we want to save atomic coordinates in a trajectory file (.dcd).\n",
        "\n",
        "After you are done, you can run the next 2 cells to equilibrate your system."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8x9Qp_dbr9HP",
        "cellView": "form"
      },
      "source": [
        "#@title ### **Parameters for MD Equilibration protocol:**\n",
        "\n",
        "# remove whitespaces\n",
        "Jobname = 'equil' #@param {type:\"string\"}\n",
        "\n",
        "Minimization_steps = \"1000\" #@param [\"1000\", \"5000\", \"10000\", \"20000\", \"50000\", \"100000\"]\n",
        "\n",
        "#@markdown Simulation time (in nanoseconds) and integration time (in femtoseconds): \n",
        "Time = \"5\" #@param {type:\"string\"}\n",
        "stride_time_eq = Time\n",
        "Integration_timestep = \"2\" #@param [\"0.5\", \"1\", \"2\", \"3\", \"4\"]\n",
        "dt_eq = Integration_timestep\n",
        "\n",
        "#@markdown Temperature (in Kelvin) and Pressure (in bar)\n",
        "Temperature = 298 #@param {type:\"string\"}\n",
        "temperature_eq = Temperature\n",
        "Pressure = 1 #@param {type:\"string\"}\n",
        "pressure_eq = Pressure\n",
        "\n",
        "#@markdown Position restraints force constant (in kJ/mol): \n",
        "Force_constant = 500 #@param {type:\"slider\", min:0, max:2000, step:100}\n",
        "\n",
        "#@markdown Frequency to write the trajectory file (in picoseconds): \n",
        "\n",
        "Write_the_trajectory = \"10\" #@param [\"10\", \"100\", \"200\", \"500\", \"1000\"]\n",
        "write_the_trajectory_eq = Write_the_trajectory\n",
        "#@markdown Frequency to write the log file (in picoseconds): \n",
        "\n",
        "Write_the_log = \"10\" #@param [\"10\", \"100\", \"200\", \"500\", \"1000\"]\n",
        "write_the_log_eq = Write_the_log\n",
        "\n",
        "\n",
        "#@markdown ---\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zoamR9iynphz",
        "cellView": "form"
      },
      "source": [
        "#@title **Runs an Equilibration MD simulation (NPT ensemble)**\n",
        "#@markdown Now, let's equilibrate our system!\n",
        "\n",
        "###########################################\n",
        "import openmm as mm\n",
        "from openmm import *\n",
        "from openmm.app import *\n",
        "from openmm.unit import *\n",
        "import pytraj as pt\n",
        "\n",
        "from sys import stdout, exit, stderr\n",
        "import os, math, fnmatch\n",
        "\n",
        "#############################################\n",
        "# Defining MD simulation parameters\n",
        "\n",
        "jobname = os.path.join(workDir, Jobname)\n",
        "coordinatefile = os.path.join(workDir, \"SYS.crd\")\n",
        "pdbfile = os.path.join(workDir, \"SYS.pdb\")\n",
        "topologyfile = os.path.join(workDir, \"SYS.prmtop\")\n",
        "\n",
        "time_ps = float(Time)*1000\n",
        "simulation_time = float(time_ps)*picosecond\t\t# in ps\n",
        "dt = int(dt_eq)*femtosecond\t\t\t\t\t\n",
        "temperature = float(temperature_eq)*kelvin\n",
        "savcrd_freq = int(write_the_trajectory_eq)*picosecond\n",
        "print_freq  = int(write_the_log_eq)*picosecond\n",
        "\n",
        "pressure\t= float(pressure_eq)*bar\n",
        "\n",
        "restraint_fc = int(Force_constant) # kJ/mol\n",
        "\n",
        "nsteps  = int(simulation_time.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "nprint  = int(print_freq.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "nsavcrd = int(savcrd_freq.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "\n",
        "#############################################\n",
        "# Defining functions to use below:\n",
        "def backup_old_log(pattern, string):\n",
        "\tresult = []\n",
        "\tfor root, dirs, files in os.walk(\"./\"):\n",
        "\t\tfor name in files:\n",
        "\t\t\tif fnmatch.fnmatch(name, pattern):\n",
        "\n",
        "\t\t\t\ttry:\n",
        "\t\t\t\t\tnumber = int(name[-2])\n",
        "\t\t\t\t\tavail = isinstance(number, int)\n",
        "\t\t\t\t\t#print(name,avail)\n",
        "\t\t\t\t\tif avail == True:\n",
        "\t\t\t\t\t\tresult.append(number)\n",
        "\t\t\t\texcept:\n",
        "\t\t\t\t\tpass\n",
        "\n",
        "\tif len(result) > 0:\n",
        "\t\tmaxnumber = max(result)\n",
        "\telse:\n",
        "\t\tmaxnumber = 0\n",
        "\n",
        "\tbackup_file = \"\\#\" + string + \".\" + str(maxnumber + 1) + \"#\"\n",
        "\tos.system(\"mv \" + string + \" \" + backup_file)\n",
        "\treturn backup_file\n",
        "\n",
        "def restraints(system, crd, fc, restraint_array):\n",
        "\n",
        "\tboxlx = system.getDefaultPeriodicBoxVectors()[0][0].value_in_unit(nanometers)\n",
        "\tboxly = system.getDefaultPeriodicBoxVectors()[1][1].value_in_unit(nanometers)\n",
        "\tboxlz = system.getDefaultPeriodicBoxVectors()[2][2].value_in_unit(nanometers)\n",
        "\n",
        "\tif fc > 0:\n",
        "\t\t# positional restraints for all heavy-atoms\n",
        "\t\tposresPROT = CustomExternalForce('k*periodicdistance(x, y, z, x0, y0, z0)^2;')\n",
        "\t\tposresPROT.addPerParticleParameter('k')\n",
        "\t\tposresPROT.addPerParticleParameter('x0')\n",
        "\t\tposresPROT.addPerParticleParameter('y0')\n",
        "\t\tposresPROT.addPerParticleParameter('z0')\n",
        "  \n",
        "\t\tfor atom1 in restraint_array:\n",
        "\t\t\tatom1 = int(atom1)\n",
        "               \n",
        "\t\t\txpos  = crd.positions[atom1].value_in_unit(nanometers)[0]\n",
        "\t\t\typos  = crd.positions[atom1].value_in_unit(nanometers)[1]\n",
        "\t\t\tzpos  = crd.positions[atom1].value_in_unit(nanometers)[2]\n",
        "\n",
        "\t\t\tposresPROT.addParticle(atom1, [fc, xpos, ypos, zpos])\n",
        "    \n",
        "\t\tsystem.addForce(posresPROT)\n",
        "\n",
        "\treturn system\n",
        "##############################################\n",
        "\n",
        "#############################################\n",
        "print(\"\\n> Simulation details:\\n\")\n",
        "print(\"\\tJob name = \" + jobname)\n",
        "print(\"\\tCoordinate file = \" + str(coordinatefile))\n",
        "print(\"\\tPDB file = \" + str(pdbfile))\n",
        "print(\"\\tTopology file = \" + str(topologyfile))\n",
        "\n",
        "print(\"\\n\\tSimulation_time = \" + str(simulation_time))\n",
        "print(\"\\tIntegration timestep = \" + str(dt))\n",
        "print(\"\\tTotal number of steps = \" +  str(nsteps))\n",
        "\n",
        "print(\"\\n\\tSave coordinates each \" + str(savcrd_freq))\n",
        "print(\"\\tPrint in log file each \" + str(print_freq))\n",
        "\n",
        "print(\"\\n\\tTemperature = \" + str(temperature))\n",
        "print(\"\\tPressure = \" + str(pressure))\n",
        "#############################################\n",
        "\n",
        "print(\"\\n> Setting the system:\\n\")\n",
        "\n",
        "print(\"\\t- Reading topology and structure file...\")\n",
        "prmtop = AmberPrmtopFile(topologyfile)\n",
        "inpcrd = AmberInpcrdFile(coordinatefile)\n",
        "\n",
        "print(\"\\t- Creating system and setting parameters...\")\n",
        "nonbondedMethod = PME\n",
        "nonbondedCutoff = 1.0*nanometers\n",
        "ewaldErrorTolerance = 0.0005\n",
        "constraints = HBonds\n",
        "rigidWater = True\n",
        "constraintTolerance = 0.000001\n",
        "friction = 1.0\n",
        "system = prmtop.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,\n",
        "                           constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance)\n",
        "\n",
        "print(\"\\t- Applying restraints. Force Constant = \" + str(Force_constant) + \"kJ/mol\")\n",
        "pt_system = pt.iterload(coordinatefile, topologyfile)\n",
        "pt_topology = pt_system.top\n",
        "restraint_array = pt.select_atoms('!(:H*) & !(:WAT) & !(:Na+) & !(:Cl-) & !(:Mg+) & !(:K+)', pt_topology)\n",
        "\n",
        "system = restraints(system, inpcrd, restraint_fc, restraint_array)\n",
        "\n",
        "print(\"\\t- Setting barostat...\")\n",
        "system.addForce(MonteCarloBarostat(pressure, temperature))\n",
        "\n",
        "print(\"\\t- Setting integrator...\")\n",
        "integrator = LangevinIntegrator(temperature, friction, dt)\n",
        "integrator.setConstraintTolerance(constraintTolerance)\n",
        "simulation = Simulation(prmtop.topology, system, integrator)\n",
        "simulation.context.setPositions(inpcrd.positions)\n",
        "if inpcrd.boxVectors is not None:\n",
        "    simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)\n",
        "\n",
        "print(\"\\t- Energy minimization: \" + str(Minimization_steps) + \" steps\")\n",
        "simulation.minimizeEnergy(tolerance=10*kilojoule/mole, maxIterations=int(Minimization_steps))\n",
        "print(\"\\t-> Potential Energy = \" + str(simulation.context.getState(getEnergy=True).getPotentialEnergy()))\n",
        "\n",
        "print(\"\\t- Setting initial velocities...\")\n",
        "simulation.context.setVelocitiesToTemperature(temperature)\n",
        "\n",
        "#############################################\n",
        "# Running Equilibration on NPT ensemble\n",
        "\n",
        "dcd_file = jobname + \".dcd\"\n",
        "log_file = jobname + \".log\"\n",
        "rst_file = jobname + \".rst\"\n",
        "prv_rst_file = jobname + \".rst\"\n",
        "pdb_file = jobname + \".pdb\"\n",
        "\n",
        "# Creating a trajectory file and reporters\n",
        "dcd = DCDReporter(dcd_file, nsavcrd)\n",
        "firstdcdstep = (nsteps) + nsavcrd\n",
        "dcd._dcd = DCDFile(dcd._out, simulation.topology, simulation.integrator.getStepSize(), firstdcdstep, nsavcrd) # charmm doesn't like first step to be 0\n",
        "\n",
        "simulation.reporters.append(dcd)\n",
        "simulation.reporters.append(StateDataReporter(stdout, nprint, step=True, speed=True, progress=True, totalSteps=nsteps, remainingTime=True, separator='\\t\\t'))\n",
        "simulation.reporters.append(StateDataReporter(log_file, nprint, step=True, kineticEnergy=True, potentialEnergy=True, totalEnergy=True, temperature=True, volume=True, speed=True))\n",
        "\n",
        "print(\"\\n> Simulating \" + str(nsteps) + \" steps...\")\n",
        "simulation.step(nsteps)\n",
        "\n",
        "simulation.reporters.clear() # remove all reporters so the next iteration don't trigger them.\n",
        "\n",
        "\n",
        "##################################\n",
        "# Writing last frame information of stride\n",
        "print(\"\\n> Writing state file (\" + str(rst_file) + \")...\")\n",
        "state = simulation.context.getState( getPositions=True, getVelocities=True )\n",
        "with open(rst_file, 'w') as f:\n",
        "\tf.write(XmlSerializer.serialize(state))\n",
        "\n",
        "last_frame = int(nsteps/nsavcrd)\n",
        "print(\"> Writing coordinate file (\" + str(pdb_file) + \", frame = \" + str(last_frame) + \")...\")\n",
        "positions = simulation.context.getState(getPositions=True).getPositions()\n",
        "PDBFile.writeFile(simulation.topology, positions, open(pdb_file, 'w'))\n",
        "\n",
        "print(\"\\n> Finished!\\n\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LXyL26HCO8Bu"
      },
      "source": [
        "---\n",
        "---\n",
        "# **Running a Production MD simulation**\n",
        "\n",
        "Finally, we will proceed with the Production simulation itself using the equilibrated system coordinates as input structure.\n",
        "\n",
        "Note that we will use here a *.rst state file* , which contains atomic velocities and positions from the last frame of the equilibration simulation, guaranteeing that our production simulation begins from a thermodynamically equilibrated system.\n",
        "\n",
        "Another important information here is the **Number_of_strides** and the **Stride_Time**. In this notebook, we simulate a defined number of *strides*, so the **simulation time = Number_of_strides*Stride_Time**. For example, we can simulate 100ns by setting *Number_of_strides=10* and *Stride_Time=10 ns*.\n",
        "\n",
        "**Important: at the end of the Production simulation, we concatenate all strides to create a complete trajectory file which can be visualized and analyzed**\n",
        "\n",
        "The idea behind this approach is to make use of the intermitent 12h/24h period that Google Colab allows us to use its GPUs. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z0JV6Zid50_o",
        "cellView": "form"
      },
      "source": [
        "#@markdown ### **Provide input file names below:** \n",
        "\n",
        "Equilibrated_PDB = 'equil.pdb' #@param {type:\"string\"}\n",
        "State_file = 'equil.rst' #@param {type:\"string\"}\n",
        "\n",
        "#@markdown ---\n",
        "#@markdown ### **Parameters for MD Production protocol:**\n",
        "\n",
        "\n",
        "# remove whitespaces\n",
        "Jobname = 'prod' #@param {type:\"string\"}\n",
        "\n",
        "#@markdown Simulation time (in nanoseconds), number of strides (integers) and integration timestep (in femtoseconds): \n",
        "Stride_Time = \"10\" #@param {type:\"string\"}\n",
        "stride_time_prod = Stride_Time\n",
        "Number_of_strides = \"1\" #@param {type:\"string\"}\n",
        "nstride = Number_of_strides\n",
        "Integration_timestep = \"2\" #@param [\"0.5\", \"1\", \"2\", \"3\", \"4\"]\n",
        "dt_prod = Integration_timestep\n",
        "\n",
        "#@markdown Temperature (in Kelvin) and Pressure (in bar)\n",
        "Temperature = 298 #@param {type:\"string\"}\n",
        "temperature_prod = Temperature\n",
        "Pressure = 1 #@param {type:\"string\"}\n",
        "pressure_prod = Pressure\n",
        "\n",
        "#@markdown Frequency to write the trajectory file (in picoseconds): \n",
        "Write_the_trajectory = \"10\" #@param [\"10\", \"100\", \"200\", \"500\", \"1000\"]\n",
        "write_the_trajectory_prod = Write_the_trajectory\n",
        "\n",
        "#@markdown Frequency to write the log file (in picoseconds): \n",
        "Write_the_log = \"10\" #@param [\"10\", \"100\", \"200\", \"500\", \"1000\"]\n",
        "write_the_log_prod = Write_the_log\n",
        "\n",
        "#@markdown ---"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0QcjKSoqHHhi",
        "cellView": "form"
      },
      "source": [
        "#@title **Runs a Production MD simulation (NPT ensemble) after equilibration**\n",
        "#\n",
        "###########################################\n",
        "import openmm as mm\n",
        "from openmm import *\n",
        "from openmm.app import *\n",
        "from openmm.unit import *\n",
        "\n",
        "from sys import stdout, exit, stderr\n",
        "import os, math, fnmatch\n",
        "\n",
        "#############################################\n",
        "# Defining MD simulation parameters\n",
        "\n",
        "jobname = os.path.join(workDir, str(Jobname))\n",
        "coordinatefile = os.path.join(workDir, \"SYS.crd\")\n",
        "pdbfile = os.path.join(workDir, Equilibrated_PDB)\n",
        "topologyfile = os.path.join(workDir, \"SYS.prmtop\")\n",
        "equil_rst_file = os.path.join(workDir, State_file)\n",
        "\n",
        "\n",
        "stride_time_ps = float(stride_time_prod)*1000\n",
        "stride_time = float(stride_time_ps)*picosecond        \n",
        "nstride = int(Number_of_strides)\n",
        "dt = int(dt_prod)*femtosecond\t\t\t\t\t\n",
        "temperature = float(temperature_prod)*kelvin\n",
        "savcrd_freq = int(write_the_trajectory_prod)*picosecond\n",
        "print_freq  = int(write_the_log_prod)*picosecond\n",
        "\n",
        "pressure\t= float(pressure_prod)*bar\n",
        "\n",
        "simulation_time = stride_time*nstride\n",
        "nsteps  = int(stride_time.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "nprint  = int(print_freq.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "nsavcrd = int(savcrd_freq.value_in_unit(picosecond)/dt.value_in_unit(picosecond))\n",
        "firststride = 1 # must be integer\n",
        "#############################################\n",
        "# Defining functions to use below:\n",
        "def backup_old_log(pattern, string):\n",
        "\tresult = []\n",
        "\tfor root, dirs, files in os.walk(\"./\"):\n",
        "\t\tfor name in files:\n",
        "\t\t\tif fnmatch.fnmatch(name, pattern):\n",
        "\n",
        "\t\t\t\ttry:\n",
        "\t\t\t\t\tnumber = int(name[-2])\n",
        "\t\t\t\t\tavail = isinstance(number, int)\n",
        "\t\t\t\t\t#print(name,avail)\n",
        "\t\t\t\t\tif avail == True:\n",
        "\t\t\t\t\t\tresult.append(number)\n",
        "\t\t\t\texcept:\n",
        "\t\t\t\t\tpass\n",
        "\n",
        "\tif len(result) > 0:\n",
        "\t\tmaxnumber = max(result)\n",
        "\telse:\n",
        "\t\tmaxnumber = 0\n",
        "\n",
        "\tbackup_file = \"\\#\" + string + \".\" + str(maxnumber + 1) + \"#\"\n",
        "\tos.system(\"mv \" + string + \" \" + backup_file)\n",
        "\treturn backup_file\n",
        "##############################################\n",
        "\n",
        "#############################################\n",
        "print(\"\\n> Simulation details:\\n\")\n",
        "print(\"\\tJob name = \" + jobname)\n",
        "print(\"\\tCoordinate file = \" + str(coordinatefile))\n",
        "print(\"\\tPDB file = \" + str(pdbfile))\n",
        "print(\"\\tTopology file = \" + str(topologyfile))\n",
        "\n",
        "print(\"\\n\\tSimulation_time = \" + str(stride_time*nstride))\n",
        "print(\"\\tIntegration timestep = \" + str(dt))\n",
        "print(\"\\tTotal number of steps = \" +  str(nsteps*nstride))\n",
        "print(\"\\tNumber of strides = \" + str(nstride) + \" (\" + str(stride_time) + \" in each stride)\")\n",
        "\n",
        "print(\"\\n\\tSave coordinates each \" + str(savcrd_freq))\n",
        "print(\"\\tSave checkpoint each \" + str(savcrd_freq))\n",
        "print(\"\\tPrint in log file each \" + str(print_freq))\n",
        "\n",
        "print(\"\\n\\tTemperature = \" + str(temperature))\n",
        "print(\"\\tPressure = \" + str(pressure))\n",
        "#############################################\n",
        "\n",
        "print(\"\\n> Setting the system:\\n\")\n",
        "\n",
        "print(\"\\t- Reading topology and structure file...\")\n",
        "prmtop = AmberPrmtopFile(topologyfile)\n",
        "inpcrd = AmberInpcrdFile(coordinatefile)\n",
        "\n",
        "\n",
        "print(\"\\t- Creating system and setting parameters...\")\n",
        "nonbondedMethod = PME\n",
        "nonbondedCutoff = 1.0*nanometers\n",
        "ewaldErrorTolerance = 0.0005\n",
        "constraints = HBonds\n",
        "rigidWater = True\n",
        "constraintTolerance = 0.000001\n",
        "friction = 1.0\n",
        "system = prmtop.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,\n",
        "                           constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance)\n",
        "\n",
        "print(\"\\t- Setting barostat...\")\n",
        "system.addForce(MonteCarloBarostat(pressure, temperature))\n",
        "\n",
        "print(\"\\t- Setting integrator...\")\n",
        "integrator = LangevinIntegrator(temperature, friction, dt)\n",
        "integrator.setConstraintTolerance(constraintTolerance)\n",
        "simulation = Simulation(prmtop.topology, system, integrator)\n",
        "simulation.context.setPositions(inpcrd.positions)\n",
        "if inpcrd.boxVectors is not None:\n",
        "\tsimulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)\n",
        "\n",
        "#############################################\n",
        "# Opening a loop of extension NSTRIDE to simulate the entire STRIDE_TIME*NSTRIDE\n",
        "for n in range(1, nstride + 1):\n",
        "\n",
        "\tprint(\"\\n\\n>>> Simulating Stride #\" + str(n) + \" <<<\")\n",
        "\n",
        "\tdcd_file = jobname + \"_\" + str(n) + \".dcd\"\n",
        "\tlog_file = jobname + \"_\" + str(n) + \".log\"\n",
        "\trst_file = jobname + \"_\" + str(n) + \".rst\"\n",
        "\tprv_rst_file = jobname + \"_\" + str(n-1) + \".rst\"\n",
        "\tpdb_file = jobname + \"_\" + str(n) + \".pdb\"\n",
        "\n",
        "\tif os.path.exists(rst_file):\n",
        "\t\tprint(\"> Stride #\" + str(n) + \" finished (\" + rst_file + \" present). Moving to next stride... <\")\n",
        "\t\tcontinue\n",
        "\n",
        "\tif n == 1:\n",
        "\t\tprint(\"\\n> Loading previous state from equilibration > \" + equil_rst_file + \" <\")\n",
        "\t\twith open(equil_rst_file, 'r') as f:\n",
        "\t\t\tsimulation.context.setState(XmlSerializer.deserialize(f.read()))\n",
        "\t\t\tcurrstep = int((n-1)*nsteps)\n",
        "\t\t\tcurrtime = currstep*dt.in_units_of(picosecond)\n",
        "\t\t\tsimulation.currentStep = currstep\n",
        "\t\t\tsimulation.context.setTime(currtime)\n",
        "\t\t\tprint(\"> Current time: \" + str(currtime) + \" (Step = \" + str(currstep) + \")\")\n",
        "\n",
        "\telse:\n",
        "\t\tprint(\"> Loading previous state from > \" + prv_rst_file + \" <\")\n",
        "\t\twith open(prv_rst_file, 'r') as f:\n",
        "\t\t\tsimulation.context.setState(XmlSerializer.deserialize(f.read()))\n",
        "\t\t\tcurrstep = int((n-1)*nsteps)\n",
        "\t\t\tcurrtime = currstep*dt.in_units_of(picosecond)\n",
        "\t\t\tsimulation.currentStep = currstep\n",
        "\t\t\tsimulation.context.setTime(currtime)\n",
        "\t\t\tprint(\"> Current time: \" + str(currtime) + \" (Step = \" + str(currstep) + \")\")\n",
        "\n",
        "\n",
        "\tdcd = DCDReporter(dcd_file, nsavcrd)\n",
        "\tfirstdcdstep = (currstep) + nsavcrd\n",
        "\tdcd._dcd = DCDFile(dcd._out, simulation.topology, simulation.integrator.getStepSize(), firstdcdstep, nsavcrd) # first step should not be 0\n",
        "\n",
        "\tsimulation.reporters.append(dcd)\n",
        "\tsimulation.reporters.append(StateDataReporter(stdout, nprint, step=True, speed=True, progress=True, totalSteps=(nsteps*nstride), remainingTime=True, separator='\\t\\t'))\n",
        "\tsimulation.reporters.append(StateDataReporter(log_file, nprint, step=True, kineticEnergy=True, potentialEnergy=True, totalEnergy=True, temperature=True, volume=True, speed=True))\n",
        "\n",
        "\tprint(\"\\n> Simulating \" + str(nsteps) + \" steps... (Stride #\" + str(n) + \")\")\n",
        "\tsimulation.step(nsteps)\n",
        "\n",
        "\tsimulation.reporters.clear() # remove all reporters so the next iteration don't trigger them.\n",
        "\n",
        "\n",
        "\t##################################\n",
        "\t# Writing last frame information of stride\n",
        "\tprint(\"\\n> Writing state file (\" + str(rst_file) + \")...\")\n",
        "\tstate = simulation.context.getState( getPositions=True, getVelocities=True )\n",
        "\twith open(rst_file, 'w') as f:\n",
        "\t\tf.write(XmlSerializer.serialize(state))\n",
        "\n",
        "\tlast_frame = int(nsteps/nsavcrd)\n",
        "\tprint(\"> Writing coordinate file (\" + str(pdb_file) + \", frame = \" + str(last_frame) + \")...\")\n",
        "\tpositions = simulation.context.getState(getPositions=True).getPositions()\n",
        "\tPDBFile.writeFile(simulation.topology, positions, open(pdb_file, 'w'))\n",
        "\n",
        "print(\"\\n> Finished!\\n\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DaLoQlJGf57o",
        "cellView": "form"
      },
      "source": [
        "#@title **Concatenate and align the trajectory**\n",
        "Skip = \"1\" #@param [\"1\", \"2\", \"5\", \"10\", \"20\", \"50\"]\n",
        "stride_traj = Skip\n",
        "Output_format = \"dcd\" #@param [\"dcd\", \"pdb\", \"trr\", \"xtc\"]\n",
        "\n",
        "#@markdown **Attention:** A high number of frames can explode the memory on Colab. You should be fine with 5000 frames or less.\n",
        "\n",
        "simulation_time_analysis = stride_time_ps*nstride\n",
        "simulation_ns = float(Stride_Time)*int(Number_of_strides)\n",
        "number_frames = int(simulation_time_analysis)/int(Write_the_trajectory)\n",
        "number_frames_analysis = number_frames/int(stride_traj)\n",
        "\n",
        "\n",
        "traj_end = os.path.join(workDir, str(Jobname) + \"_all.dcd\")\n",
        "traj_end2 = os.path.join(workDir, str(Jobname) + \"_all.\" + str(Output_format))\n",
        "template =  os.path.join(workDir, str(Jobname) + '_%s.dcd')\n",
        "\n",
        "flist = [template % str(i) for i in range(1, nstride + 1)]\n",
        "#print(flist)\n",
        "\n",
        "trajlist = pt.load(flist, pdb, stride=stride_traj)\n",
        "traj_image = trajlist.iterframe(autoimage=True, rmsfit=0)\n",
        "traj_write = pt.write_traj(traj_end, traj_image, overwrite=True)\n",
        "traj_load = pt.load(traj_end, pdb)\n",
        "traj_align = pt.align(traj_load, mask=\"!(:WAT) & !(:Na+) & !(:Cl-) & !(:Mg+) & !(:K+)\", ref=0) \n",
        "traj_write = pt.write_traj(traj_end, traj_align, overwrite=True, options='dcd')\n",
        "traj_write = pt.write_traj(traj_end2, traj_align, overwrite=True, options=Output_format)\n",
        "traj_load = pt.load(traj_end, pdb)\n",
        "print(traj_load)\n",
        "\n",
        "traj_end_check = os.path.exists(traj_end2)\n",
        "\n",
        "if traj_end_check == True:\n",
        "  print(\"Trajectory concatenated successfully! :-)\")\n",
        "else:\n",
        "  print(\"ERROR: Check your inputs! \")\n",
        "\n",
        "#mask_selection\n",
        "\n",
        "mask_protein = \"@CA\"\n",
        "pt_topology = traj_load.top\n",
        "restraint_array1 = pt.select_atoms(':0AA,0AB,0AD,0AE,0AF,0AU,0BA,0BB,0BC,0BD,0BU,0CA,0CB,0CD,0CU,0DA,0DB,0DD,0DU,0EA,0EB,0FA,0FB,0GA,0GB,0GL,0HA,0HB,0JA,0JB,0JD,0JU,0KA,0KB,0LA,0LB,0MA,0MB,0NA,0NB,0OA,0OB,0PA,0PB,0PD,0PU,0QA,0QB,0RA,0RB,0RD,0RU,0SA,0SB,0TA,0TB,0TV,0Tv,0UA,0UB,0VA,0VB,0WA,0WB,0XA,0XB,0XD,0XU,0YA,0YB,0ZA,0ZB,0aA,0aB,0aD,0aU,0bA,0bB,0bC,0bD,0bU,0cA,0cB,0cD,0cU,0dA,0dB,0dD,0dR,0dU,0eA,0eB,0fA,0fB,0gA,0gB,0gL,0hA,0hB,0jA,0jB,0jD,0jU,0kA,0kB,0lA,0lB,0mA,0mB,0nA,0nB,0oA,0oB,0pA,0pB,0pD,0pU,0qA,0qB,0rA,0rB,0rD,0rU,0sA,0sB,0tA,0tB,0tV,0tv,0uA,0uB,0vA,0vB,0wA,0wB,0xA,0xB,0xD,0xU,0yA,0yB,0zA,0zB,1AA,1AB,1AD,1AU,1BA,1BB,1BD,1BU,1CA,1CB,1CD,1CU,1DA,1DB,1DD,1DU,1EA,1EB,1FA,1FB,1GA,1GB,1HA,1HB,1JA,1JB,1JD,1JU,1KA,1KB,1LA,1LB,1MA,1MB,1NA,1NB,1OA,1OB,1PA,1PB,1PD,1PU,1QA,1QB,1RA,1RB,1RD,1RU,1TA,1TB,1TV,1Tv,1UA,1UB,1VA,1VB,1WA,1WB,1XA,1XB,1XD,1XU,1YA,1YB,1ZA,1ZB,1aA,1aB,1aD,1aU,1bA,1bB,1bD,1bU,1cA,1cB,1cD,1cU,1dA,1dB,1dD,1dU,1eA,1eB,1fA,1fB,1gA,1gB,1hA,1hB,1jA,1jB,1jD,1jU,1kA,1kB,1lA,1lB,1mA,1mB,1nA,1nB,1oA,1oB,1pA,1pB,1pD,1pU,1qA,1qB,1rA,1rB,1rD,1rU,1tA,1tB,1tV,1tv,1uA,1uB,1vA,1vB,1wA,1wB,1xA,1xB,1xD,1xU,1yA,1yB,1zA,1zB,2AA,2AB,2AD,2AE,2AF,2AU,2BA,2BB,2BD,2BU,2CA,2CB,2CD,2CU,2DA,2DB,2DD,2DU,2EA,2EB,2FA,2FB,2GA,2GB,2HA,2HB,2JA,2JB,2JD,2JU,2KA,2KB,2LA,2LB,2MA,2MB,2NA,2NB,2OA,2OB,2PA,2PB,2PD,2PU,2QA,2QB,2RA,2RB,2RD,2RU,2TA,2TB,2TV,2Tv,2UA,2UB,2XA,2XB,2XD,2XU,2ZA,2ZB,2aA,2aB,2aD,2aU,2bA,2bB,2bD,2bU,2cA,2cB,2cD,2cU,2dA,2dB,2dD,2dU,2eA,2eB,2fA,2fB,2gA,2gB,2hA,2hB,2jA,2jB,2jD,2jU,2kA,2kB,2lA,2lB,2mA,2mB,2nA,2nB,2oA,2oB,2pA,2pB,2pD,2pU,2qA,2qB,2rA,2rB,2rD,2rU,2tA,2tB,2tV,2tv,2uA,2uB,2xA,2xB,2xD,2xU,2zA,2zB,3AA,3AB,3AD,3AU,3BA,3BB,3BC,3BD,3BU,3CA,3CB,3CD,3CU,3DA,3DB,3DD,3DU,3EA,3EB,3FA,3FB,3GA,3GB,3HA,3HB,3JA,3JB,3JD,3JU,3KA,3KB,3LA,3LB,3MA,3MB,3NA,3NB,3OA,3OB,3PA,3PB,3PD,3PU,3QA,3QB,3RA,3RB,3RD,3RU,3TA,3TB,3UA,3UB,3VA,3VB,3WA,3WB,3XA,3XB,3XD,3XU,3YA,3YB,3ZA,3ZB,3aA,3aB,3aD,3aU,3bA,3bB,3bC,3bD,3bU,3cA,3cB,3cD,3cU,3dA,3dB,3dD,3dR,3dU,3eA,3eB,3fA,3fB,3gA,3gB,3hA,3hB,3jA,3jB,3jD,3jU,3kA,3kB,3lA,3lB,3mA,3mB,3nA,3nB,3oA,3oB,3pA,3pB,3pD,3pU,3qA,3qB,3rA,3rB,3rD,3rU,3tA,3tB,3uA,3uB,3vA,3vB,3wA,3wB,3xA,3xB,3xD,3xU,3yA,3yB,3zA,3zB,4AA,4AB,4AE,4AF,4BA,4BB,4BD,4BU,4CA,4CB,4CD,4CU,4DA,4DB,4EA,4EB,4FA,4FB,4GA,4GB,4GL,4HA,4HB,4JA,4JB,4JD,4JU,4KA,4KB,4LA,4LB,4MA,4MB,4NA,4NB,4OA,4OB,4PA,4PB,4PD,4PU,4QA,4QB,4RA,4RB,4SA,4SB,4TA,4TB,4TV,4Tv,4UA,4UB,4VA,4VB,4WA,4WB,4XA,4XB,4YA,4YB,4ZA,4ZB,4aA,4aB,4bA,4bB,4bD,4bU,4cA,4cB,4cD,4cU,4dA,4dB,4dR,4eA,4eB,4fA,4fB,4gA,4gB,4gL,4hA,4hB,4jA,4jB,4jD,4jU,4kA,4kB,4lA,4lB,4mA,4mB,4nA,4nB,4oA,4oB,4pA,4pB,4pD,4pU,4qA,4qB,4rA,4rB,4sA,4sB,4tA,4tB,4tV,4tv,4uA,4uB,4vA,4vB,4wA,4wB,4xA,4xB,4yA,4yB,4zA,4zB,5AD,5AU,5BA,5BB,5CA,5CB,5DD,5DU,5JA,5JB,5PA,5PB,5RD,5RU,5XD,5XU,5aD,5aU,5bA,5bB,5cA,5cB,5dD,5dU,5jA,5jB,5pA,5pB,5rD,5rU,5xD,5xU,6BD,6BU,6CD,6CU,6EA,6EB,6GA,6GB,6JD,6JU,6KA,6KB,6LA,6LB,6MA,6MB,6NA,6NB,6PD,6PU,6TA,6TB,6VA,6VB,6WA,6WB,6YA,6YB,6bD,6bU,6cD,6cU,6eA,6eB,6gA,6gB,6jD,6jU,6kA,6kB,6lA,6lB,6mA,6mB,6nA,6nB,6pD,6pU,6tA,6tB,6vA,6vB,6wA,6wB,6yA,6yB,7GL,7SA,7SB,7gL,7sA,7sB,8GL,8SA,8SB,8gL,8sA,8sB,9GL,9SA,9SB,9gL,9sA,9sB,ACX,AGL,ASA,ASB,AgL,AsA,AsB,BGL,BSA,BSB,BgL,BsA,BsB,CA2,CGL,CSA,CSB,CgL,CsA,CsB,DGL,DSA,DSB,DgL,DsA,DsB,EGL,ESA,ESB,EgL,EsA,EsB,FGL,FSA,FSB,FgL,FsA,FsB,GGL,GSA,GSB,GgL,GsA,GsB,HGL,HSA,HSB,HgL,HsA,HsB,IGL,ISA,ISB,IgL,IsA,IsB,JGL,JSA,JSB,JgL,JsA,JsB,KGL,KSA,KSB,KgL,KsA,KsB,MEX,NLN,OLS,OLT,OME,PEA,PEB,PGA,PGB,PKA,PKB,PLA,PLB,PMA,PMB,PNA,PNB,PTA,PTB,PeA,PeB,PgA,PgB,PkA,PkB,PlA,PlB,PmA,PmB,PnA,PnB,PtA,PtB,QBD,QBU,QCD,QCU,QEA,QEB,QGA,QGB,QJD,QJU,QKA,QKB,QLA,QLB,QMA,QMB,QNA,QNB,QPD,QPU,QTA,QTB,QVA,QVB,QWA,QWB,QYA,QYB,QbD,QbU,QcD,QcU,QeA,QeB,QgA,QgB,QjD,QjU,QkA,QkB,QlA,QlB,QmA,QmB,QnA,QnB,QpD,QpU,QtA,QtB,QvA,QvB,QwA,QwB,QyA,QyB,REA,REB,RGA,RGB,RKA,RKB,RLA,RLB,RMA,RMB,RNA,RNB,ROH,RTA,RTB,ReA,ReB,RgA,RgB,RkA,RkB,RlA,RlB,RmA,RmB,RnA,RnB,RtA,RtB,SEA,SEB,SGA,SGB,SKA,SKB,SLA,SLB,SMA,SMB,SNA,SNB,SO3,STA,STB,SeA,SeB,SgA,SgB,SkA,SkB,SlA,SlB,SmA,SmB,SnA,SnB,StA,StB,TAA,TAB,TBT,TDA,TDB,TEA,TEB,TFA,TFB,TGA,TGB,THA,THB,TKA,TKB,TLA,TLB,TMA,TMB,TNA,TNB,TOA,TOB,TQA,TQB,TRA,TRB,TTA,TTB,TUA,TUB,TXA,TXB,TZA,TZB,TaA,TaB,TdA,TdB,TeA,TeB,TfA,TfB,TgA,TgB,ThA,ThB,TkA,TkB,TlA,TlB,TmA,TmB,TnA,TnB,ToA,ToB,TqA,TqB,TrA,TrB,TtA,TtB,TuA,TuB,TxA,TxB,TzA,TzB,UBD,UBU,UCD,UCU,UEA,UEB,UGA,UGB,UJD,UJU,UKA,UKB,ULA,ULB,UMA,UMB,UNA,UNB,UPD,UPU,UTA,UTB,UVA,UVB,UWA,UWB,UYA,UYB,UbD,UbU,UcD,UcU,UeA,UeB,UgA,UgB,UjD,UjU,UkA,UkB,UlA,UlB,UmA,UmB,UnA,UnB,UpD,UpU,UtA,UtB,UvA,UvB,UwA,UwB,UyA,UyB,VBD,VBU,VCD,VCU,VEA,VEB,VGA,VGB,VJD,VJU,VKA,VKB,VLA,VLB,VMA,VMB,VNA,VNB,VPD,VPU,VTA,VTB,VVA,VVB,VWA,VWB,VYA,VYB,VbD,VbU,VcD,VcU,VeA,VeB,VgA,VgB,VjD,VjU,VkA,VkB,VlA,VlB,VmA,VmB,VnA,VnB,VpD,VpU,VtA,VtB,VvA,VvB,VwA,VwB,VyA,VyB,WAA,WAB,WBA,WBB,WBD,WBU,WCA,WCB,WCD,WCU,WDA,WDB,WEA,WEB,WFA,WFB,WGA,WGB,WHA,WHB,WJA,WJB,WJD,WJU,WKA,WKB,WLA,WLB,WMA,WMB,WNA,WNB,WOA,WOB,WPA,WPB,WPD,WPU,WQA,WQB,WRA,WRB,WTA,WTB,WUA,WUB,WVA,WVB,WWA,WWB,WXA,WXB,WYA,WYB,WZA,WZB,WaA,WaB,WbA,WbB,WbD,WbU,WcA,WcB,WcD,WcU,WdA,WdB,WdR,WeA,WeB,WfA,WfB,WgA,WgB,WhA,WhB,WjA,WjB,WjD,WjU,WkA,WkB,WlA,WlB,WmA,WmB,WnA,WnB,WoA,WoB,WpA,WpB,WpD,WpU,WqA,WqB,WrA,WrB,WtA,WtB,WuA,WuB,WvA,WvB,WwA,WwB,WxA,WxB,WyA,WyB,WzA,WzB,XEA,XEB,XGA,XGB,XKA,XKB,XLA,XLB,XMA,XMB,XNA,XNB,XTA,XTB,XeA,XeB,XgA,XgB,XkA,XkB,XlA,XlB,XmA,XmB,XnA,XnB,XtA,XtB,YAA,YAB,YAF,YDA,YDB,YEA,YEB,YFA,YFB,YGA,YGB,YGa,YHA,YHB,YKA,YKB,YLA,YLB,YMA,YMB,YNA,YNB,YOA,YOB,YQA,YQB,YRA,YRB,YTA,YTB,YTV,YTv,YUA,YUB,YXA,YXB,YZA,YZB,YaA,YaB,YdA,YdB,YeA,YeB,YfA,YfB,YgA,YgB,YhA,YhB,YkA,YkB,YlA,YlB,YmA,YmB,YnA,YnB,YoA,YoB,YqA,YqB,YrA,YrB,YtA,YtB,YtV,Ytv,YuA,YuB,YxA,YxB,YzA,YzB,ZAA,ZAB,ZAD,ZAU,ZDA,ZDB,ZDD,ZDU,ZEA,ZEB,ZFA,ZFB,ZGA,ZGB,ZHA,ZHB,ZKA,ZKB,ZLA,ZLB,ZMA,ZMB,ZNA,ZNB,ZOA,ZOB,ZOL,ZOL,ZQA,ZQB,ZRA,ZRB,ZRD,ZRU,ZTA,ZTB,ZUA,ZUB,ZXA,ZXB,ZXD,ZXU,ZZA,ZZB,ZaA,ZaB,ZaD,ZaU,ZdA,ZdB,ZdD,ZdU,ZeA,ZeB,ZfA,ZfB,ZgA,ZgB,ZhA,ZhB,ZkA,ZkB,ZlA,ZlB,ZmA,ZmB,ZnA,ZnB,ZoA,ZoB,ZqA,ZqB,ZrA,ZrB,ZrD,ZrU,ZtA,ZtB,ZuA,ZuB,ZxA,ZxB,ZxD,ZxU,ZzA,ZzB,HYP,NLN,OLP,OLS,OLT,CHYP,CNLN,COLP,COLS,COLT,NHYP,NNLN,NOLP,NOLS,NOLT', pt_topology)\n",
        "first_atom1 = restraint_array1[0]\n",
        "last_atom1 = restraint_array1[-1]\n",
        "mask_glycan = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "mask_protein_glycan = mask_protein + \",\" + mask_glycan\n",
        "# print(mask_protein)\n",
        "# print(mask_glycan)\n",
        "# print(mask_protein_glycan)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x_awbwdZnp6Q",
        "cellView": "form"
      },
      "source": [
        "#@title **Load, view and check the trajectory**\n",
        "#@markdown This will take a few minutes. Another coffee would be great. :-)\n",
        "\n",
        "#@markdown If you have a huge system, this step can explode the memory on Colab. Please, be careful. \n",
        "\n",
        "import warnings\n",
        "warnings.filterwarnings('ignore')\n",
        "!rm *.pdb 2> /dev/null\n",
        "\n",
        "#py3dmol functions\n",
        "class Atom(dict):\n",
        "  def __init__(self, line):\n",
        "    self[\"type\"] = line[0:6].strip()\n",
        "    self[\"idx\"] = line[6:11].strip()\n",
        "    self[\"name\"] = line[12:16].strip()\n",
        "    self[\"resname\"] = line[17:20].strip()\n",
        "    self[\"resid\"] = int(int(line[22:26]))\n",
        "    self[\"x\"] = float(line[30:38])\n",
        "    self[\"y\"] = float(line[38:46])\n",
        "    self[\"z\"] = float(line[46:54])\n",
        "    self[\"sym\"] = line[76:78].strip()\n",
        "\n",
        "  def __str__(self):\n",
        "    line = list(\" \" * 80)\n",
        "    line[0:6] = self[\"type\"].ljust(6)\n",
        "    line[6:11] = self[\"idx\"].ljust(5)\n",
        "    line[12:16] = self[\"name\"].ljust(4)\n",
        "    line[17:20] = self[\"resname\"].ljust(3)\n",
        "    line[22:26] = str(self[\"resid\"]).ljust(4)\n",
        "    line[30:38] = str(self[\"x\"]).rjust(8)\n",
        "    line[38:46] = str(self[\"y\"]).rjust(8)\n",
        "    line[46:54] = str(self[\"z\"]).rjust(8)\n",
        "    line[76:78] = self[\"sym\"].rjust(2)\n",
        "    return \"\".join(line) + \"\\n\"\n",
        "        \n",
        "class Molecule(list):\n",
        "  def __init__(self, file):\n",
        "    for line in file:\n",
        "      if \"ATOM\" in line or \"HETATM\" in line:\n",
        "        self.append(Atom(line))\n",
        "            \n",
        "    def __str__(self):\n",
        "      outstr = \"\"\n",
        "      for at in self:\n",
        "        outstr += str(at)\n",
        "      return outstr\n",
        "\n",
        "if number_frames_analysis > 10:\n",
        "  stride_animation = number_frames_analysis/10\n",
        "else:\n",
        "  stride_animation = 1\n",
        "\n",
        "u = mda.Universe(pdb, traj_end) \n",
        "\n",
        "# Write out frames for animation\n",
        "protein = u.select_atoms('not (resname WAT)')\n",
        "i = 0\n",
        "for ts in u.trajectory[0:len(u.trajectory):int(stride_animation)]: \n",
        "    if i > -1:\n",
        "        with mda.Writer('' + str(i) + '.pdb', protein.n_atoms) as W:\n",
        "            W.write(protein)\n",
        "    i = i + 1\n",
        "# Load frames as molecules\n",
        "molecules = []\n",
        "for i in range(int(len(u.trajectory)/int(stride_animation))):\n",
        "    with open('' + str(i) + '.pdb') as ifile:\n",
        "        molecules.append(Molecule(ifile))\n",
        "\n",
        "models = \"\"\n",
        "for i in range(len(molecules)):\n",
        "  models += \"MODEL \" + str(i) + \"\\n\"\n",
        "  for j,mol in enumerate(molecules[i]):\n",
        "    models += str(mol)\n",
        "  models += \"ENDMDL\\n\"\n",
        "#view.addModelsAsFrames(models)\n",
        "\n",
        "# Animation\n",
        "view = py3Dmol.view(width=800, height=600)\n",
        "view.addModelsAsFrames(models)\n",
        "for i, at in enumerate(molecules[0]):\n",
        "    default = {\"cartoon\": {'color': 'spectrum'}}\n",
        "    view.setViewStyle({'style':'outline','color':'black','width':0.1})\n",
        "    view.setStyle({'model': -1, 'serial': i+1}, at.get(\"pymol\", default))\n",
        "    view.setStyle({'model': -1, 'serial': i+1}, {'stick':{'radius':0.3}})\n",
        "    # HP = ['protein']\n",
        "    # view.setStyle({\"model\":-1,'and':[{'resn':HP}]},{'stick':{'radius':0.3}})\n",
        "view.zoomTo()\n",
        "view.animate({'loop': \"forward\"})\n",
        "view.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Emh0vU5UjgB6"
      },
      "source": [
        "---\n",
        "---\n",
        "# **Analysis**\n",
        "\n",
        "Although visualizing your trajectory can be quite useful, sometimes you also want more quantitative data.\n",
        "\n",
        "Analyses of MD trajectories vary a lot and we do not intend to cover it all here. However, one can make use of MDanalysis or PyTraj to easily analyze simulations. \n",
        "\n",
        "Below, you can find a few examples of code snippets that can help you to shed some light on your simulation behavior."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Compute hydrogen bonds between your molecule and waters**\n",
        "\n",
        "Selection = \"Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  pt_topology = traj_load.top\n",
        "  restraint_array1 = pt.select_atoms('!(:WAT) & !(:Na+) & !(:Cl-) & !(:Mg+) & !(:K+) & !:0AA,0AB,0AD,0AE,0AF,0AU,0BA,0BB,0BC,0BD,0BU,0CA,0CB,0CD,0CU,0DA,0DB,0DD,0DU,0EA,0EB,0FA,0FB,0GA,0GB,0GL,0HA,0HB,0JA,0JB,0JD,0JU,0KA,0KB,0LA,0LB,0MA,0MB,0NA,0NB,0OA,0OB,0PA,0PB,0PD,0PU,0QA,0QB,0RA,0RB,0RD,0RU,0SA,0SB,0TA,0TB,0TV,0Tv,0UA,0UB,0VA,0VB,0WA,0WB,0XA,0XB,0XD,0XU,0YA,0YB,0ZA,0ZB,0aA,0aB,0aD,0aU,0bA,0bB,0bC,0bD,0bU,0cA,0cB,0cD,0cU,0dA,0dB,0dD,0dR,0dU,0eA,0eB,0fA,0fB,0gA,0gB,0gL,0hA,0hB,0jA,0jB,0jD,0jU,0kA,0kB,0lA,0lB,0mA,0mB,0nA,0nB,0oA,0oB,0pA,0pB,0pD,0pU,0qA,0qB,0rA,0rB,0rD,0rU,0sA,0sB,0tA,0tB,0tV,0tv,0uA,0uB,0vA,0vB,0wA,0wB,0xA,0xB,0xD,0xU,0yA,0yB,0zA,0zB,1AA,1AB,1AD,1AU,1BA,1BB,1BD,1BU,1CA,1CB,1CD,1CU,1DA,1DB,1DD,1DU,1EA,1EB,1FA,1FB,1GA,1GB,1HA,1HB,1JA,1JB,1JD,1JU,1KA,1KB,1LA,1LB,1MA,1MB,1NA,1NB,1OA,1OB,1PA,1PB,1PD,1PU,1QA,1QB,1RA,1RB,1RD,1RU,1TA,1TB,1TV,1Tv,1UA,1UB,1VA,1VB,1WA,1WB,1XA,1XB,1XD,1XU,1YA,1YB,1ZA,1ZB,1aA,1aB,1aD,1aU,1bA,1bB,1bD,1bU,1cA,1cB,1cD,1cU,1dA,1dB,1dD,1dU,1eA,1eB,1fA,1fB,1gA,1gB,1hA,1hB,1jA,1jB,1jD,1jU,1kA,1kB,1lA,1lB,1mA,1mB,1nA,1nB,1oA,1oB,1pA,1pB,1pD,1pU,1qA,1qB,1rA,1rB,1rD,1rU,1tA,1tB,1tV,1tv,1uA,1uB,1vA,1vB,1wA,1wB,1xA,1xB,1xD,1xU,1yA,1yB,1zA,1zB,2AA,2AB,2AD,2AE,2AF,2AU,2BA,2BB,2BD,2BU,2CA,2CB,2CD,2CU,2DA,2DB,2DD,2DU,2EA,2EB,2FA,2FB,2GA,2GB,2HA,2HB,2JA,2JB,2JD,2JU,2KA,2KB,2LA,2LB,2MA,2MB,2NA,2NB,2OA,2OB,2PA,2PB,2PD,2PU,2QA,2QB,2RA,2RB,2RD,2RU,2TA,2TB,2TV,2Tv,2UA,2UB,2XA,2XB,2XD,2XU,2ZA,2ZB,2aA,2aB,2aD,2aU,2bA,2bB,2bD,2bU,2cA,2cB,2cD,2cU,2dA,2dB,2dD,2dU,2eA,2eB,2fA,2fB,2gA,2gB,2hA,2hB,2jA,2jB,2jD,2jU,2kA,2kB,2lA,2lB,2mA,2mB,2nA,2nB,2oA,2oB,2pA,2pB,2pD,2pU,2qA,2qB,2rA,2rB,2rD,2rU,2tA,2tB,2tV,2tv,2uA,2uB,2xA,2xB,2xD,2xU,2zA,2zB,3AA,3AB,3AD,3AU,3BA,3BB,3BC,3BD,3BU,3CA,3CB,3CD,3CU,3DA,3DB,3DD,3DU,3EA,3EB,3FA,3FB,3GA,3GB,3HA,3HB,3JA,3JB,3JD,3JU,3KA,3KB,3LA,3LB,3MA,3MB,3NA,3NB,3OA,3OB,3PA,3PB,3PD,3PU,3QA,3QB,3RA,3RB,3RD,3RU,3TA,3TB,3UA,3UB,3VA,3VB,3WA,3WB,3XA,3XB,3XD,3XU,3YA,3YB,3ZA,3ZB,3aA,3aB,3aD,3aU,3bA,3bB,3bC,3bD,3bU,3cA,3cB,3cD,3cU,3dA,3dB,3dD,3dR,3dU,3eA,3eB,3fA,3fB,3gA,3gB,3hA,3hB,3jA,3jB,3jD,3jU,3kA,3kB,3lA,3lB,3mA,3mB,3nA,3nB,3oA,3oB,3pA,3pB,3pD,3pU,3qA,3qB,3rA,3rB,3rD,3rU,3tA,3tB,3uA,3uB,3vA,3vB,3wA,3wB,3xA,3xB,3xD,3xU,3yA,3yB,3zA,3zB,4AA,4AB,4AE,4AF,4BA,4BB,4BD,4BU,4CA,4CB,4CD,4CU,4DA,4DB,4EA,4EB,4FA,4FB,4GA,4GB,4GL,4HA,4HB,4JA,4JB,4JD,4JU,4KA,4KB,4LA,4LB,4MA,4MB,4NA,4NB,4OA,4OB,4PA,4PB,4PD,4PU,4QA,4QB,4RA,4RB,4SA,4SB,4TA,4TB,4TV,4Tv,4UA,4UB,4VA,4VB,4WA,4WB,4XA,4XB,4YA,4YB,4ZA,4ZB,4aA,4aB,4bA,4bB,4bD,4bU,4cA,4cB,4cD,4cU,4dA,4dB,4dR,4eA,4eB,4fA,4fB,4gA,4gB,4gL,4hA,4hB,4jA,4jB,4jD,4jU,4kA,4kB,4lA,4lB,4mA,4mB,4nA,4nB,4oA,4oB,4pA,4pB,4pD,4pU,4qA,4qB,4rA,4rB,4sA,4sB,4tA,4tB,4tV,4tv,4uA,4uB,4vA,4vB,4wA,4wB,4xA,4xB,4yA,4yB,4zA,4zB,5AD,5AU,5BA,5BB,5CA,5CB,5DD,5DU,5JA,5JB,5PA,5PB,5RD,5RU,5XD,5XU,5aD,5aU,5bA,5bB,5cA,5cB,5dD,5dU,5jA,5jB,5pA,5pB,5rD,5rU,5xD,5xU,6BD,6BU,6CD,6CU,6EA,6EB,6GA,6GB,6JD,6JU,6KA,6KB,6LA,6LB,6MA,6MB,6NA,6NB,6PD,6PU,6TA,6TB,6VA,6VB,6WA,6WB,6YA,6YB,6bD,6bU,6cD,6cU,6eA,6eB,6gA,6gB,6jD,6jU,6kA,6kB,6lA,6lB,6mA,6mB,6nA,6nB,6pD,6pU,6tA,6tB,6vA,6vB,6wA,6wB,6yA,6yB,7GL,7SA,7SB,7gL,7sA,7sB,8GL,8SA,8SB,8gL,8sA,8sB,9GL,9SA,9SB,9gL,9sA,9sB,ACX,AGL,ASA,ASB,AgL,AsA,AsB,BGL,BSA,BSB,BgL,BsA,BsB,CA2,CGL,CSA,CSB,CgL,CsA,CsB,DGL,DSA,DSB,DgL,DsA,DsB,EGL,ESA,ESB,EgL,EsA,EsB,FGL,FSA,FSB,FgL,FsA,FsB,GGL,GSA,GSB,GgL,GsA,GsB,HGL,HSA,HSB,HgL,HsA,HsB,IGL,ISA,ISB,IgL,IsA,IsB,JGL,JSA,JSB,JgL,JsA,JsB,KGL,KSA,KSB,KgL,KsA,KsB,MEX,NLN,OLS,OLT,OME,PEA,PEB,PGA,PGB,PKA,PKB,PLA,PLB,PMA,PMB,PNA,PNB,PTA,PTB,PeA,PeB,PgA,PgB,PkA,PkB,PlA,PlB,PmA,PmB,PnA,PnB,PtA,PtB,QBD,QBU,QCD,QCU,QEA,QEB,QGA,QGB,QJD,QJU,QKA,QKB,QLA,QLB,QMA,QMB,QNA,QNB,QPD,QPU,QTA,QTB,QVA,QVB,QWA,QWB,QYA,QYB,QbD,QbU,QcD,QcU,QeA,QeB,QgA,QgB,QjD,QjU,QkA,QkB,QlA,QlB,QmA,QmB,QnA,QnB,QpD,QpU,QtA,QtB,QvA,QvB,QwA,QwB,QyA,QyB,REA,REB,RGA,RGB,RKA,RKB,RLA,RLB,RMA,RMB,RNA,RNB,ROH,RTA,RTB,ReA,ReB,RgA,RgB,RkA,RkB,RlA,RlB,RmA,RmB,RnA,RnB,RtA,RtB,SEA,SEB,SGA,SGB,SKA,SKB,SLA,SLB,SMA,SMB,SNA,SNB,SO3,STA,STB,SeA,SeB,SgA,SgB,SkA,SkB,SlA,SlB,SmA,SmB,SnA,SnB,StA,StB,TAA,TAB,TBT,TDA,TDB,TEA,TEB,TFA,TFB,TGA,TGB,THA,THB,TKA,TKB,TLA,TLB,TMA,TMB,TNA,TNB,TOA,TOB,TQA,TQB,TRA,TRB,TTA,TTB,TUA,TUB,TXA,TXB,TZA,TZB,TaA,TaB,TdA,TdB,TeA,TeB,TfA,TfB,TgA,TgB,ThA,ThB,TkA,TkB,TlA,TlB,TmA,TmB,TnA,TnB,ToA,ToB,TqA,TqB,TrA,TrB,TtA,TtB,TuA,TuB,TxA,TxB,TzA,TzB,UBD,UBU,UCD,UCU,UEA,UEB,UGA,UGB,UJD,UJU,UKA,UKB,ULA,ULB,UMA,UMB,UNA,UNB,UPD,UPU,UTA,UTB,UVA,UVB,UWA,UWB,UYA,UYB,UbD,UbU,UcD,UcU,UeA,UeB,UgA,UgB,UjD,UjU,UkA,UkB,UlA,UlB,UmA,UmB,UnA,UnB,UpD,UpU,UtA,UtB,UvA,UvB,UwA,UwB,UyA,UyB,VBD,VBU,VCD,VCU,VEA,VEB,VGA,VGB,VJD,VJU,VKA,VKB,VLA,VLB,VMA,VMB,VNA,VNB,VPD,VPU,VTA,VTB,VVA,VVB,VWA,VWB,VYA,VYB,VbD,VbU,VcD,VcU,VeA,VeB,VgA,VgB,VjD,VjU,VkA,VkB,VlA,VlB,VmA,VmB,VnA,VnB,VpD,VpU,VtA,VtB,VvA,VvB,VwA,VwB,VyA,VyB,WAA,WAB,WBA,WBB,WBD,WBU,WCA,WCB,WCD,WCU,WDA,WDB,WEA,WEB,WFA,WFB,WGA,WGB,WHA,WHB,WJA,WJB,WJD,WJU,WKA,WKB,WLA,WLB,WMA,WMB,WNA,WNB,WOA,WOB,WPA,WPB,WPD,WPU,WQA,WQB,WRA,WRB,WTA,WTB,WUA,WUB,WVA,WVB,WWA,WWB,WXA,WXB,WYA,WYB,WZA,WZB,WaA,WaB,WbA,WbB,WbD,WbU,WcA,WcB,WcD,WcU,WdA,WdB,WdR,WeA,WeB,WfA,WfB,WgA,WgB,WhA,WhB,WjA,WjB,WjD,WjU,WkA,WkB,WlA,WlB,WmA,WmB,WnA,WnB,WoA,WoB,WpA,WpB,WpD,WpU,WqA,WqB,WrA,WrB,WtA,WtB,WuA,WuB,WvA,WvB,WwA,WwB,WxA,WxB,WyA,WyB,WzA,WzB,XEA,XEB,XGA,XGB,XKA,XKB,XLA,XLB,XMA,XMB,XNA,XNB,XTA,XTB,XeA,XeB,XgA,XgB,XkA,XkB,XlA,XlB,XmA,XmB,XnA,XnB,XtA,XtB,YAA,YAB,YAF,YDA,YDB,YEA,YEB,YFA,YFB,YGA,YGB,YGa,YHA,YHB,YKA,YKB,YLA,YLB,YMA,YMB,YNA,YNB,YOA,YOB,YQA,YQB,YRA,YRB,YTA,YTB,YTV,YTv,YUA,YUB,YXA,YXB,YZA,YZB,YaA,YaB,YdA,YdB,YeA,YeB,YfA,YfB,YgA,YgB,YhA,YhB,YkA,YkB,YlA,YlB,YmA,YmB,YnA,YnB,YoA,YoB,YqA,YqB,YrA,YrB,YtA,YtB,YtV,Ytv,YuA,YuB,YxA,YxB,YzA,YzB,ZAA,ZAB,ZAD,ZAU,ZDA,ZDB,ZDD,ZDU,ZEA,ZEB,ZFA,ZFB,ZGA,ZGB,ZHA,ZHB,ZKA,ZKB,ZLA,ZLB,ZMA,ZMB,ZNA,ZNB,ZOA,ZOB,ZOL,ZOL,ZQA,ZQB,ZRA,ZRB,ZRD,ZRU,ZTA,ZTB,ZUA,ZUB,ZXA,ZXB,ZXD,ZXU,ZZA,ZZB,ZaA,ZaB,ZaD,ZaU,ZdA,ZdB,ZdD,ZdU,ZeA,ZeB,ZfA,ZfB,ZgA,ZgB,ZhA,ZhB,ZkA,ZkB,ZlA,ZlB,ZmA,ZmB,ZnA,ZnB,ZoA,ZoB,ZqA,ZqB,ZrA,ZrB,ZrD,ZrU,ZtA,ZtB,ZuA,ZuB,ZxA,ZxB,ZxD,ZxU,ZzA,ZzB,HYP,NLN,OLP,OLS,OLT,CHYP,CNLN,COLP,COLS,COLT,NHYP,NNLN,NOLP,NOLS,NOLT', pt_topology)\n",
        "  first_atom1 = restraint_array1[0]\n",
        "  last_atom1 = restraint_array1[-1]\n",
        "  mask_protein_hbond = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "  mask_end = mask_protein_hbond\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  pt_topology = traj_load.top\n",
        "  restraint_array1 = pt.select_atoms('!(:WAT) & !(:Na+) & !(:Cl-) & !(:Mg+) & !(:K+)', pt_topology)\n",
        "  first_atom1 = restraint_array1[0]\n",
        "  last_atom1 = restraint_array1[-1]\n",
        "  mask_protein_glycan_hbond = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "  mask_end = mask_protein_glycan_hbond\n",
        "else:\n",
        "  pt_topology = traj_load.top\n",
        "  restraint_array1 = pt.select_atoms(':0AA,0AB,0AD,0AE,0AF,0AU,0BA,0BB,0BC,0BD,0BU,0CA,0CB,0CD,0CU,0DA,0DB,0DD,0DU,0EA,0EB,0FA,0FB,0GA,0GB,0GL,0HA,0HB,0JA,0JB,0JD,0JU,0KA,0KB,0LA,0LB,0MA,0MB,0NA,0NB,0OA,0OB,0PA,0PB,0PD,0PU,0QA,0QB,0RA,0RB,0RD,0RU,0SA,0SB,0TA,0TB,0TV,0Tv,0UA,0UB,0VA,0VB,0WA,0WB,0XA,0XB,0XD,0XU,0YA,0YB,0ZA,0ZB,0aA,0aB,0aD,0aU,0bA,0bB,0bC,0bD,0bU,0cA,0cB,0cD,0cU,0dA,0dB,0dD,0dR,0dU,0eA,0eB,0fA,0fB,0gA,0gB,0gL,0hA,0hB,0jA,0jB,0jD,0jU,0kA,0kB,0lA,0lB,0mA,0mB,0nA,0nB,0oA,0oB,0pA,0pB,0pD,0pU,0qA,0qB,0rA,0rB,0rD,0rU,0sA,0sB,0tA,0tB,0tV,0tv,0uA,0uB,0vA,0vB,0wA,0wB,0xA,0xB,0xD,0xU,0yA,0yB,0zA,0zB,1AA,1AB,1AD,1AU,1BA,1BB,1BD,1BU,1CA,1CB,1CD,1CU,1DA,1DB,1DD,1DU,1EA,1EB,1FA,1FB,1GA,1GB,1HA,1HB,1JA,1JB,1JD,1JU,1KA,1KB,1LA,1LB,1MA,1MB,1NA,1NB,1OA,1OB,1PA,1PB,1PD,1PU,1QA,1QB,1RA,1RB,1RD,1RU,1TA,1TB,1TV,1Tv,1UA,1UB,1VA,1VB,1WA,1WB,1XA,1XB,1XD,1XU,1YA,1YB,1ZA,1ZB,1aA,1aB,1aD,1aU,1bA,1bB,1bD,1bU,1cA,1cB,1cD,1cU,1dA,1dB,1dD,1dU,1eA,1eB,1fA,1fB,1gA,1gB,1hA,1hB,1jA,1jB,1jD,1jU,1kA,1kB,1lA,1lB,1mA,1mB,1nA,1nB,1oA,1oB,1pA,1pB,1pD,1pU,1qA,1qB,1rA,1rB,1rD,1rU,1tA,1tB,1tV,1tv,1uA,1uB,1vA,1vB,1wA,1wB,1xA,1xB,1xD,1xU,1yA,1yB,1zA,1zB,2AA,2AB,2AD,2AE,2AF,2AU,2BA,2BB,2BD,2BU,2CA,2CB,2CD,2CU,2DA,2DB,2DD,2DU,2EA,2EB,2FA,2FB,2GA,2GB,2HA,2HB,2JA,2JB,2JD,2JU,2KA,2KB,2LA,2LB,2MA,2MB,2NA,2NB,2OA,2OB,2PA,2PB,2PD,2PU,2QA,2QB,2RA,2RB,2RD,2RU,2TA,2TB,2TV,2Tv,2UA,2UB,2XA,2XB,2XD,2XU,2ZA,2ZB,2aA,2aB,2aD,2aU,2bA,2bB,2bD,2bU,2cA,2cB,2cD,2cU,2dA,2dB,2dD,2dU,2eA,2eB,2fA,2fB,2gA,2gB,2hA,2hB,2jA,2jB,2jD,2jU,2kA,2kB,2lA,2lB,2mA,2mB,2nA,2nB,2oA,2oB,2pA,2pB,2pD,2pU,2qA,2qB,2rA,2rB,2rD,2rU,2tA,2tB,2tV,2tv,2uA,2uB,2xA,2xB,2xD,2xU,2zA,2zB,3AA,3AB,3AD,3AU,3BA,3BB,3BC,3BD,3BU,3CA,3CB,3CD,3CU,3DA,3DB,3DD,3DU,3EA,3EB,3FA,3FB,3GA,3GB,3HA,3HB,3JA,3JB,3JD,3JU,3KA,3KB,3LA,3LB,3MA,3MB,3NA,3NB,3OA,3OB,3PA,3PB,3PD,3PU,3QA,3QB,3RA,3RB,3RD,3RU,3TA,3TB,3UA,3UB,3VA,3VB,3WA,3WB,3XA,3XB,3XD,3XU,3YA,3YB,3ZA,3ZB,3aA,3aB,3aD,3aU,3bA,3bB,3bC,3bD,3bU,3cA,3cB,3cD,3cU,3dA,3dB,3dD,3dR,3dU,3eA,3eB,3fA,3fB,3gA,3gB,3hA,3hB,3jA,3jB,3jD,3jU,3kA,3kB,3lA,3lB,3mA,3mB,3nA,3nB,3oA,3oB,3pA,3pB,3pD,3pU,3qA,3qB,3rA,3rB,3rD,3rU,3tA,3tB,3uA,3uB,3vA,3vB,3wA,3wB,3xA,3xB,3xD,3xU,3yA,3yB,3zA,3zB,4AA,4AB,4AE,4AF,4BA,4BB,4BD,4BU,4CA,4CB,4CD,4CU,4DA,4DB,4EA,4EB,4FA,4FB,4GA,4GB,4GL,4HA,4HB,4JA,4JB,4JD,4JU,4KA,4KB,4LA,4LB,4MA,4MB,4NA,4NB,4OA,4OB,4PA,4PB,4PD,4PU,4QA,4QB,4RA,4RB,4SA,4SB,4TA,4TB,4TV,4Tv,4UA,4UB,4VA,4VB,4WA,4WB,4XA,4XB,4YA,4YB,4ZA,4ZB,4aA,4aB,4bA,4bB,4bD,4bU,4cA,4cB,4cD,4cU,4dA,4dB,4dR,4eA,4eB,4fA,4fB,4gA,4gB,4gL,4hA,4hB,4jA,4jB,4jD,4jU,4kA,4kB,4lA,4lB,4mA,4mB,4nA,4nB,4oA,4oB,4pA,4pB,4pD,4pU,4qA,4qB,4rA,4rB,4sA,4sB,4tA,4tB,4tV,4tv,4uA,4uB,4vA,4vB,4wA,4wB,4xA,4xB,4yA,4yB,4zA,4zB,5AD,5AU,5BA,5BB,5CA,5CB,5DD,5DU,5JA,5JB,5PA,5PB,5RD,5RU,5XD,5XU,5aD,5aU,5bA,5bB,5cA,5cB,5dD,5dU,5jA,5jB,5pA,5pB,5rD,5rU,5xD,5xU,6BD,6BU,6CD,6CU,6EA,6EB,6GA,6GB,6JD,6JU,6KA,6KB,6LA,6LB,6MA,6MB,6NA,6NB,6PD,6PU,6TA,6TB,6VA,6VB,6WA,6WB,6YA,6YB,6bD,6bU,6cD,6cU,6eA,6eB,6gA,6gB,6jD,6jU,6kA,6kB,6lA,6lB,6mA,6mB,6nA,6nB,6pD,6pU,6tA,6tB,6vA,6vB,6wA,6wB,6yA,6yB,7GL,7SA,7SB,7gL,7sA,7sB,8GL,8SA,8SB,8gL,8sA,8sB,9GL,9SA,9SB,9gL,9sA,9sB,ACX,AGL,ASA,ASB,AgL,AsA,AsB,BGL,BSA,BSB,BgL,BsA,BsB,CA2,CGL,CSA,CSB,CgL,CsA,CsB,DGL,DSA,DSB,DgL,DsA,DsB,EGL,ESA,ESB,EgL,EsA,EsB,FGL,FSA,FSB,FgL,FsA,FsB,GGL,GSA,GSB,GgL,GsA,GsB,HGL,HSA,HSB,HgL,HsA,HsB,IGL,ISA,ISB,IgL,IsA,IsB,JGL,JSA,JSB,JgL,JsA,JsB,KGL,KSA,KSB,KgL,KsA,KsB,MEX,NLN,OLS,OLT,OME,PEA,PEB,PGA,PGB,PKA,PKB,PLA,PLB,PMA,PMB,PNA,PNB,PTA,PTB,PeA,PeB,PgA,PgB,PkA,PkB,PlA,PlB,PmA,PmB,PnA,PnB,PtA,PtB,QBD,QBU,QCD,QCU,QEA,QEB,QGA,QGB,QJD,QJU,QKA,QKB,QLA,QLB,QMA,QMB,QNA,QNB,QPD,QPU,QTA,QTB,QVA,QVB,QWA,QWB,QYA,QYB,QbD,QbU,QcD,QcU,QeA,QeB,QgA,QgB,QjD,QjU,QkA,QkB,QlA,QlB,QmA,QmB,QnA,QnB,QpD,QpU,QtA,QtB,QvA,QvB,QwA,QwB,QyA,QyB,REA,REB,RGA,RGB,RKA,RKB,RLA,RLB,RMA,RMB,RNA,RNB,ROH,RTA,RTB,ReA,ReB,RgA,RgB,RkA,RkB,RlA,RlB,RmA,RmB,RnA,RnB,RtA,RtB,SEA,SEB,SGA,SGB,SKA,SKB,SLA,SLB,SMA,SMB,SNA,SNB,SO3,STA,STB,SeA,SeB,SgA,SgB,SkA,SkB,SlA,SlB,SmA,SmB,SnA,SnB,StA,StB,TAA,TAB,TBT,TDA,TDB,TEA,TEB,TFA,TFB,TGA,TGB,THA,THB,TKA,TKB,TLA,TLB,TMA,TMB,TNA,TNB,TOA,TOB,TQA,TQB,TRA,TRB,TTA,TTB,TUA,TUB,TXA,TXB,TZA,TZB,TaA,TaB,TdA,TdB,TeA,TeB,TfA,TfB,TgA,TgB,ThA,ThB,TkA,TkB,TlA,TlB,TmA,TmB,TnA,TnB,ToA,ToB,TqA,TqB,TrA,TrB,TtA,TtB,TuA,TuB,TxA,TxB,TzA,TzB,UBD,UBU,UCD,UCU,UEA,UEB,UGA,UGB,UJD,UJU,UKA,UKB,ULA,ULB,UMA,UMB,UNA,UNB,UPD,UPU,UTA,UTB,UVA,UVB,UWA,UWB,UYA,UYB,UbD,UbU,UcD,UcU,UeA,UeB,UgA,UgB,UjD,UjU,UkA,UkB,UlA,UlB,UmA,UmB,UnA,UnB,UpD,UpU,UtA,UtB,UvA,UvB,UwA,UwB,UyA,UyB,VBD,VBU,VCD,VCU,VEA,VEB,VGA,VGB,VJD,VJU,VKA,VKB,VLA,VLB,VMA,VMB,VNA,VNB,VPD,VPU,VTA,VTB,VVA,VVB,VWA,VWB,VYA,VYB,VbD,VbU,VcD,VcU,VeA,VeB,VgA,VgB,VjD,VjU,VkA,VkB,VlA,VlB,VmA,VmB,VnA,VnB,VpD,VpU,VtA,VtB,VvA,VvB,VwA,VwB,VyA,VyB,WAA,WAB,WBA,WBB,WBD,WBU,WCA,WCB,WCD,WCU,WDA,WDB,WEA,WEB,WFA,WFB,WGA,WGB,WHA,WHB,WJA,WJB,WJD,WJU,WKA,WKB,WLA,WLB,WMA,WMB,WNA,WNB,WOA,WOB,WPA,WPB,WPD,WPU,WQA,WQB,WRA,WRB,WTA,WTB,WUA,WUB,WVA,WVB,WWA,WWB,WXA,WXB,WYA,WYB,WZA,WZB,WaA,WaB,WbA,WbB,WbD,WbU,WcA,WcB,WcD,WcU,WdA,WdB,WdR,WeA,WeB,WfA,WfB,WgA,WgB,WhA,WhB,WjA,WjB,WjD,WjU,WkA,WkB,WlA,WlB,WmA,WmB,WnA,WnB,WoA,WoB,WpA,WpB,WpD,WpU,WqA,WqB,WrA,WrB,WtA,WtB,WuA,WuB,WvA,WvB,WwA,WwB,WxA,WxB,WyA,WyB,WzA,WzB,XEA,XEB,XGA,XGB,XKA,XKB,XLA,XLB,XMA,XMB,XNA,XNB,XTA,XTB,XeA,XeB,XgA,XgB,XkA,XkB,XlA,XlB,XmA,XmB,XnA,XnB,XtA,XtB,YAA,YAB,YAF,YDA,YDB,YEA,YEB,YFA,YFB,YGA,YGB,YGa,YHA,YHB,YKA,YKB,YLA,YLB,YMA,YMB,YNA,YNB,YOA,YOB,YQA,YQB,YRA,YRB,YTA,YTB,YTV,YTv,YUA,YUB,YXA,YXB,YZA,YZB,YaA,YaB,YdA,YdB,YeA,YeB,YfA,YfB,YgA,YgB,YhA,YhB,YkA,YkB,YlA,YlB,YmA,YmB,YnA,YnB,YoA,YoB,YqA,YqB,YrA,YrB,YtA,YtB,YtV,Ytv,YuA,YuB,YxA,YxB,YzA,YzB,ZAA,ZAB,ZAD,ZAU,ZDA,ZDB,ZDD,ZDU,ZEA,ZEB,ZFA,ZFB,ZGA,ZGB,ZHA,ZHB,ZKA,ZKB,ZLA,ZLB,ZMA,ZMB,ZNA,ZNB,ZOA,ZOB,ZOL,ZOL,ZQA,ZQB,ZRA,ZRB,ZRD,ZRU,ZTA,ZTB,ZUA,ZUB,ZXA,ZXB,ZXD,ZXU,ZZA,ZZB,ZaA,ZaB,ZaD,ZaU,ZdA,ZdB,ZdD,ZdU,ZeA,ZeB,ZfA,ZfB,ZgA,ZgB,ZhA,ZhB,ZkA,ZkB,ZlA,ZlB,ZmA,ZmB,ZnA,ZnB,ZoA,ZoB,ZqA,ZqB,ZrA,ZrB,ZrD,ZrU,ZtA,ZtB,ZuA,ZuB,ZxA,ZxB,ZxD,ZxU,ZzA,ZzB,HYP,NLN,OLP,OLS,OLT,CHYP,CNLN,COLP,COLS,COLT,NHYP,NNLN,NOLP,NOLS,NOLT', pt_topology)\n",
        "  first_atom1 = restraint_array1[0]\n",
        "  last_atom1 = restraint_array1[-1]\n",
        "  mask_glycan_hbond = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "  mask_end = mask_glycan_hbond\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'hbond_solv' #@param {type:\"string\"}\n",
        "mask_hbond1 = 'acceptormask ' + mask_end + ' donormask :WAT'\n",
        "mask_hbond2 = 'acceptormask :WAT' + ' donormask ' + mask_end\n",
        "\n",
        "\n",
        "hb1 = pt.hbond(traj_load, options=mask_hbond1, dtype='dict', distance =3.5)\n",
        "hb2 = pt.hbond(traj_load, options=mask_hbond2, dtype='dict', distance =3.5)\n",
        "\n",
        "hb_end1 = hb1['total_solute_hbonds']\n",
        "hb_end2 = hb2['total_solute_hbonds']\n",
        "hb_total = hb_end1 + hb_end2\n",
        "time = len(hb_total)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "# hb.keys()\n",
        "# print(hb_total)\n",
        "# hb_end_mean = mean(hb_total)\n",
        "# hb_end_stdev = stdev(hb_total)\n",
        "# print(\"Hydrogen Bonds Average = \" + str(\"{:.2f}\".format(hb_end_mean)) + \" \\u00B1 \" + str(\"{:.2f}\".format(hb_end_stdev)))\n",
        "\n",
        "# Plotting:\n",
        "ax = plt.plot(time_array, hb_total, alpha=0.6, color = 'black', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"Hydrogen Bonds\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(hb_total)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "metadata": {
        "cellView": "form",
        "id": "W6StVVWtyqMO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Compute intramolecular hydrogen bonds**\n",
        "\n",
        "Selection_1 = \"Protein\" #@param [\"Protein\", \"Glycan\"]\n",
        "Selection_2 = \"Glycan\" #@param [\"Protein\", \"Glycan\"]\n",
        "\n",
        "pt_topology = traj_load.top\n",
        "restraint_array1 = pt.select_atoms('!(:WAT) & !(:Na+) & !(:Cl-) & !(:Mg+) & !(:K+) & !:0AA,0AB,0AD,0AE,0AF,0AU,0BA,0BB,0BC,0BD,0BU,0CA,0CB,0CD,0CU,0DA,0DB,0DD,0DU,0EA,0EB,0FA,0FB,0GA,0GB,0GL,0HA,0HB,0JA,0JB,0JD,0JU,0KA,0KB,0LA,0LB,0MA,0MB,0NA,0NB,0OA,0OB,0PA,0PB,0PD,0PU,0QA,0QB,0RA,0RB,0RD,0RU,0SA,0SB,0TA,0TB,0TV,0Tv,0UA,0UB,0VA,0VB,0WA,0WB,0XA,0XB,0XD,0XU,0YA,0YB,0ZA,0ZB,0aA,0aB,0aD,0aU,0bA,0bB,0bC,0bD,0bU,0cA,0cB,0cD,0cU,0dA,0dB,0dD,0dR,0dU,0eA,0eB,0fA,0fB,0gA,0gB,0gL,0hA,0hB,0jA,0jB,0jD,0jU,0kA,0kB,0lA,0lB,0mA,0mB,0nA,0nB,0oA,0oB,0pA,0pB,0pD,0pU,0qA,0qB,0rA,0rB,0rD,0rU,0sA,0sB,0tA,0tB,0tV,0tv,0uA,0uB,0vA,0vB,0wA,0wB,0xA,0xB,0xD,0xU,0yA,0yB,0zA,0zB,1AA,1AB,1AD,1AU,1BA,1BB,1BD,1BU,1CA,1CB,1CD,1CU,1DA,1DB,1DD,1DU,1EA,1EB,1FA,1FB,1GA,1GB,1HA,1HB,1JA,1JB,1JD,1JU,1KA,1KB,1LA,1LB,1MA,1MB,1NA,1NB,1OA,1OB,1PA,1PB,1PD,1PU,1QA,1QB,1RA,1RB,1RD,1RU,1TA,1TB,1TV,1Tv,1UA,1UB,1VA,1VB,1WA,1WB,1XA,1XB,1XD,1XU,1YA,1YB,1ZA,1ZB,1aA,1aB,1aD,1aU,1bA,1bB,1bD,1bU,1cA,1cB,1cD,1cU,1dA,1dB,1dD,1dU,1eA,1eB,1fA,1fB,1gA,1gB,1hA,1hB,1jA,1jB,1jD,1jU,1kA,1kB,1lA,1lB,1mA,1mB,1nA,1nB,1oA,1oB,1pA,1pB,1pD,1pU,1qA,1qB,1rA,1rB,1rD,1rU,1tA,1tB,1tV,1tv,1uA,1uB,1vA,1vB,1wA,1wB,1xA,1xB,1xD,1xU,1yA,1yB,1zA,1zB,2AA,2AB,2AD,2AE,2AF,2AU,2BA,2BB,2BD,2BU,2CA,2CB,2CD,2CU,2DA,2DB,2DD,2DU,2EA,2EB,2FA,2FB,2GA,2GB,2HA,2HB,2JA,2JB,2JD,2JU,2KA,2KB,2LA,2LB,2MA,2MB,2NA,2NB,2OA,2OB,2PA,2PB,2PD,2PU,2QA,2QB,2RA,2RB,2RD,2RU,2TA,2TB,2TV,2Tv,2UA,2UB,2XA,2XB,2XD,2XU,2ZA,2ZB,2aA,2aB,2aD,2aU,2bA,2bB,2bD,2bU,2cA,2cB,2cD,2cU,2dA,2dB,2dD,2dU,2eA,2eB,2fA,2fB,2gA,2gB,2hA,2hB,2jA,2jB,2jD,2jU,2kA,2kB,2lA,2lB,2mA,2mB,2nA,2nB,2oA,2oB,2pA,2pB,2pD,2pU,2qA,2qB,2rA,2rB,2rD,2rU,2tA,2tB,2tV,2tv,2uA,2uB,2xA,2xB,2xD,2xU,2zA,2zB,3AA,3AB,3AD,3AU,3BA,3BB,3BC,3BD,3BU,3CA,3CB,3CD,3CU,3DA,3DB,3DD,3DU,3EA,3EB,3FA,3FB,3GA,3GB,3HA,3HB,3JA,3JB,3JD,3JU,3KA,3KB,3LA,3LB,3MA,3MB,3NA,3NB,3OA,3OB,3PA,3PB,3PD,3PU,3QA,3QB,3RA,3RB,3RD,3RU,3TA,3TB,3UA,3UB,3VA,3VB,3WA,3WB,3XA,3XB,3XD,3XU,3YA,3YB,3ZA,3ZB,3aA,3aB,3aD,3aU,3bA,3bB,3bC,3bD,3bU,3cA,3cB,3cD,3cU,3dA,3dB,3dD,3dR,3dU,3eA,3eB,3fA,3fB,3gA,3gB,3hA,3hB,3jA,3jB,3jD,3jU,3kA,3kB,3lA,3lB,3mA,3mB,3nA,3nB,3oA,3oB,3pA,3pB,3pD,3pU,3qA,3qB,3rA,3rB,3rD,3rU,3tA,3tB,3uA,3uB,3vA,3vB,3wA,3wB,3xA,3xB,3xD,3xU,3yA,3yB,3zA,3zB,4AA,4AB,4AE,4AF,4BA,4BB,4BD,4BU,4CA,4CB,4CD,4CU,4DA,4DB,4EA,4EB,4FA,4FB,4GA,4GB,4GL,4HA,4HB,4JA,4JB,4JD,4JU,4KA,4KB,4LA,4LB,4MA,4MB,4NA,4NB,4OA,4OB,4PA,4PB,4PD,4PU,4QA,4QB,4RA,4RB,4SA,4SB,4TA,4TB,4TV,4Tv,4UA,4UB,4VA,4VB,4WA,4WB,4XA,4XB,4YA,4YB,4ZA,4ZB,4aA,4aB,4bA,4bB,4bD,4bU,4cA,4cB,4cD,4cU,4dA,4dB,4dR,4eA,4eB,4fA,4fB,4gA,4gB,4gL,4hA,4hB,4jA,4jB,4jD,4jU,4kA,4kB,4lA,4lB,4mA,4mB,4nA,4nB,4oA,4oB,4pA,4pB,4pD,4pU,4qA,4qB,4rA,4rB,4sA,4sB,4tA,4tB,4tV,4tv,4uA,4uB,4vA,4vB,4wA,4wB,4xA,4xB,4yA,4yB,4zA,4zB,5AD,5AU,5BA,5BB,5CA,5CB,5DD,5DU,5JA,5JB,5PA,5PB,5RD,5RU,5XD,5XU,5aD,5aU,5bA,5bB,5cA,5cB,5dD,5dU,5jA,5jB,5pA,5pB,5rD,5rU,5xD,5xU,6BD,6BU,6CD,6CU,6EA,6EB,6GA,6GB,6JD,6JU,6KA,6KB,6LA,6LB,6MA,6MB,6NA,6NB,6PD,6PU,6TA,6TB,6VA,6VB,6WA,6WB,6YA,6YB,6bD,6bU,6cD,6cU,6eA,6eB,6gA,6gB,6jD,6jU,6kA,6kB,6lA,6lB,6mA,6mB,6nA,6nB,6pD,6pU,6tA,6tB,6vA,6vB,6wA,6wB,6yA,6yB,7GL,7SA,7SB,7gL,7sA,7sB,8GL,8SA,8SB,8gL,8sA,8sB,9GL,9SA,9SB,9gL,9sA,9sB,ACX,AGL,ASA,ASB,AgL,AsA,AsB,BGL,BSA,BSB,BgL,BsA,BsB,CA2,CGL,CSA,CSB,CgL,CsA,CsB,DGL,DSA,DSB,DgL,DsA,DsB,EGL,ESA,ESB,EgL,EsA,EsB,FGL,FSA,FSB,FgL,FsA,FsB,GGL,GSA,GSB,GgL,GsA,GsB,HGL,HSA,HSB,HgL,HsA,HsB,IGL,ISA,ISB,IgL,IsA,IsB,JGL,JSA,JSB,JgL,JsA,JsB,KGL,KSA,KSB,KgL,KsA,KsB,MEX,NLN,OLS,OLT,OME,PEA,PEB,PGA,PGB,PKA,PKB,PLA,PLB,PMA,PMB,PNA,PNB,PTA,PTB,PeA,PeB,PgA,PgB,PkA,PkB,PlA,PlB,PmA,PmB,PnA,PnB,PtA,PtB,QBD,QBU,QCD,QCU,QEA,QEB,QGA,QGB,QJD,QJU,QKA,QKB,QLA,QLB,QMA,QMB,QNA,QNB,QPD,QPU,QTA,QTB,QVA,QVB,QWA,QWB,QYA,QYB,QbD,QbU,QcD,QcU,QeA,QeB,QgA,QgB,QjD,QjU,QkA,QkB,QlA,QlB,QmA,QmB,QnA,QnB,QpD,QpU,QtA,QtB,QvA,QvB,QwA,QwB,QyA,QyB,REA,REB,RGA,RGB,RKA,RKB,RLA,RLB,RMA,RMB,RNA,RNB,ROH,RTA,RTB,ReA,ReB,RgA,RgB,RkA,RkB,RlA,RlB,RmA,RmB,RnA,RnB,RtA,RtB,SEA,SEB,SGA,SGB,SKA,SKB,SLA,SLB,SMA,SMB,SNA,SNB,SO3,STA,STB,SeA,SeB,SgA,SgB,SkA,SkB,SlA,SlB,SmA,SmB,SnA,SnB,StA,StB,TAA,TAB,TBT,TDA,TDB,TEA,TEB,TFA,TFB,TGA,TGB,THA,THB,TKA,TKB,TLA,TLB,TMA,TMB,TNA,TNB,TOA,TOB,TQA,TQB,TRA,TRB,TTA,TTB,TUA,TUB,TXA,TXB,TZA,TZB,TaA,TaB,TdA,TdB,TeA,TeB,TfA,TfB,TgA,TgB,ThA,ThB,TkA,TkB,TlA,TlB,TmA,TmB,TnA,TnB,ToA,ToB,TqA,TqB,TrA,TrB,TtA,TtB,TuA,TuB,TxA,TxB,TzA,TzB,UBD,UBU,UCD,UCU,UEA,UEB,UGA,UGB,UJD,UJU,UKA,UKB,ULA,ULB,UMA,UMB,UNA,UNB,UPD,UPU,UTA,UTB,UVA,UVB,UWA,UWB,UYA,UYB,UbD,UbU,UcD,UcU,UeA,UeB,UgA,UgB,UjD,UjU,UkA,UkB,UlA,UlB,UmA,UmB,UnA,UnB,UpD,UpU,UtA,UtB,UvA,UvB,UwA,UwB,UyA,UyB,VBD,VBU,VCD,VCU,VEA,VEB,VGA,VGB,VJD,VJU,VKA,VKB,VLA,VLB,VMA,VMB,VNA,VNB,VPD,VPU,VTA,VTB,VVA,VVB,VWA,VWB,VYA,VYB,VbD,VbU,VcD,VcU,VeA,VeB,VgA,VgB,VjD,VjU,VkA,VkB,VlA,VlB,VmA,VmB,VnA,VnB,VpD,VpU,VtA,VtB,VvA,VvB,VwA,VwB,VyA,VyB,WAA,WAB,WBA,WBB,WBD,WBU,WCA,WCB,WCD,WCU,WDA,WDB,WEA,WEB,WFA,WFB,WGA,WGB,WHA,WHB,WJA,WJB,WJD,WJU,WKA,WKB,WLA,WLB,WMA,WMB,WNA,WNB,WOA,WOB,WPA,WPB,WPD,WPU,WQA,WQB,WRA,WRB,WTA,WTB,WUA,WUB,WVA,WVB,WWA,WWB,WXA,WXB,WYA,WYB,WZA,WZB,WaA,WaB,WbA,WbB,WbD,WbU,WcA,WcB,WcD,WcU,WdA,WdB,WdR,WeA,WeB,WfA,WfB,WgA,WgB,WhA,WhB,WjA,WjB,WjD,WjU,WkA,WkB,WlA,WlB,WmA,WmB,WnA,WnB,WoA,WoB,WpA,WpB,WpD,WpU,WqA,WqB,WrA,WrB,WtA,WtB,WuA,WuB,WvA,WvB,WwA,WwB,WxA,WxB,WyA,WyB,WzA,WzB,XEA,XEB,XGA,XGB,XKA,XKB,XLA,XLB,XMA,XMB,XNA,XNB,XTA,XTB,XeA,XeB,XgA,XgB,XkA,XkB,XlA,XlB,XmA,XmB,XnA,XnB,XtA,XtB,YAA,YAB,YAF,YDA,YDB,YEA,YEB,YFA,YFB,YGA,YGB,YGa,YHA,YHB,YKA,YKB,YLA,YLB,YMA,YMB,YNA,YNB,YOA,YOB,YQA,YQB,YRA,YRB,YTA,YTB,YTV,YTv,YUA,YUB,YXA,YXB,YZA,YZB,YaA,YaB,YdA,YdB,YeA,YeB,YfA,YfB,YgA,YgB,YhA,YhB,YkA,YkB,YlA,YlB,YmA,YmB,YnA,YnB,YoA,YoB,YqA,YqB,YrA,YrB,YtA,YtB,YtV,Ytv,YuA,YuB,YxA,YxB,YzA,YzB,ZAA,ZAB,ZAD,ZAU,ZDA,ZDB,ZDD,ZDU,ZEA,ZEB,ZFA,ZFB,ZGA,ZGB,ZHA,ZHB,ZKA,ZKB,ZLA,ZLB,ZMA,ZMB,ZNA,ZNB,ZOA,ZOB,ZOL,ZOL,ZQA,ZQB,ZRA,ZRB,ZRD,ZRU,ZTA,ZTB,ZUA,ZUB,ZXA,ZXB,ZXD,ZXU,ZZA,ZZB,ZaA,ZaB,ZaD,ZaU,ZdA,ZdB,ZdD,ZdU,ZeA,ZeB,ZfA,ZfB,ZgA,ZgB,ZhA,ZhB,ZkA,ZkB,ZlA,ZlB,ZmA,ZmB,ZnA,ZnB,ZoA,ZoB,ZqA,ZqB,ZrA,ZrB,ZrD,ZrU,ZtA,ZtB,ZuA,ZuB,ZxA,ZxB,ZxD,ZxU,ZzA,ZzB,HYP,NLN,OLP,OLS,OLT,CHYP,CNLN,COLP,COLS,COLT,NHYP,NNLN,NOLP,NOLS,NOLT', pt_topology)\n",
        "first_atom1 = restraint_array1[0]\n",
        "last_atom1 = restraint_array1[-1]\n",
        "mask_protein_hbond = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "\n",
        "pt_topology = traj_load.top\n",
        "restraint_array1 = pt.select_atoms(':0AA,0AB,0AD,0AE,0AF,0AU,0BA,0BB,0BC,0BD,0BU,0CA,0CB,0CD,0CU,0DA,0DB,0DD,0DU,0EA,0EB,0FA,0FB,0GA,0GB,0GL,0HA,0HB,0JA,0JB,0JD,0JU,0KA,0KB,0LA,0LB,0MA,0MB,0NA,0NB,0OA,0OB,0PA,0PB,0PD,0PU,0QA,0QB,0RA,0RB,0RD,0RU,0SA,0SB,0TA,0TB,0TV,0Tv,0UA,0UB,0VA,0VB,0WA,0WB,0XA,0XB,0XD,0XU,0YA,0YB,0ZA,0ZB,0aA,0aB,0aD,0aU,0bA,0bB,0bC,0bD,0bU,0cA,0cB,0cD,0cU,0dA,0dB,0dD,0dR,0dU,0eA,0eB,0fA,0fB,0gA,0gB,0gL,0hA,0hB,0jA,0jB,0jD,0jU,0kA,0kB,0lA,0lB,0mA,0mB,0nA,0nB,0oA,0oB,0pA,0pB,0pD,0pU,0qA,0qB,0rA,0rB,0rD,0rU,0sA,0sB,0tA,0tB,0tV,0tv,0uA,0uB,0vA,0vB,0wA,0wB,0xA,0xB,0xD,0xU,0yA,0yB,0zA,0zB,1AA,1AB,1AD,1AU,1BA,1BB,1BD,1BU,1CA,1CB,1CD,1CU,1DA,1DB,1DD,1DU,1EA,1EB,1FA,1FB,1GA,1GB,1HA,1HB,1JA,1JB,1JD,1JU,1KA,1KB,1LA,1LB,1MA,1MB,1NA,1NB,1OA,1OB,1PA,1PB,1PD,1PU,1QA,1QB,1RA,1RB,1RD,1RU,1TA,1TB,1TV,1Tv,1UA,1UB,1VA,1VB,1WA,1WB,1XA,1XB,1XD,1XU,1YA,1YB,1ZA,1ZB,1aA,1aB,1aD,1aU,1bA,1bB,1bD,1bU,1cA,1cB,1cD,1cU,1dA,1dB,1dD,1dU,1eA,1eB,1fA,1fB,1gA,1gB,1hA,1hB,1jA,1jB,1jD,1jU,1kA,1kB,1lA,1lB,1mA,1mB,1nA,1nB,1oA,1oB,1pA,1pB,1pD,1pU,1qA,1qB,1rA,1rB,1rD,1rU,1tA,1tB,1tV,1tv,1uA,1uB,1vA,1vB,1wA,1wB,1xA,1xB,1xD,1xU,1yA,1yB,1zA,1zB,2AA,2AB,2AD,2AE,2AF,2AU,2BA,2BB,2BD,2BU,2CA,2CB,2CD,2CU,2DA,2DB,2DD,2DU,2EA,2EB,2FA,2FB,2GA,2GB,2HA,2HB,2JA,2JB,2JD,2JU,2KA,2KB,2LA,2LB,2MA,2MB,2NA,2NB,2OA,2OB,2PA,2PB,2PD,2PU,2QA,2QB,2RA,2RB,2RD,2RU,2TA,2TB,2TV,2Tv,2UA,2UB,2XA,2XB,2XD,2XU,2ZA,2ZB,2aA,2aB,2aD,2aU,2bA,2bB,2bD,2bU,2cA,2cB,2cD,2cU,2dA,2dB,2dD,2dU,2eA,2eB,2fA,2fB,2gA,2gB,2hA,2hB,2jA,2jB,2jD,2jU,2kA,2kB,2lA,2lB,2mA,2mB,2nA,2nB,2oA,2oB,2pA,2pB,2pD,2pU,2qA,2qB,2rA,2rB,2rD,2rU,2tA,2tB,2tV,2tv,2uA,2uB,2xA,2xB,2xD,2xU,2zA,2zB,3AA,3AB,3AD,3AU,3BA,3BB,3BC,3BD,3BU,3CA,3CB,3CD,3CU,3DA,3DB,3DD,3DU,3EA,3EB,3FA,3FB,3GA,3GB,3HA,3HB,3JA,3JB,3JD,3JU,3KA,3KB,3LA,3LB,3MA,3MB,3NA,3NB,3OA,3OB,3PA,3PB,3PD,3PU,3QA,3QB,3RA,3RB,3RD,3RU,3TA,3TB,3UA,3UB,3VA,3VB,3WA,3WB,3XA,3XB,3XD,3XU,3YA,3YB,3ZA,3ZB,3aA,3aB,3aD,3aU,3bA,3bB,3bC,3bD,3bU,3cA,3cB,3cD,3cU,3dA,3dB,3dD,3dR,3dU,3eA,3eB,3fA,3fB,3gA,3gB,3hA,3hB,3jA,3jB,3jD,3jU,3kA,3kB,3lA,3lB,3mA,3mB,3nA,3nB,3oA,3oB,3pA,3pB,3pD,3pU,3qA,3qB,3rA,3rB,3rD,3rU,3tA,3tB,3uA,3uB,3vA,3vB,3wA,3wB,3xA,3xB,3xD,3xU,3yA,3yB,3zA,3zB,4AA,4AB,4AE,4AF,4BA,4BB,4BD,4BU,4CA,4CB,4CD,4CU,4DA,4DB,4EA,4EB,4FA,4FB,4GA,4GB,4GL,4HA,4HB,4JA,4JB,4JD,4JU,4KA,4KB,4LA,4LB,4MA,4MB,4NA,4NB,4OA,4OB,4PA,4PB,4PD,4PU,4QA,4QB,4RA,4RB,4SA,4SB,4TA,4TB,4TV,4Tv,4UA,4UB,4VA,4VB,4WA,4WB,4XA,4XB,4YA,4YB,4ZA,4ZB,4aA,4aB,4bA,4bB,4bD,4bU,4cA,4cB,4cD,4cU,4dA,4dB,4dR,4eA,4eB,4fA,4fB,4gA,4gB,4gL,4hA,4hB,4jA,4jB,4jD,4jU,4kA,4kB,4lA,4lB,4mA,4mB,4nA,4nB,4oA,4oB,4pA,4pB,4pD,4pU,4qA,4qB,4rA,4rB,4sA,4sB,4tA,4tB,4tV,4tv,4uA,4uB,4vA,4vB,4wA,4wB,4xA,4xB,4yA,4yB,4zA,4zB,5AD,5AU,5BA,5BB,5CA,5CB,5DD,5DU,5JA,5JB,5PA,5PB,5RD,5RU,5XD,5XU,5aD,5aU,5bA,5bB,5cA,5cB,5dD,5dU,5jA,5jB,5pA,5pB,5rD,5rU,5xD,5xU,6BD,6BU,6CD,6CU,6EA,6EB,6GA,6GB,6JD,6JU,6KA,6KB,6LA,6LB,6MA,6MB,6NA,6NB,6PD,6PU,6TA,6TB,6VA,6VB,6WA,6WB,6YA,6YB,6bD,6bU,6cD,6cU,6eA,6eB,6gA,6gB,6jD,6jU,6kA,6kB,6lA,6lB,6mA,6mB,6nA,6nB,6pD,6pU,6tA,6tB,6vA,6vB,6wA,6wB,6yA,6yB,7GL,7SA,7SB,7gL,7sA,7sB,8GL,8SA,8SB,8gL,8sA,8sB,9GL,9SA,9SB,9gL,9sA,9sB,ACX,AGL,ASA,ASB,AgL,AsA,AsB,BGL,BSA,BSB,BgL,BsA,BsB,CA2,CGL,CSA,CSB,CgL,CsA,CsB,DGL,DSA,DSB,DgL,DsA,DsB,EGL,ESA,ESB,EgL,EsA,EsB,FGL,FSA,FSB,FgL,FsA,FsB,GGL,GSA,GSB,GgL,GsA,GsB,HGL,HSA,HSB,HgL,HsA,HsB,IGL,ISA,ISB,IgL,IsA,IsB,JGL,JSA,JSB,JgL,JsA,JsB,KGL,KSA,KSB,KgL,KsA,KsB,MEX,NLN,OLS,OLT,OME,PEA,PEB,PGA,PGB,PKA,PKB,PLA,PLB,PMA,PMB,PNA,PNB,PTA,PTB,PeA,PeB,PgA,PgB,PkA,PkB,PlA,PlB,PmA,PmB,PnA,PnB,PtA,PtB,QBD,QBU,QCD,QCU,QEA,QEB,QGA,QGB,QJD,QJU,QKA,QKB,QLA,QLB,QMA,QMB,QNA,QNB,QPD,QPU,QTA,QTB,QVA,QVB,QWA,QWB,QYA,QYB,QbD,QbU,QcD,QcU,QeA,QeB,QgA,QgB,QjD,QjU,QkA,QkB,QlA,QlB,QmA,QmB,QnA,QnB,QpD,QpU,QtA,QtB,QvA,QvB,QwA,QwB,QyA,QyB,REA,REB,RGA,RGB,RKA,RKB,RLA,RLB,RMA,RMB,RNA,RNB,ROH,RTA,RTB,ReA,ReB,RgA,RgB,RkA,RkB,RlA,RlB,RmA,RmB,RnA,RnB,RtA,RtB,SEA,SEB,SGA,SGB,SKA,SKB,SLA,SLB,SMA,SMB,SNA,SNB,SO3,STA,STB,SeA,SeB,SgA,SgB,SkA,SkB,SlA,SlB,SmA,SmB,SnA,SnB,StA,StB,TAA,TAB,TBT,TDA,TDB,TEA,TEB,TFA,TFB,TGA,TGB,THA,THB,TKA,TKB,TLA,TLB,TMA,TMB,TNA,TNB,TOA,TOB,TQA,TQB,TRA,TRB,TTA,TTB,TUA,TUB,TXA,TXB,TZA,TZB,TaA,TaB,TdA,TdB,TeA,TeB,TfA,TfB,TgA,TgB,ThA,ThB,TkA,TkB,TlA,TlB,TmA,TmB,TnA,TnB,ToA,ToB,TqA,TqB,TrA,TrB,TtA,TtB,TuA,TuB,TxA,TxB,TzA,TzB,UBD,UBU,UCD,UCU,UEA,UEB,UGA,UGB,UJD,UJU,UKA,UKB,ULA,ULB,UMA,UMB,UNA,UNB,UPD,UPU,UTA,UTB,UVA,UVB,UWA,UWB,UYA,UYB,UbD,UbU,UcD,UcU,UeA,UeB,UgA,UgB,UjD,UjU,UkA,UkB,UlA,UlB,UmA,UmB,UnA,UnB,UpD,UpU,UtA,UtB,UvA,UvB,UwA,UwB,UyA,UyB,VBD,VBU,VCD,VCU,VEA,VEB,VGA,VGB,VJD,VJU,VKA,VKB,VLA,VLB,VMA,VMB,VNA,VNB,VPD,VPU,VTA,VTB,VVA,VVB,VWA,VWB,VYA,VYB,VbD,VbU,VcD,VcU,VeA,VeB,VgA,VgB,VjD,VjU,VkA,VkB,VlA,VlB,VmA,VmB,VnA,VnB,VpD,VpU,VtA,VtB,VvA,VvB,VwA,VwB,VyA,VyB,WAA,WAB,WBA,WBB,WBD,WBU,WCA,WCB,WCD,WCU,WDA,WDB,WEA,WEB,WFA,WFB,WGA,WGB,WHA,WHB,WJA,WJB,WJD,WJU,WKA,WKB,WLA,WLB,WMA,WMB,WNA,WNB,WOA,WOB,WPA,WPB,WPD,WPU,WQA,WQB,WRA,WRB,WTA,WTB,WUA,WUB,WVA,WVB,WWA,WWB,WXA,WXB,WYA,WYB,WZA,WZB,WaA,WaB,WbA,WbB,WbD,WbU,WcA,WcB,WcD,WcU,WdA,WdB,WdR,WeA,WeB,WfA,WfB,WgA,WgB,WhA,WhB,WjA,WjB,WjD,WjU,WkA,WkB,WlA,WlB,WmA,WmB,WnA,WnB,WoA,WoB,WpA,WpB,WpD,WpU,WqA,WqB,WrA,WrB,WtA,WtB,WuA,WuB,WvA,WvB,WwA,WwB,WxA,WxB,WyA,WyB,WzA,WzB,XEA,XEB,XGA,XGB,XKA,XKB,XLA,XLB,XMA,XMB,XNA,XNB,XTA,XTB,XeA,XeB,XgA,XgB,XkA,XkB,XlA,XlB,XmA,XmB,XnA,XnB,XtA,XtB,YAA,YAB,YAF,YDA,YDB,YEA,YEB,YFA,YFB,YGA,YGB,YGa,YHA,YHB,YKA,YKB,YLA,YLB,YMA,YMB,YNA,YNB,YOA,YOB,YQA,YQB,YRA,YRB,YTA,YTB,YTV,YTv,YUA,YUB,YXA,YXB,YZA,YZB,YaA,YaB,YdA,YdB,YeA,YeB,YfA,YfB,YgA,YgB,YhA,YhB,YkA,YkB,YlA,YlB,YmA,YmB,YnA,YnB,YoA,YoB,YqA,YqB,YrA,YrB,YtA,YtB,YtV,Ytv,YuA,YuB,YxA,YxB,YzA,YzB,ZAA,ZAB,ZAD,ZAU,ZDA,ZDB,ZDD,ZDU,ZEA,ZEB,ZFA,ZFB,ZGA,ZGB,ZHA,ZHB,ZKA,ZKB,ZLA,ZLB,ZMA,ZMB,ZNA,ZNB,ZOA,ZOB,ZOL,ZOL,ZQA,ZQB,ZRA,ZRB,ZRD,ZRU,ZTA,ZTB,ZUA,ZUB,ZXA,ZXB,ZXD,ZXU,ZZA,ZZB,ZaA,ZaB,ZaD,ZaU,ZdA,ZdB,ZdD,ZdU,ZeA,ZeB,ZfA,ZfB,ZgA,ZgB,ZhA,ZhB,ZkA,ZkB,ZlA,ZlB,ZmA,ZmB,ZnA,ZnB,ZoA,ZoB,ZqA,ZqB,ZrA,ZrB,ZrD,ZrU,ZtA,ZtB,ZuA,ZuB,ZxA,ZxB,ZxD,ZxU,ZzA,ZzB,HYP,NLN,OLP,OLS,OLT,CHYP,CNLN,COLP,COLS,COLT,NHYP,NNLN,NOLP,NOLS,NOLT', pt_topology)\n",
        "first_atom1 = restraint_array1[0]\n",
        "last_atom1 = restraint_array1[-1]\n",
        "mask_glycan_hbond = \"@\" + str(first_atom1+1) + \"-\" + str(last_atom1+1)\n",
        "\n",
        "if Selection_1 == \"Protein\"  and Selection_2 == \"Protein\":\n",
        "  mask_end = mask_protein_hbond\n",
        "  mask_hbond1 = 'acceptormask ' + mask_end + ' donormask ' + mask_end\n",
        "  hb1 = pt.hbond(traj_load, options=mask_hbond1, dtype='dict', distance =3.5)\n",
        "  hb = hb1['total_solute_hbonds']\n",
        "  # print(mask_end)  \n",
        "elif Selection_1 == \"Glycan\"  and Selection_2 == \"Glycan\":\n",
        "  mask_end = mask_glycan_hbond\n",
        "  mask_hbond1 = 'acceptormask ' + mask_end + ' donormask ' + mask_end\n",
        "  hb1 = pt.hbond(traj_load, options=mask_hbond1, dtype='dict', distance =3.5)\n",
        "  hb = hb1['total_solute_hbonds']\n",
        "  # print(mask_end)  \n",
        "else:\n",
        "  mask_end1 = mask_protein_hbond\n",
        "  mask_end2 = mask_glycan_hbond\n",
        "  mask_hbond1 = 'acceptormask ' + mask_end1 + ' donormask ' + mask_end2\n",
        "  mask_hbond2 = 'acceptormask ' + mask_end2 + ' donormask ' + mask_end1\n",
        "\n",
        "  hb1 = pt.hbond(traj_load, options=mask_hbond1, dtype='dict', distance =3.5)\n",
        "  hb2 = pt.hbond(traj_load, options=mask_hbond2, dtype='dict', distance =3.5)\n",
        "\n",
        "  hb_end1 = hb1['total_solute_hbonds']\n",
        "  hb_end2 = hb2['total_solute_hbonds']\n",
        "  hb = hb_end1 + hb_end2\n",
        "  # print(mask_end1,mask_end2)\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'hbond_intra' #@param {type:\"string\"}\n",
        "\n",
        "time = len(hb)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "# hb.keys()\n",
        "# print(hb_total)\n",
        "# hb_end_mean = mean(hb_total)\n",
        "# hb_end_stdev = stdev(hb_total)\n",
        "# print(\"Hydrogen Bonds Average = \" + str(\"{:.2f}\".format(hb_end_mean)) + \" \\u00B1 \" + str(\"{:.2f}\".format(hb_end_stdev)))\n",
        "\n",
        "# Plotting:\n",
        "ax = plt.plot(time_array, hb, alpha=0.6, color = 'black', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"Hydrogen Bonds\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(hb_total)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "metadata": {
        "cellView": "form",
        "id": "vpmKEYWi0TSs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Compute distance between the atoms**\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'distance_atoms' #@param {type:\"string\"}\n",
        "#@markdown **Type the number of atoms separated by commas and without spaces (1,2,3...):**\n",
        "Selection_1 = '3' #@param {type:\"string\"}\n",
        "Selection_2 = '12' #@param {type:\"string\"}\n",
        "\n",
        "mask = \"@\" + str(Selection_1) + \" @\" + str(Selection_2)\n",
        "dist = pt.distance(traj_load, mask)\n",
        "print(\"Selected atoms = \" + Selection_1 + \", \"  + Selection_2 + \"\\n\")\n",
        "dist_mean = mean(dist)\n",
        "dist_stdev = stdev(dist)\n",
        "print(\"Distance Average = \" + str(\"{:.2f}\".format(dist_mean)) + \" \\u00B1 \" + str(\"{:.2f}\".format(dist_stdev)) + \" Å\")\n",
        "\n",
        "time = len(dist)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "\n",
        "# Plotting:\n",
        "ax = plt.plot(time_array, dist, alpha=1, color = 'magenta', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "#plt.ylim(2, 6)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"Distance [$\\AA$]\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(dist)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "metadata": {
        "cellView": "form",
        "id": "wFaxFD_I5EXF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Plot distance as a distribution**\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'distance_dist' #@param {type:\"string\"}\n",
        "\n",
        "ax = sb.kdeplot(dist, color=\"magenta\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('Distance [$\\AA$]', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax.spines['top'].set_visible(False)\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['bottom'].set_visible(True)\n",
        "ax.spines['left'].set_visible(False)\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')"
      ],
      "metadata": {
        "cellView": "form",
        "id": "lflH87Ee5L8p"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Compute specific dihedral angles**\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'dihedral_angle' #@param {type:\"string\"}\n",
        "#@markdown **Type the number of atoms:**\n",
        "Atom_1 = '1' #@param {type:\"string\"}\n",
        "Atom_2 = '2' #@param {type:\"string\"}\n",
        "Atom_3 = '3' #@param {type:\"string\"}\n",
        "Atom_4 = '4' #@param {type:\"string\"}\n",
        "\n",
        "mask = \"@\" + str(Atom_1) + \" @\" + str(Atom_2) + \" @\" + str(Atom_3) + \" @\" + str(Atom_4)\n",
        "dih = pt.dihedral(traj_load, mask)\n",
        "print(\"Selected atoms = \" + Atom_1 + \", \" + Atom_2 + \", \" + Atom_3 + \", \" + Atom_4 + \"\\n\")\n",
        "dih_mean = mean(dih)\n",
        "dih_stdev = stdev(dih)\n",
        "print(\"Dihedral Angle Average = \" + str(\"{:.2f}\".format(dih_mean)) + \" \\u00B1 \" + str(\"{:.2f}\".format(dih_stdev)) + \"°\")\n",
        "\n",
        "time = len(dih)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "\n",
        "# Plotting:\n",
        "ax = plt.plot(time_array, dih, alpha=1, color = 'orange', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "#plt.ylim(2, 6)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"Dihedral Angle ($^\\circ$)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(dih)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "metadata": {
        "cellView": "form",
        "id": "qwOLZ8Kx5Pcc"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title **Plot dihedral angle as a distribution**\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'dihedral_dist' #@param {type:\"string\"}\n",
        "\n",
        "ax = sb.kdeplot(dih, color=\"orange\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('Dihedral Angle ($^\\circ$)', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax.spines['top'].set_visible(False)\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['bottom'].set_visible(True)\n",
        "ax.spines['left'].set_visible(False)\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')"
      ],
      "metadata": {
        "cellView": "form",
        "id": "P9K5Ai7K5Ts5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wBrBMF4Puyv6",
        "cellView": "form"
      },
      "source": [
        "#@title **Compute RMSD**\n",
        "\n",
        "Selection = \"Protein+Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "\n",
        "Output_name = 'rmsd_protein' #@param {type:\"string\"}\n",
        "\n",
        "rmsd = pt.rmsd(traj_load, ref = 0, mask = mask_end)\n",
        "\n",
        "time = len(rmsd)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "\n",
        "# Plotting:\n",
        "ax = plt.plot(time_array, rmsd, alpha=0.6, color = 'blue', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "#plt.ylim(2, 6)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"RMSD [$\\AA$]\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(rmsd)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZHyMpikjuaLT",
        "cellView": "form"
      },
      "source": [
        "#@title **Plot RMSD as a distribution**\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'rmsd_dist' #@param {type:\"string\"}\n",
        "\n",
        "ax = sb.kdeplot(rmsd, color=\"blue\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('RMSD [$\\AA$]', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax.spines['top'].set_visible(False)\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['bottom'].set_visible(True)\n",
        "ax.spines['left'].set_visible(False)\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CvOFrXGXwXrV",
        "cellView": "form"
      },
      "source": [
        "#@title **Compute radius of gyration**\n",
        "\n",
        "Selection = \"Protein\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "\n",
        "Output_name = 'radius_gyration' #@param {type:\"string\"}\n",
        "\n",
        "radgyr = pt.radgyr(traj_load, mask = mask_end)\n",
        "time = len(rmsd)*int(Write_the_trajectory)/1000\n",
        "time_array = np.arange(0,time,int(Write_the_trajectory)/1000)*int(stride_traj)\n",
        "\n",
        "# Plotting:\n",
        "plt.plot(time_array, radgyr, alpha=0.6, color = 'green', linewidth = 1.0)\n",
        "plt.xlim(0, simulation_ns)\n",
        "#plt.ylim(2, 6)\n",
        "\n",
        "plt.xlabel(\"Time (ns)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"Radius of gyration ($\\AA$)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(radgyr)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2Q7FKg8Fuxr9",
        "cellView": "form"
      },
      "source": [
        "#@title **Plot radius of gyration as a distribution**\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'radius_gyration_dist' #@param {type:\"string\"}\n",
        "\n",
        "ax = sb.kdeplot(radgyr, color=\"green\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('Radius of gyration ($\\AA$)', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax.spines['top'].set_visible(False)\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['bottom'].set_visible(True)\n",
        "ax.spines['left'].set_visible(False)\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p2Y0DgwTxLWc",
        "cellView": "form"
      },
      "source": [
        "#@title **Compute RMSF**\n",
        "\n",
        "Selection = \"Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "  axis= \"Residue\"\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "  axis= \"Atom\"\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "  axis= \"Atom\"\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "\n",
        "Output_name = 'rmsf' #@param {type:\"string\"}\n",
        "\n",
        "traj_align = pt.align(traj_load, mask=mask_end, ref=0)\n",
        "rmsf = pt.rmsf(traj_align, mask = mask_end)\n",
        "# bfactor = pt.bfactors(traj_load, byres=False)\n",
        "\n",
        "# Plotting:\n",
        "plt.plot(rmsf[:,1], alpha=1.0, color = 'red', linewidth = 1.0)\n",
        "\n",
        "plt.xlabel(axis, fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel(\"RMSF ($\\AA$)\", fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.xlim(0, len(rmsf[:-1]))\n",
        "\n",
        "#plt.xticks(np.arange(min(rmsf[:1]), max(rmsf[:1])))\n",
        "plt.yticks(fontsize = 12)\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(rmsf)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JalicqqrTodW",
        "cellView": "form"
      },
      "source": [
        "#@title **2D RMSD**\n",
        "\n",
        "Selection = \"Protein+Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "\n",
        "Output_name = '2D_rmsd' #@param {type:\"string\"}\n",
        "\n",
        "last_frame = len(time_array)\n",
        "\n",
        "stride_ticks_f = (last_frame)/5\n",
        "ticks_frame = np.arange(0,(len(time_array) + float(stride_ticks_f)), float(stride_ticks_f))\n",
        "a = ticks_frame.astype(float)\n",
        "stride_ticks_t = (simulation_ns)/5\n",
        "tick_time = np.arange(0,(float(simulation_ns) + float(stride_ticks_t)), float(stride_ticks_t))\n",
        "b = tick_time.astype(float)\n",
        "\n",
        "mat1 = pt.pairwise_rmsd(traj_load, mask=mask_end, frame_indices=range(int(number_frames_analysis)))\n",
        "\n",
        "\n",
        "ax = plt.imshow(mat1, cmap = 'PRGn', origin='lower', interpolation = 'bicubic')\n",
        "plt.title('2D RMSD')\n",
        "plt.xlabel('Time (ns)', fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel('Time (ns)', fontsize = 14, fontweight = 'bold')\n",
        "# plt.xticks(fontsize = 12)\n",
        "# plt.yticks(fontsize = 12)\n",
        "plt.xticks(a, b.round(decimals=3), fontsize = 12)\n",
        "plt.yticks(a, b.round(decimals=3), fontsize = 12)\n",
        "# plt.xlim(0, a[-1])\n",
        "# plt.ylim(0, a[-1])\n",
        "\n",
        "cbar1 = plt.colorbar()\n",
        "cbar1.set_label(\"RMSD ($\\AA$)\", fontsize = 14, fontweight = 'bold')\n",
        "\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(mat1)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_mgVSbBshWFV",
        "cellView": "form"
      },
      "source": [
        "#@title **Calculate eigenvectors of Principle Component Analysis (PCA)**\n",
        "Selection = \"Protein+Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "\n",
        "data = pt.pca(traj_load, fit=True, ref=0, mask=mask_end, n_vecs=2)\n",
        "#print('projection values of each frame to first mode = {} \\n'.format(data[0][0]))\n",
        "#print('projection values of each frame to second mode = {} \\n'.format(data[0][1]))\n",
        "#print('eigvenvalues of first two modes', data[1][0])\n",
        "#print(\"\")\n",
        "#print('eigvenvectors of first two modes: \\n', data[1][1])\n",
        "\n",
        "last_frame = len(time_array)\n",
        "\n",
        "stride_ticks_f = (last_frame)/5\n",
        "ticks_frame = np.arange(0,(len(time_array) + float(stride_ticks_f)), float(stride_ticks_f))\n",
        "a = ticks_frame.astype(float)\n",
        "a2 = a.tolist()\n",
        "stride_ticks_t = (simulation_ns)/5\n",
        "tick_time = np.arange(0,(float(simulation_ns) + float(stride_ticks_t)), float(stride_ticks_t))\n",
        "b = tick_time.astype(float)\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "Output_name = 'PCA' #@param {type:\"string\"}\n",
        "\n",
        "Output_PC1 = 'PC1' #@param {type:\"string\"}\n",
        "Output_PC2 = 'PC2' #@param {type:\"string\"}\n",
        "\n",
        "%matplotlib inline\n",
        "%config InlineBackend.figure_format = 'retina'  # high resolution\n",
        "projection_data = data[0]\n",
        "plt.title(r'PCA')\n",
        "PC1 = data[0][0] \n",
        "PC2 = data[0][1]\n",
        "\n",
        "a = plt.scatter(PC1,PC2, c=range(int(number_frames_analysis)), cmap='Greens', marker='o',s=8, alpha=1)\n",
        "plt.clim(0, last_frame)\n",
        "\n",
        "plt.xlabel('PC1', fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel('PC2', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "# N = len(number_frames)\n",
        "# x2 = np.arange(N)\n",
        "\n",
        "cbar1 = plt.colorbar(a, orientation=\"vertical\")\n",
        "cbar1.set_label('Time(ns)', fontsize = 14, fontweight = 'bold')\n",
        "cbar1.set_ticks(a2)\n",
        "cbar1.set_ticklabels(b.round(decimals=3))\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "pc1=pd.DataFrame(PC1)\n",
        "pc1.to_csv(os.path.join(workDir, Output_PC1 + \".csv\"))\n",
        "pc2=pd.DataFrame(PC2)\n",
        "pc2.to_csv(os.path.join(workDir, Output_PC2 + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yce9RfNtpl-J",
        "cellView": "form"
      },
      "source": [
        "#@title **Plot Principal Component 1 (PC1) and Principal Component 2 (PC2) as a distribution**\n",
        "Output_name = 'PCA_dist' #@param {type:\"string\"}\n",
        "\n",
        "\n",
        "fig = plt.figure(figsize=(9,5))\n",
        "\n",
        "plt.subplot(1, 2, 1)\n",
        "ax = sb.kdeplot(PC1, color=\"green\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('PC1', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax.spines['top'].set_visible(False)\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['bottom'].set_visible(True)\n",
        "ax.spines['left'].set_visible(False)\n",
        "plt.subplot(1, 2, 2)\n",
        "ax2 = sb.kdeplot(PC2, color=\"purple\", shade=True, alpha=0.2, linewidth=0.5)\n",
        "plt.xlabel('PC2', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks([])\n",
        "plt.ylabel('')\n",
        "ax2.spines['top'].set_visible(False)\n",
        "ax2.spines['right'].set_visible(False)\n",
        "ax2.spines['bottom'].set_visible(True)\n",
        "ax2.spines['left'].set_visible(False)\n",
        "\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pTDb7CEfkLq1",
        "cellView": "form"
      },
      "source": [
        "#@title **Pearson's Cross Correlation (CC)**\n",
        "\n",
        "Selection = \"Protein+Glycan\" #@param [\"Protein\", \"Protein+Glycan\", \"Glycan\"]\n",
        "\n",
        "if Selection == \"Protein\":\n",
        "  mask_end = mask_protein\n",
        "elif Selection == \"Protein+Glycan\":\n",
        "  mask_end = mask_protein_glycan\n",
        "else:\n",
        "  mask_end = mask_glycan\n",
        "\n",
        "#@markdown **Provide output file names below:** \n",
        "\n",
        "Output_name = 'cross_correlation' #@param {type:\"string\"}\n",
        "\n",
        "\n",
        "traj_align = pt.align(traj_load, mask=mask_end, ref=0)\n",
        "\n",
        "mat_cc = matrix.correl(traj_align, mask_end)\n",
        "\n",
        "ax = plt.imshow(mat_cc, cmap = 'PiYG_r', interpolation = 'bicubic', vmin = -1, vmax = 1, origin='lower')\n",
        "\n",
        "plt.xlabel('Residues', fontsize = 14, fontweight = 'bold')\n",
        "plt.ylabel('Residues', fontsize = 14, fontweight = 'bold')\n",
        "plt.xticks(fontsize = 12)\n",
        "plt.yticks(fontsize = 12)\n",
        "cbar1 = plt.colorbar()\n",
        "cbar1.set_label('$CC_ij$', fontsize = 14, fontweight = 'bold')\n",
        "\n",
        "plt.savefig(os.path.join(workDir, Output_name + \".png\"), dpi=600, bbox_inches='tight')\n",
        "\n",
        "raw_data=pd.DataFrame(mat_cc)\n",
        "raw_data.to_csv(os.path.join(workDir, Output_name + \".csv\"))"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}