validation/physical_validation_of_simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking physical validity of simulations\n",
    "Michael Shirts, CU Boulder\n",
    "OSU MC/MD Summer Workshop, July 2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enable plotting in notebook\n",
    "%matplotlib notebook\n",
    "\n",
    "# change shell to the correct version\n",
    "%set_env SHELL=/bin/bash"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking whether a system of water molecules is physically consistent using GROMACS input files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing checks and loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by importing the `physical_validation` package.\n",
    "Please refer to http://physical-validation.readthedocs.io/ for the full documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import physical_validation as pv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a GROMACS parser object, which needs the location of the GROMACS executable and the location of the topology include folder as inputs. Here, we assume that `gmx` is in the PATH, and that the topology folder is in its standard location. Change this if your local installation differs from this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmx = pv.data.GromacsParser(exe='gmx_mpi',includepath='/opt/conda/share/gromacs/top/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll test simulations ran with two thermostats:\n",
    "`vr` stands for velocity-rescale, `be` for Berendsen thermostat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "algos = ['vr', 'be']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll test simulations performed in two ensembles: NVT and NPT. \n",
    "Note that in NPT, the `vr` thermostat was complemented with a Parinello-Rahman barostat, while the `be` thermostat was complemented with a Berendsen barostat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ensembles = ['NVT', 'NPT']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of simulations performed in each ensemble:\n",
    "\n",
    "* We need 2 simulations for NVT: Simulations at $T$ and $T+\\Delta T$\n",
    "\n",
    "* We use 4 simulations for NPT: $(T, P)$, $(T+\\Delta T, P)$, $(T, P+\\Delta P)$, $(T+\\Delta T, P+\\Delta P)$\n",
    "\n",
    "There are actually three tests we can perform with NPT simulations.  \n",
    "\n",
    "* Does the _volume_ (V) distribution change correctly when we vary the pressure $P$?\n",
    "* Does the _enthalpy_ $(U+PV)$ distribution change correctly when we vary the $T$? \n",
    "* Does the 2D distribution of $U$ and $V$ change correctly when we change both $T$ and $P$? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = {'NVT': 2, 'NPT': 4}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The directories the data is stored are of the form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = \"md_vr_NVT_1\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare a list we will store the parsed data objects in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_data = list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run checks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have prepared the setup variables, we will perform the checks. We will start with simulation data generated with the velocity-rescale thermostat of Bussi et al., and we will perform various checks of constant-volume ($NVT$) simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the simulation results using the GROMACS parser. This uses the `mdp` parameter file and the `top` topology file to gather information about the system and the simulation settings, and read the results from the `edr` file (trajectory of energy / volume / pressure / . . . ) and the `gro` file (position snapshot - used to read the box volume in $NVT$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(1,3):\n",
    "    d = \"systems/water/md_NVT_vr_\"+str(n)\n",
    "    simulation_data.append(gmx.get_simulation_data(\n",
    "        mdp=d + \"/mdout.mdp\",\n",
    "        edr=d + \"/system.edr\",\n",
    "        gro=d + \"/system.gro\",\n",
    "        top=\"systems/water/top/system.top\",\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "## Test the kinetic energy distribution of the simulation results.\n",
    "\n",
    "The first input is the simulation results we just read in, `strict` determines whether we test the full distribution (True) or only determine the mean and the variance of the distribution (False), `verbosity` sets the level of detail of the output  (with verbosity=0 being quiet and verbosity=3 being the most chatty), and the filename is being used to plot the resulting distribution for visual inspection (which we are also sending to the notebook).\n",
    "\n",
    "Notice that we are doing a little bit of equilibration detection and subsampling here, as discussed earlier!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('==> Kinetic energy test of NVT simulations with velocity rescale simulation ')\n",
    "print('==> Simulation 1 ')\n",
    "pv.kinetic_energy.distribution(simulation_data[0], strict=True, verbosity=2)\n",
    "print('==> Simulation 2 ')\n",
    "pv.kinetic_energy.distribution(simulation_data[1], strict=True, verbosity=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also do the non-strict test, as that can be a little more physically informative; we can see if the mean and standard deviation of the kinetic energy distributions are consistent with the temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.kinetic_energy.distribution(simulation_data[0], strict=False)\n",
    "pv.kinetic_energy.distribution(simulation_data[1], strict=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the simulations indeed have the correct temperatures, and the $T_{\\mu}$ and $T_{\\sigma}$ are consistent with each other. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to be sure, let's look at things visually.  We can send the data to the screen, and we can specify a filename to save the picture to, as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pv.kinetic_energy.distribution(simulation_data[0], strict=False, verbosity=2,filename=\"vr_NVT_ke_0\",screen=True)\n",
    "pv.kinetic_energy.distribution(simulation_data[1], strict=False, verbosity=2,filename=\"vr_NVT_ke_1\",screen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hands on exercise\n",
    "\n",
    "Check the kinetic energy distributions of the Berendsen weak-coupling thermostat (stored in directories with `be` instead of `vr` and see what you find!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next: Ensemble validation of NPT simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the two simulations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.ensemble.check(simulation_data[0], simulation_data[1], verbosity=2, filename='vr_NVT_ensemble', screen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: \n",
    "    \n",
    " Does the Berendsen thermostat obey the correct distribution of energies? Look for data in the `br` directories. \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NPT validation\n",
    "\n",
    "There are three types of NPT validation we could do: Is the enthalpy distribution correct, is the volume distribution correct, and is the joint distribution of $(U,V)$ correct?.  To do that, we need four simulations. \n",
    " - directories ending in `_1`: $T=300K$, $P=1$ atm\n",
    " - directories ending in `_2`: $T=308K$, $P=300$ atm\n",
    " - directories ending in `_3`: $T=300K$, $P=300$ atm\n",
    " - directories ending in `_4`: $T=308K$, $P=308$ atm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rezero the simulation data files\n",
    "simulation_data=list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(1,5):\n",
    "    d = \"systems/water/md_NPT_vr_\"+str(n)\n",
    "    simulation_data.append(gmx.get_simulation_data(\n",
    "        mdp=d + \"/mdout.mdp\",\n",
    "        edr=d + \"/system.edr\",\n",
    "        gro=d + \"/system.gro\",\n",
    "        top=\"systems/water/top/system.top\",\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the code automatically checks what ensemble the data was simulated in.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.ensemble.check(simulation_data[0], simulation_data[1], verbosity=2, filename='vr_NPT_ensemble', screen=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.ensemble.check(simulation_data[0], simulation_data[2], verbosity=2, filename='vr_NPT_ensemble', screen=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.ensemble.check(simulation_data[0], simulation_data[3], verbosity=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doing the same checks with just flat files of numbers.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was done with GROMACS simulations. But maybe we don't want to do GROMACS simulations!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parser = pv.data.FlatfileParser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need to manually enter information that was read automatically off of the GROMACS files. \n",
    "\n",
    "Our test system consists of 900 H$_2$O molecules whose bonds are fully constrained. Also, during the simulation, we kept the translation of the center of mass to zero, so we need to reduce the number of degrees of freedom by 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = pv.data.SystemData(\n",
    "    natoms=900*3,\n",
    "    nconstraints=900*3,\n",
    "    ndof_reduction_tra=3,\n",
    "    ndof_reduction_rot=0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The physical validation tests need some information on the units that were used in the simulation. While the strings are only used for output, the conversion respective to GROMACS units are relevant for the calculations. Please see the documentation for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "units = pv.data.UnitData(\n",
    "    kb=8.314462435405199e-3,\n",
    "    energy_str='kJ/mol',\n",
    "    energy_conversion=1.0,\n",
    "    length_str='nm',\n",
    "    length_conversion=1.0,\n",
    "    volume_str='nm^3',\n",
    "    volume_conversion=1.0,\n",
    "    temperature_str='K',\n",
    "    temperature_conversion=1.0,\n",
    "    pressure_str='bar',\n",
    "    pressure_conversion=1.0,\n",
    "    time_str='ps',\n",
    "    time_conversion=1.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the flat file example, we will only look at the simulations performed under NVT conditions using the velocity-rescale thermostat. There are two simulations: One ran at 300K, and one ran at 308K."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ensemble_1 = pv.data.EnsembleData(\n",
    "    ensemble='NVT',\n",
    "    natoms=900*3,\n",
    "    volume=3.01125**3,\n",
    "    temperature=300\n",
    ")\n",
    "ensemble_2 = pv.data.EnsembleData(\n",
    "    ensemble='NVT',\n",
    "    natoms=900*3,\n",
    "    volume=3.01125**3,\n",
    "    temperature=308\n",
    ")\n",
    "dir_1 = 'systems/water/md_NVT_vr_1'\n",
    "dir_2 = 'systems/water/md_NVT_vr_2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now read in our flat data files (1-dimensional ascii files containing energy trajectories), and create a simulation result representation usable by the physical validation tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_1 = parser.get_simulation_data(\n",
    "    units=units, ensemble=ensemble_1, system=system,\n",
    "    kinetic_ene_file=dir_1 + '/kinetic_energy.dat',\n",
    "    potential_ene_file=dir_1 + '/potential_energy.dat',\n",
    "    total_ene_file=dir_1 + '/total_energy.dat'\n",
    ")\n",
    "result_2 = parser.get_simulation_data(\n",
    "    units=units, ensemble=ensemble_2, system=system,\n",
    "    kinetic_ene_file=dir_2 + '/kinetic_energy.dat',\n",
    "    potential_ene_file=dir_2 + '/potential_energy.dat',\n",
    "    total_ene_file=dir_2 + '/total_energy.dat'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now run the tests!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('==> Kinetic energy test of simulation ' + dir_1)\n",
    "pv.kinetic_energy.distribution(result_1, strict=False, verbosity=2,\n",
    "                               filename='ke_flat_vr_NVT_1', screen=True)\n",
    "print('==> Kinetic energy test of simulation ' + dir_2)\n",
    "pv.kinetic_energy.distribution(result_2, strict=False, verbosity=2,\n",
    "                               filename='ke_flat_vr_NVT_2', screen=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('==> Potential energy test')\n",
    "pv.ensemble.check(result_1, result_2,\n",
    "                  verbosity=2, filename='pe_flat_vr_NVT', screen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking equipartition of velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load in the data.  Note that we are now loading a `.trr` file, which contains the velocities.  You would need to make sure you tell gromacs to output velocities!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a set of directories with four sets of simulation.\n",
    "* `vr` stands for the data generated with the velocity-rescale thermostat.\n",
    "* `be` stands for the data generated using the Berendsen weak-coupling thermostat.\n",
    "* `_1` denotes a simulation using a single thermostat for the entire system.\n",
    "* `_2` denotes a simulation using two separate thermostats, one for the protein and one for the solvent.\n",
    "\n",
    "Note that the topology and the trajectory have been modified to only contain the protein, as we are only interested in the equiparition of the solute, not the solvent. This reduces file size and execution time considerably."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_protein = gmx.get_simulation_data(mdp='systems/trp-cage/vr_1/protein.mdp',\n",
    "                                              top='systems/trp-cage/vr_1/trp-cage.top',\n",
    "                                              edr='systems/trp-cage/vr_1/run.edr',\n",
    "                                              trr='systems/trp-cage/vr_1/protein.trr')       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One minor change we need to do because we removed the water. In the simulations, the center of mass was artificially kept immobile to avoid the build up of numerical errors, effectively reducing the number of translational degrees of freedom of the system by 3. As we are only looking at the protein here, we're using a mass-dependent fraction of these three degrees of freedom (weight of protein: 2170.4375 amu; total weight of the system: 96390.4017 amu).\n",
    "\n",
    "If we were doing this with the full system (solvent + water), then we could just remove 3 degrees of freedom. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_protein.system.ndof_reduction_tra *= 2170.4375 / 96390.4017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.kinetic_energy.equipartition(results_protein,strict=False, filename='equipartition_vr_1', screen=True,\n",
    "                                    verbosity=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at Berendsen weak coupling!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_protein = gmx.get_simulation_data(mdp='systems/trp-cage/be_1/protein.mdp',\n",
    "                                              top='systems/trp-cage/be_1/trp-cage.top',\n",
    "                                              edr='systems/trp-cage/be_1/run.edr',\n",
    "                                              trr='systems/trp-cage/be_1/protein.trr')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_protein.system.ndof_reduction_tra *= 2170.4375 / 96390.4017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv.kinetic_energy.equipartition(results_protein,strict=False, filename='equipartition_be_1', screen=True,\n",
    "                                    verbosity=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: \n",
    "Now see what happens when we use two different thermostatting groups, one each for the solvent and the protein, with both the velocity rescale thermostat (`vr`) and the Berendsen weak-coupling thermostat (`be`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}