{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a Simulation in MDMC\n",
    "\n",
    "This tutorial continues from the [Building a Universe](building-a-universe.ipynb) tutorial, using the example universe filled with SPCE water."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `%%capture` and `%run` commands below simply executes the [Building a Universe](building-a-universe.ipynb) notebook and captures the variables into this notebook.  They are only valid if they are executed in the same folder as the [Building a Universe](building-a-universe.ipynb) notebook.  Otherwise, please copy the last section of the [Building a Universe](building-a-universe.ipynb) to set the same state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# Run Building a universe notebook and hide output\n",
    "%run \"building-a-universe.ipynb\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Simulation object\n",
    "Simulations in MDMC are run using external MD engines (e.g. LAMMPS).  First create a universe (see [Building a Universe](building-a-universe.ipynb) - an example universe can be found at the bottom of the page).  This universe object must then be passed when creating a Simulation object, along with simulation properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Simulation class\n",
    "from MDMC.MD import Simulation\n",
    "\n",
    "# Create an NPT simulation\n",
    "simulation = Simulation(universe, engine='lammps', time_step=1., temperature=300.,\n",
    "                        pressure=101325., traj_step=10, thermostat='nose',\n",
    "                        barostat='nose', t_damp=100, p_damp=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NPT simulation is created with a Nose-Hoover thermostat and a Berendsen barostat, with a temperature of 300 K, a pressure of 101325 Pa, and a time_step of 1.0 fs.  The trajectory is recorded every 10 simulation steps. The number of steps over which the temperature and pressure are damped for the thermostat and barostats is 100 and 1000 respectively.\n",
    "\n",
    "NVE and NVT simulations can be created in a similar manner by omitting one or both of the thermostat and barostat.  If both are omitted the temperature must still be specified as this initialises the atomic velocities.\n",
    "\n",
    "MDMC allows detailed control of the atomic velocities when creating the `Universe`. In the case where no velocities were provided (i.e. all `Atom` objects have the default velocity of 0) then the starting velocities are determined by the MD engine (randomly chosen from a uniform distribution, and then scaled so that the velocities are consistent with the temperature provided to the `Simulation`). If some or all of the atoms have been set with MDMC then these velocities will be scaled to the correct temperature by the MD engine. In both cases, only the velocities of atoms within the MD engine are affected, and the state of the original `Universe` is unchanged:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the simulation object, which was created from an input Universe where all atom velocities were equal to zero\n",
    "print('Velocity of first atom in MDMC universe is {}'.format(universe.atoms[0].velocity))\n",
    "print('Velocity of first atom in MD engine is {}'.format(simulation.engine.lmp.atoms[0].velocity))\n",
    "\n",
    "# In comparision, create a new simulation object where (artifcially and for demonstration purposes) one atom has non-zero velocity\n",
    "velocity = (1, 0, -1)\n",
    "print('\\nChanging the velocity of the first atom to {}\\n'.format(velocity))\n",
    "universe.atoms[0].velocity = velocity\n",
    "simulation_2 = Simulation(universe, engine='lammps', time_step=1., temperature=300.,\n",
    "                          pressure=101325., traj_step=10, thermostat='nose',\n",
    "                          barostat='nose', t_damp=100, p_damp=1000)\n",
    "print('Velocity of first atom in MDMC universe is {}'.format(universe.atoms[0].velocity))\n",
    "print('Velocity of first atom in MD engine is {}'.format(simulation_2.engine.lmp.atoms[0].velocity))\n",
    "\n",
    "# Reset atom velocity back to zero\n",
    "universe.atoms[0].velocity = (0, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy minimisation and running a simulation\n",
    "The universe energy can be minimised by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Minimise the system during a 100-step MD run, minimising every 10 steps\n",
    "simulation.minimize(100, minimize_every=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation can be equilibrated by (this will take ~30s):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.run(1000, equilibration=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trajectories are not stored during equilibration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation can be run by (this will take ~60s):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the simulation for 2000 steps\n",
    "simulation.run(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trajectory\n",
    "An MDMC CompactTrajectory object can be created following a simulation run using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = simulation.trajectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The times of all of the steps of trajectory can be accessed with the trajectory.times attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory.times[1] - trajectory.times[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variation between MD engines\n",
    "In theory, all simulations on a MDMC Universe should be able to use any MD engine, although in practice this is limited by whether a particular MD engine supports a specific feature and if it has been implemented in the MD engine interface.  If a feature is not supported by or implemented for a specific MD engine, MDMC will raise a NotImplementedError."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.11.2 64-bit",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}