{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Equilibrating an MD box of acetonitrile\n\n# Goals\n\nIn this tutorial we learn how to perform a thermal equilibration of a box\nof acetonitrile molecules using ASE. We will:\n\n* build a coarse-grained CH\u2083\u2013C\u2261N triatomic model.\n* set up a periodic box at the experimental density (298 K).\n* apply rigid-molecule constraints.\n* use the ACN force field.\n* equilibrate with Langevin dynamics.\n* scale small box of 27 molecules to 216 molecules.\n\n# ACN model\n\nThe acetonitrile force field implemented in ASE (:mod:`ase.calculators.acn`)\nis an interaction potential between three-site linear molecules. The atoms\nof the methyl group are treated as a single site centered on the methyl\ncarbon (hydrogens are not explicit). Therefore:\n\n* assign the **methyl mass** to the outer carbon (``m_me``),\n* use the atomic sequence **Me\u2013C\u2013N** repeated for all molecules,\n* keep molecules **rigid** during MD with :class:`FixLinearTriatomic`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport ase.units as units\nfrom ase import Atoms\nfrom ase.calculators.acn import ACN, m_me, r_cn, r_mec\nfrom ase.constraints import FixLinearTriatomic\nfrom ase.io import Trajectory\nfrom ase.md import Langevin\nfrom ase.visualize.plot import plot_atoms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 1: molecule\nBuild one CH3\u2013C\u2261N as a linear triatomic \u201cC\u2013C\u2013N\u201d. The first carbon is\nthe methyl site; we assign it the CH3 mass. Rotate slightly to avoid\nperfect alignment with the cell axes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pos = [[0, 0, -r_mec], [0, 0, 0], [0, 0, r_cn]]\natoms = Atoms('CCN', positions=pos)\natoms.rotate(30, 'x')\n\nmasses = atoms.get_masses()\nmasses[0] = m_me\natoms.set_masses(masses)\n\nmol = atoms.copy()\nmol.set_pbc(False)\nmol.set_cell([20, 20, 20])\nmol.center()\n\nfig, ax = plt.subplots(figsize=(4, 4))\nplot_atoms(\n    mol,\n    ax=ax,\n    rotation='30x,30y,0z',\n    show_unit_cell=0,\n    radii=0.75,\n)\nax.set_axis_off()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 2: Set up small box of 27-molecules\nMatch density 0.776 g/cm^3 at 298 K. Compute cubic box length L from\nmass and density. Build 3\u00d73\u00d73 supercell and enable PBC.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "density = 0.776 / 1e24  # g / \u00c5^3\nL = ((masses.sum() / units.mol) / density) ** (1 / 3.0)\n\natoms.set_cell((L, L, L))\natoms.center()\natoms = atoms.repeat((3, 3, 3))\natoms.set_pbc(True)\n\nbox27 = atoms.copy()\nfig, ax = plt.subplots(figsize=(5, 5))\nplot_atoms(\n    box27,\n    ax=ax,\n    rotation='35x,35y,0z',\n    show_unit_cell=2,\n    radii=0.75,\n)\nax.set_axis_off()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 3: Set constraints\nKeep each \u201cC\u2013C\u2013N\u201d rigid during MD using FixLinearTriatomic.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nm = 27\ntriples = [(3 * i, 3 * i + 1, 3 * i + 2) for i in range(nm)]\natoms.constraints = FixLinearTriatomic(triples=triples)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 4: MD run for 27-molecules system\nAssign ACN with cutoff = half the smallest box edge. Langevin MD at\n300 K, 1 fs timestep. Save a frame every step.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2)\n\ntag = 'acn_27mol_300K'\nmd = Langevin(\n    atoms,\n    1 * units.fs,\n    temperature_K=300,\n    friction=0.01,\n    logfile=tag + '.log',\n)\ntraj = Trajectory(tag + '.traj', 'w', atoms)\nmd.attach(traj.write, interval=10)\nmd.run(5000)  # 5 ps @ 1 fs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 5: scale system size to 216 molecules\nRepeat 2\u00d72\u00d72 to reach 216 molecules. Reapply constraints and update\nthe ACN cutoff for the new cell.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms.set_constraint()\natoms = atoms.repeat((2, 2, 2))\n\nnm = 216\ntriples = [(3 * i, 3 * i + 1, 3 * i + 2) for i in range(nm)]\natoms.constraints = FixLinearTriatomic(triples=triples)\n\natoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 6: MD run for 216-molecules system\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tag = 'acn_216mol_300K'\nmd = Langevin(\n    atoms,\n    2 * units.fs,\n    temperature_K=300,\n    friction=0.01,\n    logfile=tag + '.log',\n)\ntraj = Trajectory(tag + '.traj', 'w', atoms)\nmd.attach(traj.write, interval=10)\n\ntimes_ps, epots, ekins, etots, temps = [], [], [], [], []\nsample_interval = 10  # sample every 10 MD steps for lighter plots\n\n\ndef sample():\n    # Time in ps (same as MDLogger: dyn.get_time() / (1000 * units.fs))\n    t_ps = md.get_time() / (1000.0 * units.fs)\n    ep = atoms.get_potential_energy()  # eV total\n    ek = atoms.get_kinetic_energy()  # eV total\n    T = atoms.get_temperature()  # K\n    times_ps.append(t_ps)\n    epots.append(ep)\n    ekins.append(ek)\n    etots.append(ep + ek)\n    temps.append(T)\n\n\n# initial sample at t=0\nsample()\nmd.attach(sample, interval=sample_interval)\nmd.run(1000)  # 6 ps @ 2 fs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot Instantaneous temperature vs time.\nDoes the system equilibrated well?\nWhat is the average temperature? Should we run longer simulations?\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(6, 4))\nax.plot(times_ps, temps, label='T (K)')\nax.set_xlabel('Time (ps)')\nax.set_ylabel('Temperature (K)')\nax.legend(loc='best')\nax.grid(True, linewidth=0.5, alpha=0.5)\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Next steps\n* View trajectories::\n\n    ase gui acn_27mol_300K.traj\n\n    ase gui acn_216mol_300K.traj\n\n* Plot other thermodynamic quantities (p.e., k.e., and total energy).\n\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
}