.. contents::
   :local:

------------
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 <Dynamics#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 <Dynamics#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 <Observer#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]]"