.. 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 `_ 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 add methods from module By using the ``-r`` option you can repeat the atoms R times the -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 `` 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]]"