{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QMMM workflow using GROMACS and VOTCA-XTP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is this tutorial about\n",
    "In this tutorial, we will learn how to set and perform excited state calculation using the Votca XTP library. We will use methane as our QM region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "* You will need to install **VOTCA** using the instructions described [here](https://github.com/votca/votca/blob/master/share/doc/INSTALL.rst)\n",
    "* Once the installation is completed you need to activate the VOTCA enviroment by running the `VOTCARC.bash` script that has been installed at the bin subfolder for the path that you have provided for the installation step above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interacting with the XTP command line interface\n",
    "The XTP package offers the following command line interface that the user can interact with:\n",
    "* [xtp_map](https://www.votca.org/xtp/xtp_map.html)\n",
    "* [xtp_parallel](https://www.votca.org/xtp/xtp_parallel.html)\n",
    "* [xtp_run](https://www.votca.org/xtp/xtp_run.html)\n",
    "* [xtp_tools](https://www.votca.org/xtp/xtp_tools.html)\n",
    "\n",
    "Run the following command to view the help message of `xtp_tools`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_tools -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note\n",
    "> * In Jupyter the `!` symbol means: *run the following command as a standard unix command*\n",
    "> * In Jupyter the command `%env` set an environmental variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the environment\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove previous hdf5 file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the topology from the Gromacs file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "runs the mapping from MD coordinates to segments and creates an [hdf5 file](https://www.hdfgroup.org/solutions/hdf5/). You can explore the generated `state.hdf5` file with e.g. hdf5itebrowser. In Python, you can use the [h5py library](https://www.h5py.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_map -t MD_FILES/topol.tpr -c MD_FILES/conf.gro -s system.xml -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Check the mapping\n",
    "\n",
    "Let us first output `.pdb` files for the segments, qmmolecules and classical segments in order to check the mapping. Use `xtp_run -d mapchecker` to see all options `mapchecker` calculator takes. We use the `-c` option to change one option on the commandline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [mapchecker section of the manual](https://votca.github.io/xtp/mapchecker.html) you can find a table with the `mapchecker` input variables and their corresponding defaults. Finally, the following command run the check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_run -e mapchecker -c map_file=system.xml -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neighborlist Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following step is to determine the neighbouring pairs for exciton transport. See the [neighborlist options](https://www.votca.org/xtp/neighborlist.html) for further information. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can run the calculation using 4 threads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_run -e neighborlist -c exciton_cutoff=0.5 constant=0.6 -f state.hdf5 -t 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read reorganization energies\n",
    "In this step we will read the in site reorganization energies and store them in the `state.hdf5` file. We just need to copy the input file and execute the calculation. The side energies have to be calculated by the user beforehand and put into an xml file. We added them to `system.xml`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_run -e einternal -c energies_file=system.xml -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute site energy\n",
    "In this step we will perform some *QMMM* calculations to compute the site energies. The `qmmm_mm.xml` file contains some predefined settings to perform the *MM* calculations. Let us first copy these settings into the state file. Instead of using the `-c` option we now use the `-o` option to read in options from an xml file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e qmmm -o qmmm_mm.xml -f state.hdf5 -j \"write\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous command generates a `qmmm_mm_jobs.xml` containing 4000 *MM* jobs to compute, if you examine that file, it should look something like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<jobs>\n",
    "  <job>\n",
    "    <id>0</id>\n",
    "    <tag>Methane_0:n</tag>\n",
    "    <input>\n",
    "      <site_energies>0:n</site_energies>\n",
    "      <regions>\n",
    "        <region>\n",
    "          <id>0</id>\n",
    "          <segments>0:n</segments>\n",
    "        </region>\n",
    "      </regions>\n",
    "    </input>\n",
    "    <status>AVAILABLE</status>\n",
    "  </job>    \n",
    "...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run just the first 4 jobs by settings all jobs `status` to `COMPLETE` except for the first four. This can be easily done with [sed](https://www.gnu.org/software/sed/manual/sed.html) as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -i \"s/AVAILABLE/COMPLETE/g\" qmmm_mm_jobs.xml\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_mm_jobs.xml\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_mm_jobs.xml\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_mm_jobs.xml\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_mm_jobs.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the jobs and save the results in the state file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e qmmm -o qmmm_mm.xml -f state.hdf5 -x 2 -j \"run\"\n",
    "!xtp_parallel -e qmmm -o qmmm_mm.xml -f state.hdf5 -j \"read\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Site energy and pair energy analysis\n",
    "In this step we generate an histogram and compute the correlation function of site energies and pair energy differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_run -e eanalyze -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now see a set of files prefixed with `eanalyze` containing the histrogram and correlation functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls eanalyze*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QM energy calculation\n",
    "Our next task is to perform the qm calculations for each segment that we have stored in the hdf5 file. The calculations take place in 3 stages: write the jobs to a file, perform the computation and finally save the results to the state file. We provided a small options file to make the computation cheaper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat eqm.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the GWBSE mode to `G0W0`,  the `ranges` to `full` and the `basisset` and `auxbasisset` to `3-21G` and `aux-def2-svp`. For more information, check the [eqm calculator options](https://votca.github.io/xtp/eqm.html). For the sake of computational time let just compute the `gw` approximation and the `singlet`. You can also request the `triplet` or `all`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will write the job in a file and enable only the first 2,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e eqm -o eqm.xml -f state.hdf5 -s 0 -j \"write\"\n",
    "!sed -i \"s/AVAILABLE/COMPLETE/g\" eqm.jobs\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' eqm.jobs\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' eqm.jobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let run these 2 jobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e eqm -o eqm.xml -f state.hdf5 -x 2 -s 0 -j run -q 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QM calculation for pairs\n",
    "In the following step we will run QM calculations for each pair in the hdf5 file. As the calculations on the previous step, we will first write the jobs in a file, then run them and finally store the results in the state file. First, we need to copy the input to our local folder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the previous section, we set the GWBSE mode to `G0W0`, the `ranges` to `full` and the `basisset` and `auxbasisset` to `3-21G` and `aux-def2-svp`. But we compute only the `gw` approximation, as the BSE is formed in the coupling step only once and we do not have to diagonalize it. For more information, check the [iqm calculator options](https://www.votca.org/xtp/iqm.html). We only compute the `singlet` couplings. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running the calculations, we need to specify in the `iqm` input which states to read into the jobfile for each segment type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat iqm.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's write the jobs to the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e iqm -o iqm.xml -f state.hdf5 -s 0 -j \"write\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the jobs that we just write down, let's make available only the first job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -i \"s/AVAILABLE/COMPLETE/g\" iqm.jobs\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' iqm.jobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run and store the jobs results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e iqm -o iqm.xml -f state.hdf5 -x 2 -s 0 -j run -q 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we read the results into the state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e iqm -o iqm.xml -f state.hdf5 -j \"read\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coupling\n",
    "We can now compute the classical coupling of transition in the aformentioned three stages,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to change in the `iexcitoncl` input the name `map_file` option and add the state. check all the available of the [iexcitoncl calculator](https://www.votca.org/xtp/iexcitoncl.html). We do this via the commandline using the `-c` option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e iexcitoncl -c map_file=system.xml states=Methane:n2s1 -f state.hdf5 -j \"write\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head -n 15 exciton.jobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run and save the jobs. For demo purposes we will run only the first job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -i \"s/AVAILABLE/COMPLETE/g\" exciton.jobs\n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' exciton.jobs\n",
    "!xtp_parallel -e iexcitoncl -c map_file=system.xml states=Methane:n2s1 -f state.hdf5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e iexcitoncl -c map_file=system.xml -f state.hdf5 -j \"read\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coupling analysis\n",
    "Using the coupling computed in the previous steps, we will generate an histogram for the squared couplings in logarithmic scale,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_run -e ianalyze -c states=e,h,s -f state.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QMMM calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally let us run a proper qmmm calculation using the qmmm calculator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat qmmm.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!xtp_parallel -e qmmm -o qmmm.xml -f state.hdf5 -j \"write\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets run just the first job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -i \"s/AVAILABLE/COMPLETE/g\" qmmm_jobs.xml                                                                                                                                                            \n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_jobs.xml   \n",
    "!sed -i '0,/COMPLETE/s/COMPLETE/AVAILABLE/' qmmm_jobs.xml  \n",
    "!xtp_parallel -e qmmm -o qmmm.xml -x 2 -f state.hdf5 -j run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, save the results. We could read them in but that is a bit pointless. Maybe check out how to turn a checkpoint file into an or orbfile (look at the scripts) and visualise it with the `gencube` tool. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!xtp_parallel -e qmmm -o OPTIONFILES/qmmm.xml -f state.hdf5 -j \"read\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
