Equilibrating an MD box of acetonitrile

In this tutorial we see how to perform a thermal equilibration of an MD box of classical acetonitrile molecules using the Langevin module and the implementation of an acetonitrile force field in ASE.

The acetonitrile force field implemented in ASE (ase.calculators.acn) is an interaction potential between three-site linear molecules, in which the atoms of the methyl group are treated as a single site centered on the methyl carbon, i.e. hydrogens are not considered explicitly. For this reason, while setting up a box of acetonitrile one has to assign the mass of a methyl to the outer carbon atom. The calculator requires the atomic sequence to be MeCN … MeCN or NCMeNCMe … NCMe, where Me represents the methyl site.

As for the TIPnP models, the acetonitrile potential works with rigid molecules. However, due to the linearity of the acetonitrile molecular model, we cannot fix the geometry by constraining all interatomic distances using FixBondLengths, as is done for TIPnP water. Instead, we must use the class FixLinearTriatomic

The MD procedure we use for the equilibration closely follows the one presented in the tutorial Equilibrating a TIPnP Water Box.

import numpy as np

import ase.units as units
from ase import Atoms
from ase.calculators.acn import ACN, m_me, r_cn, r_mec
from ase.constraints import FixLinearTriatomic
from ase.io import Trajectory
from ase.md import Langevin

pos = [[0, 0, -r_mec],
       [0, 0, 0],
       [0, 0, r_cn]]
atoms = Atoms('CCN', positions=pos)
atoms.rotate(30, 'x')

# First C of each molecule needs to have the mass of a methyl group
masses = atoms.get_masses()
masses[::3] = m_me
atoms.set_masses(masses)

# Determine side length of a box with the density of acetonitrile at 298 K
# Density in g/Ang3 (https://pubs.acs.org/doi/10.1021/je00001a006)
d = 0.776 / 1e24
L = ((masses.sum() / units.mol) / d)**(1 / 3.)
# Set up box of 27 acetonitrile molecules
atoms.set_cell((L, L, L))
atoms.center()
atoms = atoms.repeat((3, 3, 3))
atoms.set_pbc(True)

# Set constraints for rigid triatomic molecules
nm = 27
atoms.constraints = FixLinearTriatomic(
    triples=[(3 * i, 3 * i + 1, 3 * i + 2)
             for i in range(nm)])

tag = 'acn_27mol_300K'
atoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2)

# Create Langevin object
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(5000)

# Repeat box and equilibrate further
atoms.set_constraint()
atoms = atoms.repeat((2, 2, 2))
nm = 216
atoms.constraints = FixLinearTriatomic(
    triples=[(3 * i, 3 * i + 1, 3 * i + 2)
             for i in range(nm)])

tag = 'acn_216mol_300K'
atoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2)

# Create Langevin object
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(3000)