Trajectory_analysis_tools.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [
        "vv5cq2Wv7UhB",
        "lZM0vkKdy6v6",
        "reax83yQAG-9"
      ],
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#### Trajectory analysis tools <font color=\"DarkSeaGreen\">/ GROMACS-on-Colab</font> [<img src=\"https://github.githubassets.com/favicons/favicon.png\" width=\"16\">](https://github.com/bioinfkaustin/gromacs-on-colab)\n",
        "<small>Suite: [`Build_to_Google_Drive.ipynb`](https://colab.research.google.com/github/bioinfkaustin/gromacs-on-colab/blob/main/Build_to_Google_Drive.ipynb) | [`GROMACS_for_CHARMM-GUI.ipynb`](https://colab.research.google.com/github/bioinfkaustin/gromacs-on-colab/blob/main/GROMACS_for_CHARMM-GUI.ipynb) | [`GROMACS_for_production.ipynb`](https://colab.research.google.com/github/bioinfkaustin/gromacs-on-colab/blob/main/GROMACS_for_production.ipynb) | `Trajectory_analysis_tools.ipynb`</small>"
      ],
      "metadata": {
        "id": "4zwSgdGDP9cV"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Documentation\n",
        "Please click ***↳ cells hidden*** below to show the documentation for this notebook."
      ],
      "metadata": {
        "id": "vv5cq2Wv7UhB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### License\n",
        "\n",
        "> This notebook as a work of software is licensed under the terms of the [AGPL-3.0](https://opensource.org/licenses/AGPL-3.0) or later."
      ],
      "metadata": {
        "id": "5dTAvj38Dal8"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### About this software\n",
        "\n",
        "> This notebook allows a number of **analysis tools** to be run on a simulated trajectory. A description of each tool is provided in its respective code cell.\n",
        ">\n",
        "> The specified **project folder** must contain the final output of a production simulation (e.g. `sim_reference.xtc`), on which the analyses will be performed."
      ],
      "metadata": {
        "id": "Jb0ESikELkRH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### Installation\n",
        "\n",
        "> The installation notebook, [`Build_to_Google_Drive.ipynb`](https://colab.research.google.com/github/bioinfkaustin/gromacs-on-colab/blob/main/Build_to_Google_Drive.ipynb), must be run before using this notebook."
      ],
      "metadata": {
        "id": "Cl4YtmH7mICa"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Configuration"
      ],
      "metadata": {
        "id": "1pWciuVhQCYo"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Pm9Dk5FBDoF4",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "\n",
        "import os\n",
        "import re\n",
        "import shutil\n",
        "\n",
        "#@markdown Specify the location of the **GROMACS project folder** to analyse. It should contain the outputs from a production run, e.g. `sim.gro`, `sim.trr` and/or `sim.xtc`, and so on.\n",
        "project_folder = \"{GoogleDrive}/GROMACS/7FBF_FABPH_vs_octanoic_acid\" #@param {type: \"string\"}\n",
        "# default: {GoogleDrive}/GROMACS/7FBF_FABPH_vs_octanoic_acid\n",
        "\n",
        "#@markdown Provide the unique filename prefix used for the simulation.\n",
        "output_prefix = \"sim\" #@param {type: \"string\"}\n",
        "# default: sim\n",
        "\n",
        "#@markdown **After filling in this form, run the notebook by clicking *Runtime -> Run all* in the toolbar.**\n",
        "\n",
        "\n",
        "#\n",
        "# Google Drive\n",
        "#\n",
        "\n",
        "if not os.path.isdir(\"/content/drive/MyDrive\") and project_folder.startswith(\"{GoogleDrive}\"):\n",
        "  from google.colab import drive\n",
        "  drive.mount(\"/content/drive\")\n",
        "\n",
        "if not os.path.isdir(\"/content/drive/MyDrive\"):\n",
        "  raise RuntimeError(\"Error: could not connect to Google Drive\")\n",
        "\n",
        "\n",
        "#\n",
        "# Methods for parsing and validation\n",
        "#\n",
        "\n",
        "def _path(s, exists=False):\n",
        "  if \"{GoogleDrive}\" in s and not s.startswith(\"{GoogleDrive}\"):  raise ValueError(f\"Error: {{GoogleDrive}} is a path prefix, but appears later: {s}\")\n",
        "  s = s.format(GoogleDrive=\"/content/drive/MyDrive\")\n",
        "  #     ^^^ raises KeyError if any {...} placeholder is present except {GoogleDrive}\n",
        "\n",
        "  if exists  and not os.path.isdir(s):  raise FileNotFoundError(f\"Error: folder not found: {s}\")\n",
        "  return os.path.abspath(s)\n",
        "\n",
        "def parse(s, mandatory=False, path=False, exists=False):\n",
        "  s = s.strip()\n",
        "  # error on e.g. \"\"\n",
        "  if mandatory and not s:\n",
        "    raise ValueError(\"Error: mandatory setting without value\")\n",
        "  # e.g. {GoogleDrive}/archive.tgz -> /content/drive/MyDrive/archive.tgz\n",
        "  if path:\n",
        "    s = _path(s, exists=exists)\n",
        "  return s\n",
        "\n",
        "\n",
        "#\n",
        "# Validate the input values\n",
        "#\n",
        "\n",
        "# project_folder\n",
        "\n",
        "project_folder = parse(project_folder, mandatory=True, path=True, exists=True)\n",
        "\n",
        "\n",
        "# output_prefix\n",
        "\n",
        "if not output_prefix:\n",
        "  raise RuntimeError(\"Error: an output prefix must be provided\")\n",
        "\n",
        "if not re.match(r\"^[0-9a-zA-Z]+$\", output_prefix):\n",
        "  raise RuntimeError(f\"Error: special characters are not allowed in output prefix, but got: {output_prefix}\")\n",
        "\n",
        "\n",
        "#\n",
        "# Make sure that the notebook is in the start folder\n",
        "#\n",
        "\n",
        "if \"START\" not in os.environ or not os.environ[\"START\"]:\n",
        "  %env START={os.getcwd()}\n",
        "else:\n",
        "  %cd {os.environ[\"START\"]}\n",
        "\n",
        "\n",
        "#\n",
        "# Use a clean scratch directory for the rest of the run\n",
        "#\n",
        "\n",
        "try:\n",
        "  shutil.rmtree(\"scratch\")\n",
        "except FileNotFoundError:\n",
        "  pass\n",
        "os.makedirs(\"scratch\")\n",
        "%cd \"scratch\"\n",
        "\n",
        "\n",
        "#\n",
        "# Verbosity\n",
        "#\n",
        "\n",
        "#@markdown \\\n",
        "#@markdown By default, only relevant output is displayed. To see the full output, including any error messages, enable the `verbose` option below.\n",
        "verbose=False #@param {type: \"boolean\"}\n",
        "if verbose:\n",
        "  %env VERBOSE=true\n",
        "else:\n",
        "  %env VERBOSE=false"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Installation"
      ],
      "metadata": {
        "id": "lZM0vkKdy6v6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@markdown In the following cells, applications are downloaded from a **persistent cache** in your Google Drive.\n",
        "#@markdown\n",
        "#@markdown This cell sets up the cache folder.\n",
        "\n",
        "storage = \"/content/drive/MyDrive/gromacs-on-colab\"\n",
        "%env STORAGE={storage}"
      ],
      "metadata": {
        "cellView": "form",
        "id": "bhcpS51_dHd4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash\n",
        "#@markdown **GROMACS** is installed from cache.\n",
        "\n",
        "if [[ -d \"/usr/local/gromacs\" ]]; then\n",
        "  exit 0  # already installed\n",
        "fi\n",
        "\n",
        "gromacs_vers=\"2023.2\" #@param {type: \"string\"}\n",
        "cache_gromacs=\"${STORAGE}/gromacs-${gromacs_vers}.tar.gz\"\n",
        "\n",
        "if [[ -s \"${cache_gromacs}\" ]]; then\n",
        "  tar -xzf \"${cache_gromacs}\" -C \"/usr/local\"\n",
        "else\n",
        "  echo \"Error: GROMACS installation not found\" >&2\n",
        "  echo \"(Have you installed GROMACS to your Google Drive?)\" >&2\n",
        "  exit 1\n",
        "fi"
      ],
      "metadata": {
        "cellView": "form",
        "id": "WV8aP1j0Ra0E"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash\n",
        "#@markdown The CHARMM36 forcefield is downloaded from cache.\n",
        "\n",
        "if [[ -d \"${START}/charmm36.ff\" ]]; then\n",
        "  exit 0  # already installed\n",
        "fi\n",
        "\n",
        "charmm36_vers=\"jul2022\" #@param {type: \"string\"}\n",
        "cache_charmm36=\"${STORAGE}/charmm36-${charmm36_vers}.tar.gz\"\n",
        "\n",
        "if [[ -s \"${cache_charmm36}\" ]]; then\n",
        "  tar -xzf \"${cache_charmm36}\" -C \"${START}\"\n",
        "else\n",
        "  echo \"Error: CHARMM36 forcefield installation not found\" >&2\n",
        "  echo \"(Have you installed the CHARMM36 forcefield to your Google Drive?)\" >&2\n",
        "  exit 1\n",
        "fi"
      ],
      "metadata": {
        "cellView": "form",
        "id": "QWr9d9Pg5-R5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Input"
      ],
      "metadata": {
        "id": "reax83yQAG-9"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$project_folder\" \"$output_prefix\"\n",
        "project_folder=\"$1\"\n",
        "output_prefix=\"$2\"\n",
        "#@markdown Extract the trajectory from the project folder.\n",
        "\n",
        "if [[ ! -z \"$(ls -A)\" ]]; then\n",
        "  exit 0  # already extracted\n",
        "fi\n",
        "\n",
        "if [[ ! -d \"${project_folder}\" ]]; then\n",
        "  echo \"Error: folder not found: ${project_folder}\" >&2\n",
        "  exit 1\n",
        "fi\n",
        "\n",
        "pushd \"${project_folder}\"\n",
        "\n",
        "cp \"grompp.mdp\"  \"conf.gro\"  \"restraint.gro\"  \"index.ndx\"  \"topol.top\"  \"${START}/scratch/\" 2> /dev/null\n",
        "\n",
        "top_level_dir=\"$(realpath .)\"\n",
        "function get_includes_recursively {\n",
        "  f=\"$1\"\n",
        "  sed -E \"/^#include/!d; s/^#include +['\\\"]//; s/['\\\"]$//\" \"${f}\" | while read -r g; do\n",
        "    d=\"$(dirname \"${g}\")\"\n",
        "    b=\"$(basename \"${g}\")\"\n",
        "    pushd \"${d}\" > /dev/null\n",
        "    get_includes_recursively \"${b}\"\n",
        "    echo \"$(realpath --relative-to=\"${top_level_dir}\" \"${b}\")\"\n",
        "    popd > /dev/null\n",
        "  done\n",
        "}\n",
        "get_includes_recursively \"topol.top\" | while read -r f; do\n",
        "  cp --parents \"${f}\" \"${START}/scratch/\"\n",
        "done\n",
        "\n",
        "if [[ -s \"${output_prefix}.part0001.tar.gz\" ]]; then  # part-file data\n",
        "  latest_part=\"$(ls | egrep '^'\"${output_prefix}\"'\\.part[0-9]+\\.tar.gz$' | sort -V | tail -n1 | egrep -o 'part[0-9]+')\"\n",
        "  tar -xzf \"${output_prefix}.${latest_part}.tar.gz\" -C \"${START}/scratch\" 2> /dev/null\n",
        "else  # backwards compatible with appended-file data from previous versions of the notebook\n",
        "  cp \"${output_prefix}\".* \"${START}/scratch/\" 2> /dev/null\n",
        "fi\n",
        "\n",
        "cp \"${output_prefix}_reference.xtc\" \"${START}/scratch/\" 2> /dev/null\n",
        "\n",
        "popd\n",
        "\n",
        "if [[ ! -s \"grompp.mdp\" || ! -s \"conf.gro\" || ! -s \"topol.top\" ]]; then\n",
        "  echo \"Error: essential files not found: grompp.mdp, conf.gro, topol.top\" >&2\n",
        "  exit 1\n",
        "fi\n",
        "\n",
        "if [[ ! -s \"${output_prefix}.tpr\" ]]; then\n",
        "  echo \"Error: simulation parameters not found: ${output_prefix}.tpr\" >&2\n",
        "  exit 1\n",
        "fi\n",
        "\n",
        "# The output trajectory is called ..._reference.xtc and must be present.\n",
        "# For backwards compatibility, it can be made if missing, but only if the data is in the appended-file format used by previous versions of this toolkit.\n",
        "if [[ ! -s \"${output_prefix}_reference.xtc\" ]]; then\n",
        "  source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "  [[ -s \"${output_prefix}.trr\" ]] && trr_frames=\"$(gmx check -f \"${output_prefix}.trr\" 2>&1 | awk '$1 ~ /^Coords/ { print $2 }')\" || trr_frames=0\n",
        "  [[ -s \"${output_prefix}.xtc\" ]] && xtc_frames=\"$(gmx check -f \"${output_prefix}.xtc\" 2>&1 | awk '$1 ~ /^Coords/ { print $2 }')\" || xtc_frames=0\n",
        "  if (( $trr_frames == 0 && $xtc_frames == 0 )); then  # neither found, therefore error in part-file data\n",
        "    echo \"Error: output trajectory not found: ${output_prefix}_reference.xtc\" >&2\n",
        "    exit 1\n",
        "  elif (( $trr_frames >= $xtc_frames )); then  # this is appended-file data and the trr has more frames, or same frames with better precision\n",
        "    xtcext=\"trr\"\n",
        "    ndecarg=\"-ndec 4\"\n",
        "  else  # this is appended-file data and the xtc has more frames\n",
        "    xtcext=\"xtc\"\n",
        "    ndecarg=\"\"\n",
        "  fi\n",
        "  gmx trjconv -f \"${output_prefix}.${xtcext}\" -s \"${output_prefix}.tpr\" -pbc nojump $ndecarg -o \"nojump.xtc\" <<< \"0\"\n",
        "  gmx trjconv -f \"nojump.xtc\" -s \"${output_prefix}.tpr\" -fit progressive $ndecarg -o \"progressive.xtc\" < <(echo \"C-alpha\"; echo \"0\")\n",
        "  gmx trjconv -f \"progressive.xtc\" -s \"${output_prefix}.tpr\" -pbc mol $ndecarg -o \"${output_prefix}_reference.xtc\" <<< \"0\"\n",
        "fi"
      ],
      "metadata": {
        "cellView": "form",
        "id": "KhWcXUPPzoNk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Analyses"
      ],
      "metadata": {
        "id": "RjORKXToACjc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@markdown This cell deliberately throws an error, to interrupt *Run all*.\n",
        "#@markdown\n",
        "#@markdown **In the cells below, press the <u>play button</u> (left) to use the specified tool.**\n",
        "raise StopIteration(\"Please continue by running the desired cells manually.\")\n",
        "\n"
      ],
      "metadata": {
        "cellView": "form",
        "id": "nwaG8Uh6DO0R"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "oyzcEN77agrl"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$project_folder\"\n",
        "project_folder=\"$1\"\n",
        "\n",
        "#@markdown **Save analyses to Google Drive.** After using the tools below, run this cell to save outputs (e.g. graphs) to the project folder.\n",
        "#@markdown\n",
        "#@markdown (The upload may take a moment. Please verify the files are in your Google Drive before disconnecting from the session.)\n",
        "#@markdown\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* `dat.*`</small>\n",
        "\n",
        "cp \"dat.\"* \"${project_folder}/\""
      ],
      "metadata": {
        "cellView": "form",
        "id": "u6HS-fy8a2bq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "wDjIJRSIa1Gq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$output_prefix\"\n",
        "exec 3>&1\n",
        "if $VERBOSE; then\n",
        "  exec 2>&1  # fix stderr and stdout appearing out of order\n",
        "else\n",
        "  exec >/dev/null 2>/dev/null  # silence\n",
        "fi\n",
        "output_prefix=\"$1\"\n",
        "\n",
        "#@markdown **Compute a centroid,** a representative frame within a time range.\n",
        "#@markdown\n",
        "#@markdown The frames in the time range are clustered by the RMSD of the specified group. In the largest cluster, the frame with the smallest mean RMSD to the other frames is the \"centroid.\"\n",
        "rmsd_group=\"LIG\" #@param {type: \"string\"}\n",
        "start_ps=0 #@param {type: \"integer\"}\n",
        "end_ps=1000 #@param {type: \"integer\"}\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* optionally, `index.ndx` \\\n",
        "#@markdown *Output:* `dat.{output_prefix}_centroid_{rmsd_group}_{start_ps}-{end_ps}.gro` and the equivalent `.pdb`</small>\n",
        "\n",
        "source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "[[ -s \"index.ndx\" ]] && ndxarg=\"-n index.ndx\" || ndxarg=\"\"\n",
        "\n",
        "\n",
        "# Calculate the all-against-all matrix of RMSDs for the time range\n",
        "( echo \"C-alpha\"; echo \"${rmsd_group}\" ) | gmx rms -nomw -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" -b $start_ps -e $end_ps $ndxarg -m \"rmsd.xpm\"\n",
        "\n",
        "# Use the RMSDs to cluster the frames in the time range\n",
        "( echo \"System\"; echo \"System\" ) | gmx cluster -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" -b $start_ps -e $end_ps -dm \"rmsd.xpm\" -cl \"clusters.gro\" -clid \"clust-id.xvg\" -sz \"clust-size.xvg\"\n",
        "\n",
        "# Get the centroid of the largest cluster\n",
        "target_frame=$(paste <(fgrep \"frame t=\" \"clusters.gro\") <(sed '/^[@#]/d' \"clust-size.xvg\") | sort -n -k5 | tail -n1 | awk '{ print $3 }')\n",
        "name=\"${output_prefix}_centroid_${rmsd_group}_${start_ps}-${end_ps}\"\n",
        "gmx trjconv -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" -dump $target_frame -o \"dat.${name}.gro\" <<< \"0\"\n",
        "\n",
        "gmx trjconv -s \"${output_prefix}.tpr\" -f \"dat.${name}.gro\" -o \"dat.${name}.pdb\" <<< \"0\"\n",
        "\n",
        "rm \"rmsd.xpm\" \"clusters.gro\" \"clust-id.xvg\" \"clust-size.xvg\"\n",
        "\n",
        "for i in $(seq 5); do echo \"\"; done\n",
        "{\n",
        "  echo \"***\"\n",
        "  echo \"Made:  dat.${name}.gro\"\n",
        "  echo \"Made:  dat.${name}.pdb\"\n",
        "  echo \"***\"\n",
        "} >&3"
      ],
      "metadata": {
        "id": "xXKZaZh068o3",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "Xj4_U0QkanKh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$output_prefix\"\n",
        "output_prefix=\"$1\"\n",
        "exec 3>&1\n",
        "if $VERBOSE; then\n",
        "  exec 2>&1  # fix stderr and stdout appearing out of order\n",
        "else\n",
        "  exec >/dev/null 2>/dev/null  # silence\n",
        "fi\n",
        "\n",
        "#@markdown **Compute RMSD over time data** for the trajectory, with respect to some specified frame of reference.\n",
        "#@markdown\n",
        "#@markdown A frame of reference can either be a **time** in ps, e.g. 0, or a **time range of a computed centroid**, e.g. 0-1000.\n",
        "#@markdown\n",
        "#@markdown (For convenience, you may specify multiple RMSD groups and/or frames of reference, comma separated, and all will be calculated.)\n",
        "rmsd_group=\"LIG\" #@param {type: \"string\"}\n",
        "reference_ps=\"0, 0-1000\" #@param {type: \"string\"}\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* optionally, `index.ndx`; and optionally, for centroid input, `dat.{output_prefix}_centroid_{rmsd_group}_{start_ps}-{end_ps}.gro` \\\\\n",
        "#@markdown *Output:* `dat.{output_prefix}_rmsd_{rmsd_group}_{reference_ps}.csv`</small>\n",
        "\n",
        "source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "[[ -s \"index.ndx\" ]] && ndxarg=\"-n index.ndx\" || ndxarg=\"\"\n",
        "\n",
        "log=\"$(mktemp)\"\n",
        "\n",
        "rmsd_groups=\"$(tr -d \" \" <<< \"${rmsd_group}\" | tr \",\" \" \")\"\n",
        "for rmsd_group in $rmsd_groups; do\n",
        "\n",
        "  ranges_ps=\"$(tr -d \" \" <<< \"${reference_ps}\" | tr \",\" \" \")\"\n",
        "  for range_ps in $ranges_ps; do\n",
        "\n",
        "    if [[ -z \"${range_ps}\" || $range_ps == 0 ]]; then  # first frame\n",
        "      if [[ -s \"conf.gro\" ]]; then\n",
        "        input=\"conf.gro\"\n",
        "      else\n",
        "        range_ps=0  # extract the first frame later\n",
        "      fi\n",
        "    fi\n",
        "\n",
        "    if [[ -z \"${range_ps}\" ]]; then  # already handled\n",
        "      :\n",
        "    elif [[ \"${range_ps}\" =~ - ]]; then  # a centroid already calculated\n",
        "      input=\"dat.${output_prefix}_centroid_${rmsd_group}_${range_ps}.gro\"\n",
        "      if [[ ! -s \"${input}\" ]]; then\n",
        "        echo \"Error: centroid not found for range ${range_ps}: ${input}\" >&2\n",
        "        exit 1\n",
        "      fi\n",
        "    else  # some specified frame to extract\n",
        "      frame_ps=$(($range_ps))\n",
        "      if [[ -z \"${frame_ps}\" ]]; then\n",
        "        echo \"Error: invalid frame time: ${frame_ps}\" >&2\n",
        "        exit 1\n",
        "      fi\n",
        "      input=\"dat.${output_prefix}_frame_${frame_ps}.gro\"\n",
        "      if [[ ! -s \"${input}\" ]]; then\n",
        "        gmx trjconv -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" -dump $frame_ps -o \"${input}\" <<< \"0\"\n",
        "\n",
        "        echo \"Made:  ${input}\" >> \"${log}\"\n",
        "      fi\n",
        "    fi\n",
        "\n",
        "    ( echo \"C-alpha\"; echo \"${rmsd_group}\"; ) | gmx rms -nomw -s \"${input}\" $ndxarg -f \"${output_prefix}_reference.xtc\" -o \"rmsd.xvg\"\n",
        "\n",
        "    name=\"${output_prefix}_rmsd_${rmsd_group}_${range_ps}\"\n",
        "    {\n",
        "      echo \"Time_(ns),RMSD_(A)\"\n",
        "      sed '/^[@#]/d' \"rmsd.xvg\" | awk '{ print $1/1000\",\"$2*10 }'\n",
        "    } > \"dat.${name}.csv\"\n",
        "\n",
        "    echo \"Made:  dat.${name}.csv\" >> \"${log}\"\n",
        "\n",
        "    rm \"rmsd.xvg\"\n",
        "\n",
        "  done\n",
        "\n",
        "done\n",
        "\n",
        "for i in $(seq 5); do echo \"\"; done\n",
        "{\n",
        "  echo \"***\"\n",
        "  cat \"${log}\"\n",
        "  echo \"***\"\n",
        "} >&3\n",
        "rm \"${log}\""
      ],
      "metadata": {
        "id": "9NvIMpWA6Wfd",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "wS5lr8gdaoXE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@markdown **Plot computed RMSD data.** Tabular output from an RMSD analysis (above) may be visualised as a plot.\n",
        "#@markdown\n",
        "#@markdown Multiple RMSD groups and/or frames of reference may be specified as a comma separated list, in which case multiple plots will be made on the same axes.\n",
        "rmsd_groups = \"LIG\" #@param {type: \"string\"}\n",
        "reference_ps = \"0, 0-1000\" #@param {type: \"string\"}\n",
        "#@markdown Optionally, plot a threshold line (on the $y$-axis) at a particular RMSD value. (0 to disable.)\n",
        "threshold_A = 5.5 #@param {type: \"number\"}\n",
        "#@markdown Optionally, cut off the plot at a particular maximum time. (0 to disable.)\n",
        "max_t_ns = 0 #@param {type: \"number\"}\n",
        "# defaults:    \"LIG\"    \"0, 0-1000\"    5.5    0\n",
        "# \"\"\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* `dat.{output_prefix}_rmsd_{rmsd_group}_{reference_ps}.csv` \\\\\n",
        "#@markdown *Output:* `dat.{output_prefix}_rmsd_{rmsd_groups}_{ranges_ps}.png`</small>\n",
        "\n",
        "rmsd_groups = rmsd_groups.replace(\" \", \"\").split(\",\")\n",
        "ranges_ps = reference_ps.replace(\" \", \"\").split(\",\")\n",
        "\n",
        "\n",
        "#\n",
        "# Get the data from the various RMSD files\n",
        "#\n",
        "\n",
        "import pandas as pd\n",
        "\n",
        "t_ = \"Time_(ns)\"\n",
        "y_ = \"RMSD_(A)\"\n",
        "\n",
        "groups, ranges, labels, tables = list(), list(), list(), list()\n",
        "for g in rmsd_groups:\n",
        "  for r in ranges_ps:\n",
        "    try:\n",
        "      df = pd.read_csv(f\"dat.{output_prefix}_rmsd_{g}_{r}.csv\")\n",
        "    except:\n",
        "      continue\n",
        "\n",
        "    groups.append(g)\n",
        "    ranges.append(r)\n",
        "    if \"-\" in r:\n",
        "      start = float(r.split(\"-\")[0]) / 1000.\n",
        "      labels.append(f\"vs. {r} ns centroid: {g}\")\n",
        "    else:\n",
        "      start = float(r) / 1000.\n",
        "      labels.append(f\"vs. {r} ns frame: {g}\")\n",
        "\n",
        "    # Apply $t$ range\n",
        "    df = df[(df[t_] >= start)]\n",
        "    if max_t_ns > 0.:\n",
        "      df = df[(df[t_] <= max_t_ns)]\n",
        "    tables.append(df)\n",
        "\n",
        "\n",
        "#\n",
        "# Plot\n",
        "#\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "from matplotlib.ticker import AutoLocator\n",
        "\n",
        "def moving_average(df):\n",
        "  return df.rolling(21, center=True).mean()\n",
        "\n",
        "plt.figure()\n",
        "\n",
        "with sns.color_palette(\"deep\", n_colors=len(tables)):\n",
        "  for l, df in zip(labels, tables):\n",
        "    plt.plot(df[t_], df[y_], alpha=0.1)\n",
        "\n",
        "  for l, df in zip(labels, tables):  # loop again to use same colours\n",
        "    smooth = moving_average(df)\n",
        "    plt.plot(smooth[t_], smooth[y_], label=l)\n",
        "\n",
        "  plt.legend()\n",
        "\n",
        "  if threshold_A != 0.:\n",
        "    plt.axhline(y=threshold_A, color=\"black\", linestyle=\":\")\n",
        "\n",
        "# Title and axis labels\n",
        "plt.title(os.path.basename(project_folder).replace(\"_\", \" \").replace(\"vs \", \"vs. \"))\n",
        "plt.xlabel(\"Time (ns)\")\n",
        "plt.ylabel(f\"RMSD of group (A)\")\n",
        "\n",
        "# Set x-axis ticks\n",
        "ticker = AutoLocator()\n",
        "ticker.set_params(steps=[1, 2, 3, 5, 10] if round(df[t_].values[-1]) % 3 == 0 else [1, 2, 5, 10], min_n_ticks=6)\n",
        "ax = plt.gca()\n",
        "ax.xaxis.set_major_locator(ticker)\n",
        "\n",
        "# Set y-axis limit\n",
        "if any(max(df[y_]) > 16. for df in tables):\n",
        "  ax.set_ylim([0, 16.])\n",
        "\n",
        "\n",
        "#\n",
        "# Output\n",
        "#\n",
        "\n",
        "groups_name = \"_\".join(rmsd_groups)\n",
        "ranges_name = \"_\".join(ranges_ps)\n",
        "name = f\"{output_prefix}_rmsd_{groups_name}_{ranges_name}\"\n",
        "plt.gcf().savefig(f\"dat.{name}.png\")\n",
        "plt.show()\n",
        "\n",
        "print(\"***\")\n",
        "print(f\"Made:  dat.{name}.png\")\n",
        "print(\"***\")"
      ],
      "metadata": {
        "id": "QIVKLd7EAdve",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "7IE7dNDyapkE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$output_prefix\"\n",
        "exec 3>&1\n",
        "if $VERBOSE; then\n",
        "  exec 2>&1  # fix stderr and stdout appearing out of order\n",
        "else\n",
        "  exec >/dev/null 2>/dev/null  # silence\n",
        "fi\n",
        "output_prefix=\"$1\"\n",
        "\n",
        "#@markdown **Compute atom pair distance over time data** for the trajectory.\n",
        "\n",
        "#@markdown The atom numbers (IDs) of the atoms to measure the distance between. These are found in e.g. the `conf.gro` file for the system.\n",
        "pair_label=\"A-B\" #@param {type: \"string\"}\n",
        "atom_a=1234 #@param {type: \"integer\"}\n",
        "atom_b=5678 #@param {type: \"integer\"}\n",
        "\n",
        "#@markdown (Please note that GROMACS and VMD use atom numbers which are different by $1$. If you use the VMD picker tool to inspect an atom and it says its number is $1000$, please type $1001$ above.)\n",
        "#@markdown\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Output:* `dat.{output_prefix}_distance_{pair_label}_0.csv`</small>\n",
        "\n",
        "source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "\n",
        "cat > \"dist.ndx\" <<EOS\n",
        "[ Distance ]\n",
        "${atom_a} ${atom_b}\n",
        "EOS\n",
        "\n",
        "gmx distance -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" -n \"dist.ndx\" -oall \"dist.xvg\" <<< \"Distance\"\n",
        "\n",
        "name=\"${output_prefix}_distance_${pair_label}_0\"\n",
        "{\n",
        "  echo \"Time_(ns),RMSD_(A)\"\n",
        "  sed '/^[@#]/d' \"dist.xvg\" | awk '{ print $1/1000\",\"$2*10 }'\n",
        "} > \"dat.${name}.csv\"\n",
        "\n",
        "for i in $(seq 5); do echo \"\"; done\n",
        "{\n",
        "  echo \"***\"\n",
        "  echo \"Information: compared these atoms:\"\n",
        "  awk \"NR == ${atom_a} + 2\" \"conf.gro\"\n",
        "  awk \"NR == ${atom_b} + 2\" \"conf.gro\"\n",
        "  echo \"Made:  dat.${name}.csv\"\n",
        "  echo \"***\"\n",
        "} >&3"
      ],
      "metadata": {
        "id": "clJR7IY6oo8R",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "_Iyx5YRwasqH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@markdown **Plot computed atom pair data.** Tabular output from an atom pair distance analysis (above) may be visualised as a plot.\n",
        "#@markdown\n",
        "#@markdown Multiple atom pair labels may be specified as a comma separated list, in which case multiple plots will be made on the same axes.\n",
        "pair_labels = \"A-B\" #@param {type: \"string\"}\n",
        "#@markdown Optionally, plot a threshold line (on the $y$-axis) at a particular distance value. (0 to disable.)\n",
        "threshold_A = 0. #@param {type: \"number\"}\n",
        "#@markdown Optionally, cut off the plot at a particular maximum time. (0 to disable.)\n",
        "max_t_ns = 0 #@param {type: \"number\"}\n",
        "# defaults:    \"A-B\"    0.    0\n",
        "# \"\"\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* `dat.{output_prefix}_distance_{pair_labels}_0.csv` \\\\\n",
        "#@markdown *Output:* `dat.{output_prefix}_rmsd_{pair_labels}_0.png`</small>\n",
        "\n",
        "pair_labels = pair_labels.replace(\" \", \"\").split(\",\")\n",
        "\n",
        "\n",
        "#\n",
        "# Get the data from the various distance files\n",
        "#\n",
        "\n",
        "import pandas as pd\n",
        "\n",
        "t_ = \"Time_(ns)\"\n",
        "y_ = \"RMSD_(A)\"\n",
        "\n",
        "labels, tables = list(), list()\n",
        "for p in pair_labels:\n",
        "  try:\n",
        "    df = pd.read_csv(f\"dat.{output_prefix}_distance_{p}_0.csv\")\n",
        "  except:\n",
        "    continue\n",
        "\n",
        "  labels.append(p)\n",
        "\n",
        "  # Apply $t$ range\n",
        "  if max_t_ns > 0.:\n",
        "    df = df[(df[t_] <= max_t_ns)]\n",
        "  tables.append(df)\n",
        "\n",
        "\n",
        "#\n",
        "# Plot\n",
        "#\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "from matplotlib.ticker import AutoLocator\n",
        "\n",
        "def moving_average(df):\n",
        "  return df.rolling(21, center=True).mean()\n",
        "\n",
        "plt.figure()\n",
        "\n",
        "with sns.color_palette(\"deep\", n_colors=len(tables)):\n",
        "  for l, df in zip(labels, tables):\n",
        "    plt.plot(df[t_], df[y_], alpha=0.1)\n",
        "\n",
        "  for l, df in zip(labels, tables):\n",
        "    smooth = moving_average(df)\n",
        "    plt.plot(smooth[t_], smooth[y_], label=l)\n",
        "\n",
        "  plt.legend()\n",
        "\n",
        "  if threshold_A != 0.:\n",
        "    plt.axhline(y=threshold_A, color=\"black\", linestyle=\":\")\n",
        "\n",
        "# Title and axis labels\n",
        "plt.title(os.path.basename(project_folder).replace(\"_\", \" \").replace(\"vs \", \"vs. \"))\n",
        "plt.xlabel(\"Time (ns)\")\n",
        "plt.ylabel(f\"Distance of pair of atoms (A)\")\n",
        "\n",
        "# Set x-axis ticks\n",
        "ticker = AutoLocator()\n",
        "ticker.set_params(steps=[1, 2, 3, 5, 10] if round(df[t_].values[-1]) % 3 == 0 else [1, 2, 5, 10], min_n_ticks=6)\n",
        "ax = plt.gca()\n",
        "ax.xaxis.set_major_locator(ticker)\n",
        "\n",
        "\n",
        "#\n",
        "# Output\n",
        "#\n",
        "\n",
        "pairs_name = \"_\".join(pair_labels)\n",
        "name = f\"{output_prefix}_distance_{pairs_name}_0\"\n",
        "plt.gcf().savefig(f\"dat.{name}.png\")\n",
        "plt.show()\n",
        "\n",
        "print(\"***\")\n",
        "print(f\"Made:  dat.{name}.png\")\n",
        "print(\"***\")"
      ],
      "metadata": {
        "cellView": "form",
        "id": "MU2Qrbtwj7aB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "V-nsZ2JCatqY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$output_prefix\"\n",
        "output_prefix=\"$1\"\n",
        "exec 3>&1\n",
        "if $VERBOSE; then\n",
        "  exec 2>&1  # fix stderr and stdout appearing out of order\n",
        "else\n",
        "  exec >/dev/null 2>/dev/null  # silence\n",
        "fi\n",
        "set +H  # turn off history expansion so phrases that look like \"!...\" don't break\n",
        "\n",
        "#@markdown **Compute an energy analysis** for a **protein-ligand complex**, by rewatching an existing simulation. The calculated energies are the averages for the specified time range.\n",
        "#@markdown\n",
        "#@markdown (This may take a minute.)\n",
        "energy_group=\"LIG\" #@param {type: \"string\"}\n",
        "start_ps=0 #@param {type: \"integer\"}\n",
        "end_ps=1000 #@param {type: \"integer\"}\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* `index.ndx` (is created if missing)\n",
        "\n",
        "source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "\n",
        "\n",
        "#\n",
        "# Energy monitoring rerun of simulation with protein\n",
        "#\n",
        "\n",
        "if [[ ! -s \"index.ndx\" ]]; then\n",
        "  gmx make_ndx -f \"conf.gro\" -o \"index.ndx\" <<< \"q\"\n",
        "fi\n",
        "\n",
        "if ! fgrep -q \"Surround\" \"index.ndx\"; then\n",
        "  group_id=$(fgrep '[' \"index.ndx\" | wc -l)\n",
        "  { echo \"! \\\"${energy_group}\\\"\"; echo \"name ${group_id} Surround\"; echo \"q\"; } | gmx make_ndx -f \"conf.gro\" -n \"index.ndx\" -o \"index.ndx\"\n",
        "fi\n",
        "\n",
        "( cat \"grompp.mdp\"; echo \"energygrps = Surround ${energy_group}\" ) > \"lie_protein.mdp\"\n",
        "[[ ! -s \"restraint.gro\" ]] && cp \"conf.gro\" \"restraint.gro\"\n",
        "gmx grompp -f \"lie_protein.mdp\" -o \"lie_protein.tpr\" -c \"conf.gro\" -p \"topol.top\" -r \"restraint.gro\" -n \"index.ndx\" -maxwarn 999\n",
        "\n",
        "rerun=\"${output_prefix}_reference.xtc\"\n",
        "gmx trjconv -s \"${output_prefix}.tpr\" -f \"${rerun}\" -e $end_ps -o \"rerun.trr\" <<< \"0\"\n",
        "\n",
        "gmx mdrun -deffnm \"lie_protein\" -rerun \"rerun.trr\"  ||  exit $?\n",
        "\n",
        "\n",
        "#\n",
        "# Output the energy\n",
        "#\n",
        "\n",
        "read -r protein_qq protein_lj < <({\n",
        "  echo \"Coul-SR:Surround-${energy_group}\";\n",
        "  echo \"LJ-SR:Surround-${energy_group}\";\n",
        "  echo;\n",
        "} | gmx energy -f \"lie_protein.edr\" -b $start_ps -e $end_ps 2>/dev/null | tail -n2 | awk '{ print $2 }' | tr \"\\n\" \" \")\n",
        "\n",
        "protein_total=\"$(perl -le \"print(${protein_qq} + ${protein_lj})\")\"\n",
        "\n",
        "for i in $(seq 5); do echo \"\"; done\n",
        "{\n",
        "  echo \"***\"\n",
        "  echo \"Total interaction energy (ligand with protein):  ${protein_total}\"\n",
        "  echo \"Electrostatics term:  ${protein_qq}\"\n",
        "  echo \"Van der Waals term:  ${protein_lj}\"\n",
        "  echo \"***\"\n",
        "} >&3"
      ],
      "metadata": {
        "id": "QaOVYBg63RWU",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "  \\"
      ],
      "metadata": {
        "id": "Chj-Zo_bphsn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash -s \"$output_prefix\"\n",
        "output_prefix=\"$1\"\n",
        "exec 3>&1\n",
        "if $VERBOSE; then\n",
        "  exec 2>&1  # fix stderr and stdout appearing out of order\n",
        "else\n",
        "  exec >/dev/null 2>/dev/null  # silence\n",
        "fi\n",
        "set +H  # turn off history expansion so phrases that look like \"!...\" don't break\n",
        "\n",
        "#@markdown **Compute an energy analysis** for the **ligand in only solvent**, by computing a new simulation in a box of water.\n",
        "#@markdown\n",
        "#@markdown (This may take ten minutes.)\n",
        "energy_group=\"LIG\" #@param {type: \"string\"}\n",
        "#@markdown <font size=5> </font>\n",
        "#@markdown <small>*Input:* optionally, `index.ndx`\n",
        "\n",
        "source \"/usr/local/gromacs/bin/GMXRC.bash\"\n",
        "[[ -s \"index.ndx\" ]] && ndxarg=\"-n index.ndx\" || ndxarg=\"\"\n",
        "\n",
        "\n",
        "#\n",
        "# Simulate ligand in solvent\n",
        "#\n",
        "\n",
        "gmx trjconv -s \"${output_prefix}.tpr\" -f \"${output_prefix}_reference.xtc\" $ndxarg -dump 1000 -o \"ligand.gro\" <<< \"${energy_group}\"\n",
        "\n",
        "gmx editconf -f \"ligand.gro\" -o \"out.gro\" -c -d 1.0 -bt \"dodecahedron\"\n",
        "mv \"out.gro\" \"ligand.gro\"\n",
        "\n",
        "sed -E '/\\[\\s*molecules\\s*\\]/q' \"topol.top\" > \"ligand.top\"\n",
        "echo \"${energy_group} 1\" >> \"ligand.top\"\n",
        "\n",
        "# Put the full CHARMM36 forcefield first, deleting other forcefield.itps, to ensure TIP3 and ions are included\n",
        "# (This makes the assumption that the input system uses CHARMM36)\n",
        "awk '! x && $0 ~ /^#include/ { x=1; print \"#include \\\"'\"${START}\"'/charmm36.ff/forcefield.itp\\\"\" }; $0 !~ /forcefield.itp/ { print }' \"ligand.top\" > \"ligand.top.tmp\"\n",
        "mv \"ligand.top.tmp\" \"ligand.top\"\n",
        "\n",
        "# Solvate with TIP3\n",
        "sed 's/SOL /TIP3/g' \"/usr/local/gromacs/share/gromacs/top/spc216.gro\" > \"TIP3.gro\"\n",
        "gmx solvate -cp \"ligand.gro\" -cs \"TIP3.gro\" -o \"out.gro\" -p \"ligand.top\"\n",
        "mv \"out.gro\" \"ligand.gro\"\n",
        "\n",
        "\n",
        "#\n",
        "# Minimisation\n",
        "#\n",
        "\n",
        "cat > \"ligand_min.mdp\" <<EOS\n",
        "; See GROMACS_for_CHARMM-GUI.ipynb\n",
        "\n",
        "integrator            =  steep\n",
        "emtol                 =  1000.0\n",
        "nsteps                =  5000\n",
        "\n",
        "cutoff-scheme         =  Verlet\n",
        "nstlist               =  10\n",
        "coulombtype           =  PME\n",
        "rcoulomb              =  1.2\n",
        "vdwtype               =  Cut-off\n",
        "vdw-modifier          =  Force-switch\n",
        "rvdw-switch           =  1.0\n",
        "rvdw                  =  1.2\n",
        "\n",
        "constraints           =  h-bonds\n",
        "constraint-algorithm  =  LINCS\n",
        "EOS\n",
        "\n",
        "gmx grompp -f \"ligand_min.mdp\" -o \"ligand_min.tpr\" -c \"ligand.gro\" -p \"ligand.top\" -maxwarn 999\n",
        "\n",
        "gmx genion -s \"ligand_min.tpr\" -o \"ligand_with_ions.gro\" -p \"ligand.top\" -pname \"POT\" -nname \"CLA\" -neutral <<< \"TIP3\"\n",
        "mv \"ligand_min.gro\" \"ligand_min.gro.bak\"\n",
        "mv \"ligand_with_ions.gro\" \"ligand_min.gro\"\n",
        "\n",
        "gmx grompp -f \"ligand_min.mdp\" -o \"ligand_min.tpr\" -c \"ligand_min.gro\" -p \"ligand.top\" -maxwarn 999\n",
        "\n",
        "gmx mdrun -v -stepout 1000 -deffnm \"ligand_min\"  ||  exit $?\n",
        "\n",
        "\n",
        "#\n",
        "# Equilibration\n",
        "#\n",
        "\n",
        "cat > \"ligand_eq.mdp\" <<EOS\n",
        "; See GROMACS_for_CHARMM-GUI.ipynb\n",
        "\n",
        "integrator  =  md\n",
        "dt          =  0.001     ; 1 fs\n",
        "nsteps      =  100000    ; 100 ps simulation time\n",
        "comm-mode   =  Linear\n",
        "\n",
        "cutoff-scheme  =  Verlet\n",
        "nstlist        =  20\n",
        "coulombtype    =  PME\n",
        "rcoulomb       =  1.2\n",
        "vdwtype        =  Cut-off\n",
        "vdw-modifier   =  Force-switch\n",
        "rvdw-switch    =  1.0\n",
        "rvdw           =  1.2\n",
        "\n",
        "tcoupl   = Berendsen\n",
        "tc-grps  = System\n",
        "tau-t    = 1.0\n",
        "ref-t    = 303.15\n",
        "\n",
        "pcoupl            = Berendsen\n",
        "pcoupltype        = isotropic\n",
        "tau-p             = 5.0\n",
        "compressibility   = 4.5e-5\n",
        "ref-p             = 1.0\n",
        "refcoord-scaling  = com\n",
        "\n",
        "constraints           =  h-bonds\n",
        "constraint-algorithm  =  LINCS\n",
        "continuation          =  yes\n",
        "EOS\n",
        "\n",
        "gmx grompp -f \"ligand_eq.mdp\" -o \"ligand_eq.tpr\" -c \"ligand_min.gro\" -p \"ligand.top\" -maxwarn 999\n",
        "\n",
        "gmx mdrun -v -stepout 1000 -deffnm \"ligand_eq\"  ||  exit $?\n",
        "\n",
        "\n",
        "#\n",
        "# Production\n",
        "#\n",
        "\n",
        "cat > \"ligand_prod.mdp\" <<EOS\n",
        "; See GROMACS_for_CHARMM-GUI.ipynb\n",
        "\n",
        "integrator  =  md\n",
        "dt          =  0.002     ; 2 fs\n",
        "nsteps      =  50000     ; 100 ps simulation time\n",
        "comm-mode   =  Linear\n",
        "\n",
        "nstxout  =  5000    ; 10 frames total\n",
        "nstvout  =  5000\n",
        "\n",
        "cutoff-scheme  =  Verlet\n",
        "nstlist        =  20\n",
        "coulombtype    =  PME\n",
        "rcoulomb       =  1.2\n",
        "vdwtype        =  Cut-off\n",
        "vdw-modifier   =  Force-switch\n",
        "rvdw-switch    =  1.0\n",
        "rvdw           =  1.2\n",
        "\n",
        "tcoupl   =  V-Rescale\n",
        "tc-grps  =  System\n",
        "tau-t    =  1.0\n",
        "ref-t    =  303.15\n",
        "\n",
        "pcoupl           =  C-Rescale\n",
        "pcoupltype       =  isotropic\n",
        "tau-p            =  5.0\n",
        "compressibility  =  4.5e-5\n",
        "ref-p            =  1.0\n",
        "\n",
        "constraints           =  h-bonds\n",
        "constraint-algorithm  =  LINCS\n",
        "continuation          =  yes\n",
        "EOS\n",
        "\n",
        "gmx grompp -f \"ligand_prod.mdp\" -o \"ligand_prod.tpr\" -c \"ligand_eq.gro\" -p \"ligand.top\" -maxwarn 999\n",
        "\n",
        "gmx mdrun -v -stepout 1000 -deffnm \"ligand_prod\"  ||  exit $?\n",
        "\n",
        "\n",
        "#\n",
        "# Energy monitoring rerun of simulation with ligand\n",
        "#\n",
        "\n",
        "# TODO FIXME\n",
        "# This is a quick implementation that doesn't work properly if the energy_group given is not a default group.\n",
        "# It's probably better to do it the other way around -- Surround = Water + Ions, ${energy_group} = ~Surround\n",
        "\n",
        "echo \"q\" | gmx make_ndx -f \"ligand_min.gro\" -o \"ligand.ndx\"\n",
        "group_id=$(fgrep '[' \"ligand.ndx\" | wc -l)\n",
        "{ echo \"! \\\"${energy_group}\\\"\"; echo \"name ${group_id} Surround\"; echo \"q\"; } | gmx make_ndx -f \"ligand_min.gro\" -n \"ligand.ndx\" -o \"ligand.ndx\"\n",
        "\n",
        "( cat \"ligand_prod.mdp\"; echo \"energygrps = ${energy_group} Surround\" ) > \"lie_ligand.mdp\"\n",
        "gmx grompp -f \"lie_ligand.mdp\" -o \"lie_ligand.tpr\" -c \"ligand_min.gro\" -p \"ligand.top\" -n \"ligand.ndx\" -maxwarn 999\n",
        "\n",
        "gmx mdrun -deffnm \"lie_ligand\" -rerun \"ligand_prod.trr\"  ||  exit $?\n",
        "\n",
        "\n",
        "#\n",
        "# Output the energy\n",
        "#\n",
        "\n",
        "read -r ligand_qq ligand_lj < <({\n",
        "  echo \"Coul-SR:${energy_group}-Surround\";\n",
        "  echo \"LJ-SR:${energy_group}-Surround\";\n",
        "  echo;\n",
        "} | gmx energy -f \"lie_ligand.edr\" 2>/dev/null | tail -n2 | awk '{ print $2 }' | tr \"\\n\" \" \")\n",
        "\n",
        "ligand_total=\"$(perl -le \"print(${ligand_qq} + ${ligand_lj})\")\"\n",
        "\n",
        "for i in $(seq 5); do echo \"\"; done\n",
        "{\n",
        "  echo \"***\"\n",
        "  echo \"Total interaction energy (ligand with solvent):  ${ligand_total}\"\n",
        "  echo \"Electrostatics term:  ${ligand_qq}\"\n",
        "  echo \"Van der Waals term:  ${ligand_lj}\"\n",
        "  echo \"***\"\n",
        "} >&3"
      ],
      "metadata": {
        "cellView": "form",
        "id": "s5mVqG53pewS"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}