{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a Simulation in MDMC"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create an example universe, filled with water with an SPCE forcefield. If you'd like to learn more about any of this, please read the following:\n",
    "\n",
    "- [Creating Atomic Configurations](creating-atomic-configurations.ipynb)\n",
    "- [Defining Molecule Interactions](defining-molecule-interactions.ipynb)\n",
    "- [Applying a Forcefield](applying-a-forcefield.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import *\n",
    "from MDMC.MD.force_fields.three_site_water import ThreeSiteWater, add_three_site_water_ff\n",
    "\n",
    "universe = Universe(dimensions=24.83602653)\n",
    "universe.fill(ThreeSiteWater(), num_density=0.03356718472021752)\n",
    "add_three_site_water_ff(universe, cutoff=10.0, ewald=1e-4)"
   ]
  },
  {
   "attachments": {},
   "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, as above.  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='openmm', time_step=1., temperature=300., traj_step=10)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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(f'Velocity of first atom in MDMC universe is {universe.atoms[0].velocity}')\n",
    "state = simulation.engine.openmm_simulation.context.getState(velocities=True)\n",
    "print(f'Velocity of first atom in MD engine is {state.getVelocities()[0]}')\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(f'Changing the velocity of the first atom to {velocity}')\n",
    "universe.atoms[0].velocity = velocity\n",
    "simulation_2 = Simulation(universe, engine='openmm', time_step=1., temperature=300.,\n",
    "                          pressure=101325., traj_step=10, thermostat='nose',\n",
    "                          barostat='nose', t_damp=100, p_damp=1000, openmm_platform=\"CPU\")\n",
    "print(f'Velocity of first atom in MDMC universe is {universe.atoms[0].velocity}')\n",
    "state = simulation_2.engine.openmm_simulation.context.getState(velocities=True)\n",
    "print(f'Velocity of first atom in MD engine is {state.getVelocities()[0]}')\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)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the equilibration phase, statistics and trajectories are not captured. If you have not explicitly specified a thermostat or barostat in your calculation, a Berendsen thermostat will be used for the duration of the equilibration. This Berendsen thermostat will not be used during the production run."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation can be run like so. (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 (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.11.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}