{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d6ace33-bd89-49ac-a148-06c99d7d071c",
   "metadata": {},
   "source": [
    "# Getting Started with CHEG231MD\n",
    "\n",
    "Eric M. Furst\\\n",
    "October 2024\n",
    "\n",
    "Test run of the simple MD simulation.\n",
    "\n",
    "The temperature will depend on density. At any density, change temp by \n",
    "scaling max velocity by sqrt(Tnew/Told) \n",
    "\n",
    "$$ \\langle v^2 \\rangle = v_\\mathrm{max}^2/2 $$\n",
    "\n",
    "Without a thermostat, both the temperature and pressure equilibrate.\n",
    "This makes calculating an isotherm difficult, but it shows the\n",
    "combined contributions of potential and kinetic energy\n",
    "\n",
    "rho and vmax values that give results close to T+ = 2:\\\n",
    "0.7 5.8\\\n",
    "0.6 5.47\\\n",
    "0.4 5.04"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7f0016e-5ba0-45b9-9d0e-df7bde2c901f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import CHEG231MD as md\n",
    "\n",
    "# create a simulation with # particles, density, and max speed\n",
    "# N^(1/3), rho+, max vel+\n",
    "\n",
    "#sim = md.MDSimulation(8,0.4,5.04)\n",
    "#sim = md.MDSimulation(8,0.6,5.47)\n",
    "sim = md.MDSimulation(8,0.7,5.8)\n",
    "\n",
    "# initialize variables and lists\n",
    "pres=[]\n",
    "temp=[]\n",
    "time=[]\n",
    "Pavg = 0\n",
    "Tavg = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a14e5dc1-7601-4972-911e-0f42eb32a6c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run time steps of the simulation\n",
    "for timestep in range(1,1001):\n",
    "\n",
    "    # advance simulation one timestep\n",
    "    sim.move()\n",
    "    \n",
    "    # only average after an equilibration time\n",
    "    if timestep > 100:\n",
    "        # production stage\n",
    "        Pavg += (sim.P-Pavg)/(timestep+1)\n",
    "        Tavg += (sim.T-Tavg)/(timestep+1)\n",
    "    else:\n",
    "        # equilibration stage\n",
    "        Pavg = sim.P\n",
    "        Tavg = sim.T\n",
    "\n",
    "   # print timestep, ke, pe, total e, Tavg, Pavg        \n",
    "    if timestep%50==0: # only print every 50 timesteps\n",
    "        # timestep, <ke>, <pe>, <e>, T, P, Tavg, Pavg\n",
    "        print(\"%4d  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\" % \n",
    "              (timestep,sim.ke/sim.N,sim.pe/sim.N,\n",
    "                (sim.ke+sim.pe)/sim.N,sim.T,sim.P,Tavg,Pavg))\n",
    "\n",
    "    pres.append(sim.P)\n",
    "    temp.append(sim.T)\n",
    "    time.append(timestep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3993efe-4dc7-4b09-9630-90201f989487",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A different way to calculate the average pressure and temperature\n",
    "# with standard deviations\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "print(\"Pressure:    {:.2f}+/-{:.2f}\"\n",
    "      .format(np.mean(pres[200:]),np.std(pres[200:])))\n",
    "print(\"Temperature: {:.2f}+/-{:.2f}\"\n",
    "      .format(np.mean(temp[200:]),np.std(temp[200:])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77d0de84-eb9b-4ea5-b5d3-4780335f3d7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot pressure and temperature\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(time[0:],pres[0:],label=\"P*\")\n",
    "ax.plot(time[0:],temp[0:],label=\"T*\")\n",
    "ax.set_xlabel(\"time step\")\n",
    "ax.legend(frameon=False)\n",
    "ax.set_xscale('linear')\n",
    "plt.title(\"Dimensionless instantaneous pressure and temperature\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}