Molecular dynamics

Note

These examples can be used without Asap installed, then the ase.EMT calculator (implemented in Python) is used, but nearly superhuman patience is required.

Here we demonstrate now simple molecular dynamics is performed. A crystal is set up, the atoms are given momenta corresponding to a temperature of 300K, then Newtons second law is integrated numerically with a time step of 5 fs (a good choice for copper).

"""Demonstrates molecular dynamics with constant energy."""

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units

# Use Asap for a huge performance increase if it is installed
use_asap = False

if use_asap:
    from asap3 import EMT
    size = 10
else:
    from ase.calculators.emt import EMT
    size = 3

# Set up a crystal
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol='Cu',
                          size=(size, size, size),
                          pbc=True)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.calc = EMT()

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, temperature_K=300)

# We want to run MD with constant energy using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.


def printenergy(a):
    """Function to print the potential, kinetic and total energy"""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))


# Now run the dynamics
printenergy(atoms)
for i in range(20):
    dyn.run(10)
    printenergy(atoms)

Note how the total energy is conserved, but the kinetic energy quickly drops to half the expected value. Why?

Instead of printing within a loop, it is possible to use an “observer” to observe the atoms and do the printing (or more sophisticated analysis).

"""Demonstrates molecular dynamics with constant energy."""

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units

# Use Asap for a huge performance increase if it is installed
use_asap = True

if use_asap:
    from asap3 import EMT
    size = 10
else:
    from ase.calculators.emt import EMT
    size = 3

# Set up a crystal
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol="Cu",
                          size=(size, size, size),
                          pbc=True)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.calc = EMT()

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, temperature_K=300)

# We want to run MD with constant energy using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.


def printenergy(a=atoms):  # store a reference to atoms in the definition.
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))


# Now run the dynamics
dyn.attach(printenergy, interval=10)
printenergy()
dyn.run(200)

Constant temperature MD

Often, you want to control the temperature of an MD simulation. This can be done with the Langevin dynamics module. In the previous examples, replace the line dyn = VelocityVerlet(...) with:

dyn = Langevin(atoms, 5*units.fs, T*units.kB, 0.002)

where T is the desired temperature in Kelvin. You also need to import Langevin, see the class below.

The Langevin dynamics will then slowly adjust the total energy of the system so the temperature approaches the desired one.

As a slightly less boring example, let us use this to melt a chunk of copper by starting the simulation without any momentum of the atoms (no kinetic energy), and with a desired temperature above the melting point. We will also save information about the atoms in a trajectory file called moldyn3.traj.

"""Demonstrates molecular dynamics with constant temperature."""
from ase.lattice.cubic import FaceCenteredCubic
from ase.md.langevin import Langevin
from ase.io.trajectory import Trajectory
from ase import units

from asap3 import EMT  # Way too slow with ase.EMT !
size = 10

T = 1500  # Kelvin

# Set up a crystal
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol="Cu",
                          size=(size, size, size),
                          pbc=False)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.calc = EMT()

# We want to run MD with constant energy using the Langevin algorithm
# with a time step of 5 fs, the temperature T and the friction
# coefficient to 0.02 atomic units.
dyn = Langevin(atoms, 5 * units.fs, T * units.kB, 0.002)


def printenergy(a=atoms):  # store a reference to atoms in the definition.
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))


dyn.attach(printenergy, interval=50)

# We also want to save the positions of all atoms after every 100th time step.
traj = Trajectory('moldyn3.traj', 'w', atoms)
dyn.attach(traj.write, interval=50)

# Now run the dynamics
printenergy()
dyn.run(5000)

After running the simulation, you can study the result with the command

ase gui moldyn3.traj

Try plotting the kinetic energy. You will not see a well-defined melting point due to finite size effects (including surface melting), but you will probably see an almost flat region where the inside of the system melts. The outermost layers melt at a lower temperature.

Note

The Langevin dynamics will by default keep the position and momentum of the center of mass unperturbed. This is another improvement over just setting momenta corresponding to a temperature, as we did before.

Isolated particle MD

When simulating isolated particles with MD, it is sometimes preferable to set random momenta corresponding to a specific temperature and let the system evolve freely. With a relatively high temperature, the is however a risk that the collection of atoms will drift out of the simulation box because the randomized momenta gave the center of mass a small but non-zero velocity too.

Let us see what happens when we propagate a nanoparticle for a long time:

"""Demonstrates molecular dynamics for isolated particles."""
from ase.cluster.cubic import FaceCenteredCubic
from ase.optimize import QuasiNewton
from ase.md.velocitydistribution import (MaxwellBoltzmannDistribution,
                                         Stationary, ZeroRotation)
from ase.md.verlet import VelocityVerlet
from ase import units

# Use Asap for a huge performance increase if it is installed
use_asap = True

if use_asap:
    from asap3 import EMT
    size = 4
else:
    from ase.calculators.emt import EMT
    size = 2

# Set up a nanoparticle
atoms = FaceCenteredCubic('Cu',
                          surfaces=[[1, 0, 0], [1, 1, 0], [1, 1, 1]],
                          layers=(size, size, size),
                          vacuum=4)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.calc = EMT()

# Do a quick relaxation of the cluster
qn = QuasiNewton(atoms)
qn.run(0.001, 10)

# Set the momenta corresponding to T=1200K
MaxwellBoltzmannDistribution(atoms, temperature_K=1200)
Stationary(atoms)  # zero linear momentum
ZeroRotation(atoms)  # zero angular momentum

# We want to run MD using the VelocityVerlet algorithm.

# Save trajectory:
dyn = VelocityVerlet(atoms, 5 * units.fs, trajectory='moldyn4.traj')


def printenergy(a=atoms):  # store a reference to atoms in the definition.
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))


dyn.attach(printenergy, interval=10)

# Now run the dynamics
printenergy()
dyn.run(2000)

After running the simulation, use ASE’s GUI to compare the results with how it looks if you comment out either the line that says \(Stationary(atoms)\), \(ZeroRotation(atoms)\) or both.

ase gui moldyn4.traj

Try playing the movie with a high frame rate and set frame skipping to a low number. Can you spot the subtle difference?