Skip to main content
`Ctrl`+`K`

# 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.0)
    # 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)