A simple example of an Asap simulation

The following script sets up a 500-atom cluster of atoms in the shape of a small crystallite, and then performs a constant-temperature molecular dynamics simulation at 800K using the Langevin algorithm. After a while, the temperature is increased to 2000K. The potential, kinetic and total energies are printed every 100 timestep, and the atoms are written to a trajectory file every 2000 timesteps:

# First a documentation string to describe the example
"""An example script demonstrating an NVT MD simulation with Asap.

This example runs a molecular dynamics simulation on a small copper
cluster in the NVT ensemble, i.e. with constant particle number,
constant volume and constant temperature.  It uses Langevin dynamics
to get the constant temperature.
"""

# Import ASAP and relatives
from asap3 import EMT, units, Trajectory
from asap3.md.langevin import Langevin
from ase.lattice.cubic import FaceCenteredCubic
from ase.visualize import view
from asap3.md.velocitydistribution import MaxwellBoltzmannDistribution

# Import other essential modules
import numpy as np

# Set up an fcc crystal of copper, 1372 atoms.
atoms = FaceCenteredCubic(size=(7,7,7), symbol="Cu", pbc=False)

# Now the standard EMT Calculator is attached
atoms.calc = EMT()

# Make an object doing Langevin dynamics at a temperature of 800 K
dyn = Langevin(atoms, timestep=5*units.fs, temperature_K=800,
               friction=0.005)

# Set the momenta corresponding to T=1600K.  The temperature will
# quickly drop to half of that as the energy is distributed evenly
# among the kinetic and potential energy.
MaxwellBoltzmannDistribution(atoms, temperature_K=1600)

# Make a trajectory
traj = Trajectory('MD_Cluster.traj', "w", atoms)
dyn.attach(traj, interval=500) # Automatically writes the initial configuration

# Print the energies
_step = 0
def printenergy(a):
    global _step
    n = len(a)
    ekin = a.get_kinetic_energy() / n
    epot = a.get_potential_energy() / n
    print ("%4d: E_kin = %-9.5f  E_pot = %-9.5f  E_tot = %-9.5f  T = %.1f K" %
           (_step, ekin, epot, ekin+epot, atoms.get_temperature()))
    _step += 1

printenergy(atoms)
view(atoms)  #If computer is busy, the plot will appear with some delay.

# Now do the dynamics, doing 5000 timesteps, writing energies every 50 steps
dyn.attach(printenergy, 50, atoms)
dyn.run(5000)
view(atoms)

# Now increase the temperature to 2000 K and continue
dyn.set_temperature(temperature_K=2000)
dyn.run(5000)
view(atoms)

You can also download MD_Cluster3.py