Trajectories

Trajectories are useful for storing the positions of the atoms in a long molecular dynamics run or in a structural optimization run. Here is a typical use case:

>>> from ASE.Trajectories.NetCDFTrajectory import NetCDFTrajectory
>>> from ASE.Dynamics.VelocityVerlet import VelocityVerlet
>>> from ASE.Calculators.PairPotential import PairPotential
>>> from ASE.IO.NetCDF import ReadNetCDF
>>> water = ReadNetCDF('water.nc')
>>> water.SetCalculator(PairPotential())
>>> dyn = VelocityVerlet(water, dt=0.05)
>>> traj = NetCDFTrajectory('traj.nc', water, interval=5)
>>> dyn.Attach(traj)
>>> dyn.Run(steps=100)

In this example, we are doing molecular dynamics using the Velocity Verlet algorithm. The line

>>> traj = NetCDFTrajectory('traj.nc', water, interval=5)

creates a trajectory object. The configurations are written to a netCDF file ('traj.nc').

The dyn.Attach(traj) command tells the dynamics object to notify the trajectory object after every move of the atoms (the trajectory is an observer, and the dynamics is a subject - see the observer pattern). The trajectory object was told to only write every fifth configuration to file (the interval=5 argument).

The dynamics object will call the Update() method of the trajectory - this method can also be called manually.

Opening a netCDF trajectory file

If no ListOfAtoms object is supplied in the NetCDFTrajectory constructor, the netCDF file is opened in read mode. If a ListOfAtoms object is pressent, write mode is used, and a new file is created (an existing file will be renamed to have a .bak extension). There is an optional mode argument with two possible values: 'new' will try to create a new netCDF file and report an error fi the file exists, and 'append' will append to an existing file.

Getting information out of a trajectory

For getting information out of a trajectory, the following methods can be used:

SetFrameNumber(frame):
Set the current frame number.
GetFrameNumber():
Get the current frame number.
GetListOfAtoms(frame=None):
Create a new ListOfAtoms object for configuration number frame (defaults to the current frame number).
Get(name, frame=None):
Get the named data for configuration number frame (defaults to the current frame number).

Examples:

>>> diff = (traj.Get('CartesianPositions', -1) -  # last
...         traj.Get('CartesianPositions', 0))    # first
>>> traj.Get('AtomicNumbers')
array([8, 1, 1])
>>> water2 = traj.GetListOfAtoms()

A netCDF trajectory can be opened in read only mode like this:

>>> traj = NetCDFTrajectory('traj.nc')

By default, the following quantities are allways put in a netCDF trajectory: 'AtomicNumbers', 'CartesianPositions', 'UnitCell' and 'Tags'. Those of 'PotentialEnergy', 'CartesianForces' and 'Stress' that are available are also put in. The tags and the atomic numbers are written once only. The default behavior can be changed with the Add(name) and Delete(name) methods:

>>> t = NetCDFTrajectory('traj.nc', water)
>>> t.Delete('CartesianForces')
>>> t.Add('UnitCell', once=True)

This trajectory will have no forces and the unit cell is only written once.

Note

The length and energy units used in the Python session that ceated a netCDF trajectory will always be included in the file. When using a trajectory in a Python session with different units, conversion will take place automatically.

Conversion has not been tested yet!

Adding user defined data

The Add(name, ...) method can be used to add almost anything to a trajectory. The add method has the following optional arguments:

data:
A callable returning the data in a numarray (or the actual data, in case the data does not change)
shape:
The shape of the numarray. Use the string 'natoms' for the number of atoms.
typecode:
The typecode for the numarray (num.Int, num.Float, ...).
once:
A boolean deciding whether the data should be written once or for every update.
units:
A tuple of two integers...

For the missing optional arguments, the Add() method will try to make good guesses for default values:

  • If name is one of the built-in quantities of the ListOfAtoms object, then everything is known.
  • If data is given as a numarray or a callable returning a numarray, then the shape and the typecode can be extracted. However, this will call the callable - this can be avoided by supplying the shape and typecode.
  • once will default to False - write every time.
  • units will default to (0, 0) - a pure number.

Suppose you want to include one integer for each atom, once:

>>> data = num.array([1, 2, 3])
>>> t = NetCDFTrajectory('traj.nc', water)
>>> t.Add('Stuff', data, once=True)

Or if you have a function that calculates something from a molecule (something = calculate_something(water)):

>>> def data():
...     return calculate_something(water)
...
>>> t.Add('Something', data)

Note

Float32 <--> Float64 does not work!

The plottrajectory tool

plottrajectory will open a rasmol window for wieving the atomic structure, a Gnuplot window for showing xy data and a Tk window for changing the current frame number. By default the Gnuplot window will show the total energy as a functions of the trajectory frame number.

The plottrajectory tool is called from the command line and the syntax is:

plottrajectory [-r R1 R2 R3] [-u usermodule]  [-p "datamethods"] trajectory

options

--version show programs version number and exit
-h, --help show this help message and exit
-r R1R2R3 repeat R1, R2, R3 times along the three axes
-p listofplots adding a gnuplot window
-u <usermodule>
 add methods from module

By using the -r option you can repeat the atoms R<x> times the <x>-axis.

By using the -p option you can plot xy data using the gnuplot window, default being (frame number,total energy) plot. The syntax is:

-p "[list,list,..]"

list will here correspond to a separate gnuplot window. Each list can have one or more elements that will be added to the gnuplot window, some examples:

[[x]] plot the method x against the frame number
[[x,y]] plot the method y against the method y.
[[[x,y],[x,y]]] two plots in the same gnuplot window.
[[x,y],[x,y]] two gnuplot windows.

By using the -u <usermodule> option you can import methods defined in your own module, the method must be functions or callable classes, taking a listofatoms object as argument.

An example could be the file distance.py calculating the distance between two atoms:

import Numeric as num
class Distance:
    """Returns the distance between atom number
       ``a`` and ``b``."""
    def __init__(self, a, b):
        self.a = a
        self.b = b
    def __call__(self, atoms):
        pos = atoms.GetCartesianPositions()
        d = pos[self.a] - pos[self.b]
        return num.sqrt(num.dot(d, d))
    def __repr__(self):
        return 'Distance bewteen atom %d and %d' % \
               (self.a, self.b)

Now you can test the method using:

plottrajectory -p "[[Distance(0,1)]]" -u distance.py mytraj.nc

If you think that your method could be of use to others, please consider sending it to the mailing list, and we will add it to the utilities allready defined.

At the moment the following methods are defined:

Geometry:
Position(atom-number,coordinate): correspond to the ListOfAtoms
method atoms.GetCartesianPositions()[atom-number, coodinate]
Force(atom-number,coordinate): correspond to the ListOfAtoms
method atoms.GetCartesianForces()[atom-number,coordinate]

MaxForce: The magnitude of then largest force acting on any atom.

Volume: Volume of the unitcell.

Distance(atom1,atom2): The distance between atom1 and atom2.

Molecular-Dynamics:

Temperature: The temperature of the atoms.

PotentialEnergy: The total potential energy of the atoms.

So if you want to plot the Total energy as a function of some distance between two atoms you can use:

-p "[[Distance(0,1),TotalEnergy]]"

or if you want to plot the temperature and total energy in two plots you can do:

-p "[[TotalEnergy], [Temperature]]"

ase2: Trajectories (last edited 2010-10-20 09:11:15 by localhost)