{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Equilibrating a TIPnP Water Box\n\nThis tutorial shows how to use the TIP3P and TIP4P force fields in\nASE.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>Due to GitLab limit we have to cut the simulation steps in the tutorial.\n  Please adjust the number of steps to your own system.\n  We recommend multiply the number of steps we provided by 1000 for\n  a realistic use.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import ase.units as units\nfrom ase.build import molecule\nfrom ase.calculators.tip3p import TIP3P\nfrom ase.constraints import FixBondLengths\nfrom ase.io.trajectory import Trajectory\nfrom ase.md import Langevin"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first create a starting point of the simulaiton.\nWe will create a water box at 20 \u00b0C density.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms = molecule('H2O')\ndensity = 0.9982  # denisty of water in g/cm^3 at 20\u00b0C\nbox_length = ((atoms.get_masses().sum() / units.mol) / (density * 1e-24)) ** (\n    1 / 3\n)  # box length in \u00c5\natoms.set_cell((box_length, box_length, box_length))\natoms.center()\n# Repeat the water molecule we just created to end up with a PBC cell\natoms *= (3, 3, 3)\natoms.set_pbc(True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualise the starting box\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nfrom ase.visualize.plot import plot_atoms\n\nfig, ax = plt.subplots()\nplot_atoms(atoms, ax)\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since the TIPnP type water interpotentials are for rigid\nmolecules, there are no intramolecular force terms, and we need to\nconstrain all internal degrees of freedom. For this, we're\nusing the RATTLE-type constraints of the `FixBondLengths` class to\nconstrain all internal atomic distances (O-H1, O-H2, and H1-H2) for\neach molecule.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms.constraints = FixBondLengths(\n    [(3 * i + j, 3 * i + (j + 1) % 3) for i in range(3**3) for j in [0, 1, 2]]\n)\n# RATTLE-type constraints on O-H1, O-H2, H1-H2.\ntag = 'tip3p_27mol_equil'\natoms.calc = TIP3P(rc=4.5)  # set the calculator to be the TIP3P force field"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For efficiency, we first equillibrate a smaller box, and then repeat that\nonce more for the final equillibration. However, the potentials are not\nparallelized, and are mainly included for testing and for use with QM/MM\ntasks, so expect to let it run for some time.\nFor illustration, we will equillibrate our system with the Langevin\nthermostat.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Equillibrate in a small box first.\nmd = Langevin(\n    atoms,\n    1 * units.fs,\n    temperature_K=293.15,  # 20 \u00b0C\n    friction=0.01,\n    logfile=tag + '.log',\n)\n\ntraj = Trajectory(tag + '.traj', 'w', atoms)\nmd.attach(traj.write, interval=1)\nmd.run(4)  # please use 4000 to better equilibrate\n\n# Repeat box and equilibrate further.\ntag = 'tip3p_216mol_equil'\natoms.set_constraint()  # repeat not compatible with FixBondLengths currently.\natoms *= (2, 2, 2)\natoms.constraints = FixBondLengths(\n    [\n        (3 * i + j, 3 * i + (j + 1) % 3)\n        for i in range(len(atoms) // 3)\n        for j in [0, 1, 2]\n    ]\n)\natoms.calc = TIP3P(rc=7.0)\nmd = Langevin(\n    atoms,\n    2 * units.fs,\n    temperature_K=293.15,  # 20 \u00b0C\n    friction=0.01,\n    logfile=tag + '.log',\n)\n\ntraj = Trajectory(tag + '.traj', 'w', atoms)\nmd.attach(traj.write, interval=1)\nmd.run(2)  # please use 2000 to better equilibrate"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p></p></div>\n\n The temperature calculated by ASE is assuming all degrees of freedom\n are available to the system. Since the constraints have removed the 3\n vibrational modes from each water, the shown temperature will be 2/3\n of the actual value.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The procedure for the TIP4P force field is the same, with the following\nexception: the atomic sequence **must** be OHH, OHH, ... .\nSo to perform the same task using TIP4P, you simply have to import\nthat calculator instead:\n``from ase.calculators.tip4p import TIP4P, rOH, angleHOH``\n\nMore info about the TIP4P potential: :mod:`ase.calculators.tip4p`\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
}