{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Molecular dynamics\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>These examples *can* be run without ``asap3`` installed. In that case,\n  ASE\u2019s Python implementation of the EMT calculator can be used instead, but it\n  is much slower.</p></div>\n\n## Goal\n\nIn this tutorial, we will learn how to perform basic molecular dynamics (MD)\nsimulations using ASE.\n\nThe key objectives are:\n\n- Understand how to set up a crystal structure (Cu atoms on an FCC lattice).\n- Initialize velocities from Maxwell\u2013Boltzmann distribution corresponding to a\n  chosen temperature.\n- Integrate Newton\u2019s equations of motion using Velocity-Verlet algorithm and we\n  monitor the temperature using Langevin thermostat.\n- Monitor and analyze thermodynamic quantities (potential energy, kinetic\n  energy, total energy, temperature).\n- Save trajectories and visualize atomic motion with ASE\u2019s GUI.\n- Explore MD in different scenarios:\n  - Constant energy MD (NVE ensemble)\n  - Constant temperature MD (NVT ensemble)\n  - Isolated nanoparticle simulations\n\nBy the end of this tutorial, you should be able to set up your own MD\nsimulations, monitor energy conservation, and visualize system evolution.\n\n## Part 1: Basic Molecular Dynamics Simulation\n\nWe start by creating a copper crystal, assigning random velocities\ncorresponding to Maxwell Boltzmann Distribution at 300 K, and running dynamics\nin the NVE ensemble (constant energy).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\n# choose one of the following implementations of EMT:\n# included in ase\n# from ase.calculators.emt import EMT\n# faster performance\nfrom asap3 import EMT\n\nfrom ase import units\nfrom ase.cluster.cubic import FaceCenteredCubic as ClusterFCC\nfrom ase.io.trajectory import Trajectory\nfrom ase.lattice.cubic import FaceCenteredCubic as LatticeFCC\nfrom ase.md.langevin import Langevin  # for later NPT simulations\nfrom ase.md.velocitydistribution import (\n    MaxwellBoltzmannDistribution,\n    Stationary,\n    ZeroRotation,\n)\nfrom ase.md.verlet import VelocityVerlet\nfrom ase.optimize import QuasiNewton\nfrom ase.visualize.plot import plot_atoms\n\n# Set up initial positions of Cu atoms on Fcc crystal lattice\nsize = 10\natoms = LatticeFCC(\n    directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],\n    symbol='Cu',\n    size=(size, size, size),\n    pbc=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Before setting up the MD simulation, we take a look at the initial structure:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(5, 5))\nplot_atoms(atoms, ax, rotation=('45x,45y,0z'), show_unit_cell=2, radii=0.75)\nax.set_axis_off()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's run the MD simulation and monitor the kinetic and potential energy\nof the whole system:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Describe the interatomic interactions with the Effective Medium Theory (EMT)\natoms.calc = EMT()\n\n# Set the initial velocities corresponding to T=300K from Maxwell Boltzmann\n# Distribution\nMaxwellBoltzmannDistribution(atoms, temperature_K=300)\n\n# We use Velocity Verlet algorithm to integrate the Newton's equations.\ntimestep_fs = 5\ndyn = VelocityVerlet(atoms, timestep_fs * units.fs)  # 5 fs time step.\n\n\ndef printenergy(a):\n    \"\"\"\n    Function to print the thermodynamical properties i.e potential energy,\n    kinetic energy and total energy\n    \"\"\"\n    epot = a.get_potential_energy()\n    ekin = a.get_kinetic_energy()\n    temp = a.get_temperature()\n    print(\n        f'Energy per atom: Epot ={epot:6.3f}eV  Ekin = {ekin:.3f}eV '\n        f'(T={temp:.3f}K) Etot = {epot + ekin:.3f}eV'\n    )\n\n\n# Now run the dynamics\nprint('running a NVE simulation of fcc Cu')\nprintenergy(atoms)\n# init lists to for energy vs time data\ntime_ps, epot, ekin = [], [], []\nmdind = 0\nsteps_per_block = 10\nfor i in range(20):\n    dyn.run(steps_per_block)\n    mdind += steps_per_block\n    printenergy(atoms)\n    # save the energies of the current MD step\n    time_ps.append(mdind * timestep_fs / 1000.0)\n    epot.append(atoms.get_potential_energy())\n    ekin.append(atoms.get_kinetic_energy())\n\netot = np.array(epot) + np.array(ekin)\n\n# Plot energies vs time\nfig, ax = plt.subplots(figsize=(6, 4))\nax.plot(time_ps, epot, label='Potential energy')\nax.plot(time_ps, ekin, label='Kinetic energy')\nax.plot(time_ps, etot, label='Total energy')\nax.set_xlabel('Time (ps)')\nax.set_ylabel('Energy (eV)')\nax.legend(loc='best')\nax.grid(True, linewidth=0.5, alpha=0.5)\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note how the total energy is conserved, but the kinetic energy quickly\ndrops to half the expected value. Why?\n\nWhat you learned here:\n\n- How to set up a basic MD run.\n- How to monitor the energy over time.\n- That total energy is approximately conserved in NVE simulations, what is\n  the error in total energy?\n\nExercise: Tune the time step from 5fs to 10fs and 50fs, what changes do you\nobserve in total energy?\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Part 2: Constant temperature MD\n\nIn many cases, you want to control temperature (NVT ensemble). This\ncan be done using a thermostat, like -- in this tutorial -- Langevin\nthermostat.\nCompared to the previous example, we replace the line\n``dyn = VelocityVerlet(...)`` with\n\n::\n\n  dyn = Langevin(atoms, timestep=5 * units.fs, temperature_K=T,\n  friction=0.02)\n\nwhere ``T`` is the desired temperature in Kelvin. For that we also imported\nthe Langevin in the beginning.\n\nThe Langevin dynamics will then slowly adjust the total energy of the\nsystem so the temperature approaches the desired one.\n\nAs a slightly less boring example, let us use this to melt a chunk of\ncopper by starting the simulation without any momentum of the atoms\n(no kinetic energy), and with a desired temperature above the melting\npoint. We will also save information about the atoms in a trajectory\nfile called ``moldyn3.traj``.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>It is recommended to use the ``asap3`` implementation of the ``EMT``\n  calculator here, because its performance benefits over the ``ase``\n  implementation.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "size = 10\nT = 1500  # Kelvin\n\n# Set up a crystal\natoms = LatticeFCC(\n    directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],\n    symbol='Cu',\n    size=(size, size, size),\n    pbc=False,\n)\n\n# Describe the interatomic interactions with the Effective Medium Theory\natoms.calc = EMT()\n\n# We want to run MD with constant energy using the Langevin algorithm\n# with a time step of 5 fs, the temperature T and the friction\n# coefficient to 0.02 atomic units.\ntimestep_fs = 5\ndyn = Langevin(\n    atoms, timestep=timestep_fs * units.fs, temperature_K=T, friction=0.02\n)\n\n# We also want to save the positions of all atoms after every 100th time step.\ntraj = Trajectory('fccCu_NPT.traj', 'w', atoms)\n\n# Now run the dynamics\nprint('running a NVT simulation of fcc Cu')\nprintenergy(atoms)\ntime_ps, temperature = [], []\nmdind = 0\nsteps_per_block = 10\nfor i in range(200):\n    dyn.run(steps_per_block)\n    mdind += steps_per_block\n    printenergy(atoms)\n    # save the temperature of the current MD step\n    time_ps.append(mdind * timestep_fs / 1000.0)\n    temperature.append(atoms.get_temperature())\n\n# Plot temperatures vs time\nfig, ax = plt.subplots(figsize=(6, 4))\nax.plot(time_ps, temperature)\nax.set_xlabel('Time (ps)')\nax.set_ylabel('Temperature (K)')\nax.grid(True, linewidth=0.5, alpha=0.5)\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After running the simulation, you can study the result with the\ncommand\n\n::\n\n  ase gui fccCu_NPT.traj\n\nTry plotting the kinetic energy. Like in the temperature vs time plot you\nwill *not* see a well-defined melting point due to finite size effects\n(including surface melting),\nbut you will probably see an almost flat region where the inside of\nthe system melts. The outermost layers melt at a lower temperature.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The Langevin dynamics will by default keep the position and momentum\n  of the center of mass unperturbed. This is another improvement over\n  just setting momenta corresponding to a temperature, as we did before.</p></div>\n\n\n## Part 3: Isolated particle MD\n\nWhen simulating isolated particles with MD, it is sometimes preferable\nto set random momenta corresponding to a specific temperature and let the\nsystem evolve freely. With a relatively high temperature, the is however\na risk that the collection of atoms will drift out of the simulation box\nbecause the randomized momenta gave the center of mass a small but\nnon-zero velocity too.\n\nLet us see what happens when we propagate a nanoparticle:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "size = 4\natoms = ClusterFCC(\n    'Cu',\n    surfaces=[[1, 0, 0], [1, 1, 0], [1, 1, 1]],\n    layers=(size, size, size),\n    vacuum=4,\n)\n# asap3 requires a non-zero cell even if pbc are not applied\natoms.cell = [40] * 3\natoms.set_pbc(False)  # isolated cluster (explicit, for clarity)\n\n# Describe the interatomic interactions with the Effective Medium Theory\natoms.calc = EMT()\n\n# Quick relaxation of the cluster\nqn = QuasiNewton(atoms)\nqn.run(fmax=0.001, steps=10)\n\n# Set the momenta corresponding to T=1200 K\nMaxwellBoltzmannDistribution(atoms, temperature_K=1200)\nStationary(atoms)  # zero linear momentum\nZeroRotation(atoms)  # zero angular momentum\n\n# Run MD using the Velocity Verlet algorithm and save trajectory\ndyn = VelocityVerlet(atoms, 5 * units.fs, trajectory='nanoparticleCu_NVE.traj')\n\nprint('running a NVE simulation of a Cu nanoparticle')\nprintenergy(atoms)\nsteps_per_block = 10\nfor i in range(200):\n    dyn.run(steps_per_block)\n    printenergy(atoms)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After running the simulation, use `ase-gui` to compare the resulting\ntrajectory with how it looks if you comment out either the line that says\n``Stationary(atoms)``, ``ZeroRotation(atoms)`` or both:\n\n::\n\n  ase gui nanoparticleCu_NVE.traj\n\nTry playing the movie with a high frame rate and set frame skipping to a\nlow number. Can you spot the subtle difference?\n\n"
      ]
    }
  ],
  "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.13.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}