Equilibrating a TIPnP Water Box

This tutorial shows how to use the TIP3P and TIP4P force fields in ASE.

Since the TIPnP type water interpotentials are for rigid molecules, there are no intramolecular force terms, and we need to constrain all internal degrees of freedom. For this, we’re using the RATTLE-type constraints of the The FixBondLengths class class to constrain all internal atomic distances (O-H1, O-H2, and H1-H2) for each molecule.

The box is equillibrated with the Langevin thermostat.

For efficiency, we first equillibrate a smaller box, and then repeat that once more for the final equillibration. However, the potentials are not parallelized, and are mainly included for testing and for use with QM/MM tasks, so expect to let it run for some time.

The following is for TIP3P:

from ase import Atoms
from ase.constraints import FixBondLengths
from ase.calculators.tip3p import TIP3P, rOH, angleHOH
from ase.md import Langevin
import ase.units as units
from ase.io.trajectory import Trajectory
import numpy as np


# Set up water box at 20 deg C density
x = angleHOH * np.pi / 180 / 2
pos = [[0, 0, 0],
       [0, rOH * np.cos(x), rOH * np.sin(x)],
       [0, rOH * np.cos(x), -rOH * np.sin(x)]]
atoms = Atoms('OH2', positions=pos)

vol = ((18.01528 / 6.022140857e23) / (0.9982 / 1e24))**(1 / 3.)
atoms.set_cell((vol, vol, vol))
atoms.center()

atoms = atoms.repeat((3, 3, 3))
atoms.set_pbc(True)

# RATTLE-type constraints on O-H1, O-H2, H1-H2.
atoms.constraints = FixBondLengths([(3 * i + j, 3 * i + (j + 1) % 3)
                                    for i in range(3**3)
                                    for j in [0, 1, 2]])

tag = 'tip3p_27mol_equil'
atoms.calc = TIP3P(rc=4.5)
md = Langevin(atoms, 1 * units.fs, temperature=300 * units.kB,
              friction=0.01, logfile=tag + '.log')

traj = Trajectory(tag + '.traj', 'w', atoms)
md.attach(traj.write, interval=1)
md.run(4000)

# Repeat box and equilibrate further.
tag = 'tip3p_216mol_equil'
atoms.set_constraint()  # repeat not compatible with FixBondLengths currently.
atoms = atoms.repeat((2, 2, 2))
atoms.constraints = FixBondLengths([(3 * i + j, 3 * i + (j + 1) % 3)
                                    for i in range(len(atoms) / 3)
                                    for j in [0, 1, 2]])
atoms.calc = TIP3P(rc=7.)
md = Langevin(atoms, 2 * units.fs, temperature=300 * units.kB,
              friction=0.01, logfile=tag + '.log')

traj = Trajectory(tag + '.traj', 'w', atoms)
md.attach(traj.write, interval=1)
md.run(2000)

Note

The temperature calculated by ASE is assuming all degrees of freedom are available to the system. Since the constraints have removed the 3 vibrational modes from each water, the shown temperature will be 2/3 of the actual value.

The procedure for the TIP4P force field is the same, with the following exception: the atomic sequence must be OHH, OHH, … .

So to perform the same task using TIP4P, you simply have to import that calculator instead:

from ase.calculators.tip4p import TIP4P, rOH, angleHOH

More info about the TIP4P potential: ase.calculators.tip4p