Dynamics object

There are two types of dynamics objects: Objects for structure optimization and objects for various kinds of molecular dynamics. Both types of dynamics works with a ListOfAtoms (possibly through a filter) attached to a calculator. Let us look at an example: Molecular dynamics with the water molecule from the last section.

>>> from ASE.Dynamics.VelocityVerlet import VelocityVerlet
>>> dyn = VelocityVerlet(water, dt=0.05)
>>> for i in range(10):
...     pot = water.GetPotentialEnergy()
...     kin = water.GetKineticEnergy()
...     print pot + kin, pot, kin
...     dyn.Run(steps=10)
...
-9.66399870713 -9.66399870713 0.0
-9.66400132881 -9.66413187948 0.000130550671208
-9.66400554131 -9.66434484148 0.000339300168651
-9.66400543329 -9.66433839178 0.000332958486978
-9.66400111904 -9.66412094698 0.000119827942723
-9.66399869087 -9.66399910839 4.17518564929e-07
-9.6640016008 -9.66414515729 0.000143556487225
-9.66400571344 -9.6643526226 0.000346909165177
-9.66400522718 -9.66432906126 0.000323834081303
-9.66400086336 -9.66410826704 0.000107403681751

The dyn object gets the water molecule and a time step for the integration of Newtons equation. The Run(steps=10) call takes 10 steps in the Velocity Verlet algorithm.

All dynamics objects have a Run() method. This is because they are all derived from a common base class called Dynamics. The Dynamics class defines the following interface:

Run(steps=10**10, maxtime=1e1000):
Take steps steps, but stop after maxtime seconds.
Step():
Take one step.

Structure optimization

Dynamics objects used for structure optimization are all derived from a Minimizer base class, that again is derived from the Dynamics class. A minimizer will only use a subset of the ListOfAtoms interface: GetGeneralizedForces(), GetGeneralizedPositions(), SetGeneralizedPositions() and perhaps GetPotentialEnergy().

A minimizer has the Run() and Step() methods and also the following two methods:

Converged():
Return True if converged, othervice return False.
Converge(maxsteps=10**10, maxtime=1e1000):
Run until convergence, but no more than maxsteps steps and stop after maxtime seconds.

The convergence criteria can be set when constucting the minimizer usen one of these keywords:

fmax=maximum_force:
The absolute magnitude of the force should be less than maximum_force for all atoms.
converged=function:
The function argument should be a callable that returns True or False meaning converged or not converged.

(see next paragraph for examples).

MDMin

MDMin is a minimizer using the simple but efficient Molecular Dynamics minimization algorithm, sometimes known as QuickMin.

Minimizing the energy of the water molecule can be done using the MDMin minimizer:

>>> m = MDMin(water, fmax=0.05)
>>> m.Converge()
>>> print water.GetCartesianForces()

or equivalently, but unnecessarily complicated:

>>> def f():
...     forces = water.GetCartesianForces()
...     for F in forces:
...         if F[0]**2 + F[1]**2 + f[2]**2 > 0.05**2:
...             return False
...     return True
...
>>> m = MDMin(water, converged=f)
>>> m.Converge()

The MDmin algorithm is a modification of the usual VelocityVerlet molecular dynamics algorithm. Newtons second law is solved numerically, but after each time step the dot product between the forces and the momenta is checked. If it is zero, the system has just passed through a (local) minimum in the potential energy, the kinetic energy is large and about to decrease again. At this point, the momentum is set to zero. Unlike a "real" molecular dynamics, the masses of the atoms are not used, instead all masses are set to one.

The MDmin algorithm exists in two flavors, one where each atom is tested and stopped individually (QuickMinAtomByAtom in the old ASE), and one where all coordinates are treated as one long vector, and all momenta are set to zero if the dotproduct between the momentum vector and force vector (both of length 3N) is zero (QuickMinAllCoordinates in the old ASE). This module implements the latter version.

Although the algorithm is primitive, it performs very well because it takes advantage of the physics of the problem. Once the system is so near the minimum that the potential energy surface is approximately quadratic it becomes advantageous to switch to a minimization method with quadratic convergence, such as Conjugate Gradient or Quasi Newton.

Conjugate Gradient

Another algorithm for minimization is the well-known Conjugate Gradient algorithm (see, for example, Numerical Recipes in C, by Press, Teukolsky, Vetterling and Flannery).

>>> cg = ConjugateGradient(water, fmax=0.05)
>>> cg.Converge(maxsteps=100)

The method used to perform 1D (line-) minimizations along a given search direction is specificable by passing an appropriate LineMinimizer object to the constructor.

>>> br = BrentLineMinimizer(0.1) # length scale for initial step size
>>> cg = ConjugateGradient(water, lineMin=br,fmax=0.05)
>>> cg.Converge(maxsteps=100)

There are currently two LineMinimizer classes. BrentLineMinimizer uses the Brent algorithm described in Numerical Recipes. LM1, written by Jens Jorgen Mortensen, uses a more agressive approach to reduce the number of function/gradient evaluations per line search. LM1 is the defalt LineMinimizer.

Quasi Newton

A standard reference to this method is given in reference [1].

This implemetation is a translation of a fortran implementation of this method by Frank Jensen. Thanks to Frank Jensen for allowing us to use his code and for help in translating it into python.

The QuasiNewton class implements a minimizer, plus a saddle point/transition state search algorithm.

The current quasi-Newton algorithm uses direct diagonalization of the Hessian matrix, and this will limit the number of dynamics atoms that can be handled.

The QuasiNewton method can be used like this:

from ASE.Dynamics.QuasiNewton import QuasiNewton
qn = QuasiNewton(water, fmax=0.05)
qn.Converge(maxsteps=100)

The algorithm has the following parameters, in addtion to the standard method of a mimimizer:

radius:
The initial maximum length of any displacement vector, default is 0.05*sqrt(number of atoms). radius will be adjusted during the mimimization. It is often a good idea specifying a small initial value for radius, this will allow the algorithm to get a better estimate for the Hessian, before any large steps are attempted.
maxradius:
The maximum value radius can reach, default 0.05*sqrt(number of atoms).
hessian:
Using this keyword an initial hessian can be specified, otherwise a diagonal hessian is used as a staring guess.
diagonal:
The default hessian is a diagonal matrix, with elements given by diagonal, default is 50. For the initial steps this will correspond to the inverse of a step length, so setting a high value for diagonal`, will give rather small steps until the Hessian is updated. On the other hand setting a low value will give large initial steps. If the steps are too large, and the energy is increased, the steps is rejected, and the current trust radius will be reduced ontil the step is downhill.
hessianupdate:

Three methods for updating the Hessian matrix is defined:

  • BFGS (default)
  • Powell
  • Bofill

See more details in reference [1].

logfilename:

Set the filename for the log file, default qn.log. The keyword QN can be used for an easy overview of the minimization:

grep QN qn.log

This will look like:

QN:  quasi-Newton minimization started Fri Jul 21 09:50:53 2006
QN:    iter    time        energy      max force  radius
QN:      1    9:50:53    -463.257949     3.041     0.122
QN:      2    9:53:00    -463.423123     0.851     0.122
QN:      3    9:54:59    -463.443752     0.185     0.087
QN:      4    9:56:44    -463.445170     0.042     0.122

The last column is the current trust radius, which are updated during the minimization.

restartfile:
For each step QuasiNewton will write a file that saves the complete state of the algorithm, including Hessian matrix, current radius, allowing it to be restarted efficiently. If restartfile is set to an existing file, the algorithm will continue using the information in the file. Note this should only be done using a restart file from the previous iterations, otherwise the algorithm will be confused. If restartfile is set to a non existing file, the algorithm will use this name for the restart file. Default name for the restart file is qn-restart.pickle.
transitionstate:
If transitionstate=True the algorithm will search for a transition state. See the section transition state search.

In addtion the QuasiNewton class has the follwing methods:

SetHessian:
Set an intial Hessian matrix.
GetHessian:
Return the current Hessian matrix.
WriteHessian(filename):
Save hesian to file
ReadHessian(filename):
Read Hessian from file.

For each update the current Hessian matrix is output to file hessian.pickle. This file can be used, by using the ReadHessian command as an input to another similar calculation.

Having a good guess for the Hessian matrix can improve the rate of convergence. This is illustrated using the PairPotential:

from ASE.Calculators.PairPotential import PairPotential
from ASE.Dynamics.QuasiNewton import QuasiNewton
from ASE.IO.NetCDF import ReadNetCDF


for n in range(2):

    atoms = ReadNetCDF('slab.nc')
    atoms.SetCalculator(PairPotential())
    fmax = 0.01
    qn = QuasiNewton(atoms,fmax=fmax)
    if n==1: qn.ReadHessian('hessian1.pickle')
    qn.Converge()
    print "max force", max(abs(atoms.GetCartesianForces().flat))
    print "number of function evaluations", qn.GetNumberFuncEvals()
    qn.SaveHessian('hessian1.pickle')

This will produce the following output:

max force 0.00953150406961
number of function evaluations 17
max force 0.00435079765281
number of function evaluations 5

So in this case the number of function evaluation is reduce from 17 to 5. In theory, for a harmonic potential or for a potential in the harmonic region, only one step should be neccesary if the exact Hessian is known. Often one might be able to use the Hessian from a similar calculation, eg. if one is mapping out high symmetry positions for a molecule on a surface.

For the often used procedure of extending a bond-length, miminizing the remaning degress of freedom, along the way it can be a good idea to re-use the Hessian matrix from one bond-length to the next, this can significantly speed up convergence, see example.

Mimimizing the Nudged Elastic Band

A Nudged Elastic Band filter band can be minimized using the quasi-Newton algorithm:

from ASE.Dynamics.QuasiNewton import QuasiNewton
qn = QuasiNewton(band, fmax=0.05)
qn.Converge()

The quasi-Newton algorithm will mimimize the band as one listofatoms like object, so that the size of the Hessian matrix will be 3*N*M, N is the number of dynamics atoms and M is the number of images. This will limit the applicablity of the present inplementation to systems around 50 dynamics atoms and 10 images. More details of the method will follow shortly.

Warning

The quasi-Newton algorithm should not be used for the Climbing Nudged Elastic Band. The quasi-Newton algorithm minimize the band energy, and will not allow the climbing image to increase the band energy.

Tip

The tool nebfit can be used to check the progress of the algorithmm, and to make sure that the images are evenly distributed along the reaction path.

Molecular dynamics

There are currently three types of objects for molecular dynamics: Velocity Verlet, Langevin and NPT Dynamics. These classes are derived from the MolecularDynamics class (derived from the Dynamics class), and they have the Run() and Step() as well as a SetTimeStep() method. The time-step can also be set using the dt keyword when constructing the molecular dynamics object.

The following methods from the ListOfAtoms interface are used when doing molecular dynamics:

  • GetCartesianForces()
  • GetCartesianPositions()
  • SetCartesianPositions()
  • GetCartesianMomenta()
  • SetCartesianMomenta()
  • GetMasses()
  • (what about GetPotentialEnergy() ?????).

examples ... how to choose the time step, ...

Velocity Verlet

Velocity Verlet is the standard algorithm for molecular dynamics with constant number of particles, constant volume, and constant total energy (the NVE or microcanonical ensemble). It works by numerically integrating Newton's second law, using the well-tested velocity Verlet algorithm. This algorithm has been shown to have good properties for molecular dynamics, including good long-term stability.

The VelocityVerlet object is constructed with two arguments, the list of atoms and the timestep. Good values for the timestep depend on the potential, choosing it too small wastes computer time, choosing it too big causes the energy of the system to grow rapidly. A good value for copper with the Effective Medium Potential is 5 fs.

Langevin

Langevin dynamics is one possible way to do molecular dynamics with constant temperature (the NVT or canonical ensemble). The atoms are coupled to a heat bath with the desired temperature through. This is done through a modification to the equation of motion: a small friction and a fluctuating force are added to Newton's second law.

Note that the temperature is specified in energy units.

>>> from ASE.IO.NetCDF import ReadNetCDF
>>> from ASE.Calculators.PairPotential import PairPotential
>>> from ASE.Dynamics.Langevin import Langevin
>>> atoms = ReadNetCDF('copper_test.nc')
>>> atoms.SetCalculator(PairPotential())
>>> dyn = Langevin(atoms, dt=0.1, temperature=0.1, friction=0.01)
>>> for i in range(100):
...     print atoms.GetKineticEnergy() / len(atoms)
...     dyn.Run(50)
...
0.0
0.00850477560303
0.0166539513107
0.0198288079643
0.0299205976147
0.0307109339604
0.0363229195405
   [ ... ]
0.14992986571
0.1457703258
0.15382755386
0.14028708292
0.15933636908
>>>
NPT Dynamics

Dynamics with constant pressure (or optionally, constant stress) and constant temperature (NPT or N,stress,T ensemble). It uses the combination of Nosé-Hoover and Parrinello-Rahman dynamics proposed by Melchionna et al. [1] and later modified by Melchionna [2]. The differential equations are integrated using a centered difference method [3]. Details of the implementation are available in the document NPTdynamics.tex, distributed with the module.

The dynamics object is called with the following parameters:

atoms:
The list of atoms.
timestep:
The timestep in units matching eV, Å, u.
temperature:
The desired temperature in eV.
externalstress:
The external stress in eV/Å^3. Either a symmetric 3x3 tensor, a 6-vector representing the same, or a scalar representing the pressure. Note that the stress is positive in tension whereas the pressure is positive in compression: giving a scalar p is equivalent to giving the tensor (-p. -p, -p, 0, 0, 0).
ttime:
Characteristic timescale of the thermostat. Set to None to disable the thermostat.
pfactor:
A constant in the barostat differential equation. If a characteristic barostat timescale of ptime is desired, set pfactor to ptime^2 * B (where B is the Bulk Modulus). Set to None to disable the barostat. Typical metallic bulk moduli are of the order of 100 GPa or 0.6 eV/Å^3.
mask=None:
Optional argument. A tuple of three integers (0 or 1), indicating if the system can change size along the three Cartesian axes. Set to (1,1,1) or None to allow a fully flexible computational box. Set to (1,1,0) to disallow elongations along the z-axis etc.

Useful parameter values:

  • The same timestep can be used as in Verlet dynamics, i.e. 5 fs is fine for bulk copper.
  • The ttime and pfactor are quite critical[4], too small values may cause instabilites and/or wrong fluctuations in T / p. Too large values cause an oscillation which is slow to die. Good values for the characteristic times seem to be 25 fs for ttime, and 75 fs for ptime (used to calculate pfactor), at least for bulk copper with 15000-200000 atoms. But this is not well tested, it is IMPORTANT to monitor the temperature and stress/pressure fluctuations.

It has the following methods:

Run(n):
Perform n timesteps.
Initialize():
Estimates the dynamic variables for time=-1 to start the algorithm. This is automatically called before the first timestep.
SetStress():
Set the external stress. Use with care. It is preferable to set the right value when creating the object.
SetMask():
Change the mask. Use with care, as you may "freeze" a fluctuation in the strain rate.
GetGibbsFreeEnergy():
Gibbs free energy is supposed to be preserved by this dynamics. This is mainly intended as a diagnostic tool.

References:

[1] S. Melchionna, G. Ciccotti and B. L. Holian, Molecular Physics 78, p. 533 (1993).

[2] S. Melchionna, Physical Review E 61, p. 6165 (2000).

[3] B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover, Physical Review A 41, p. 4552 (1990).

[4] F. D. Di Tolla and M. Ronchetti, Physical Review E 48, p. 1726 (1993).

Adaptive Nudged Elastic Band Approach

The AdaptiveNudgedElasticBand dynamics implements the Adaptive Nudged Elastic Band Approch, see reference [1]. ANEBA is a variation of the Nudged Elastic Band method (NEB). In ANEBA a small number of images are used to bracket the transition state and then adaptively increase the number of images around this tranition state, read more in reference [1].

Setup the ANEB calculation like this:

from ASE.Dynamics.AdaptiveNudgedElasticBand import AdaptiveNudgedElasticBand
aneb = AdaptiveNudgedElasticBand(configs,prefix='al100',linear_interpolation=True)
aneb.Converge()

configs is a list of ListOfAtoms objects or a list of filters. Exactly 5 elements should be in configs. The option linear_interpolation is used to tell ANEB to make a initial linear interpolation between the first and the final configuration in configs.

Note that ANEB is implemented as a dynamics, using the NEB filter.

A restart script will look like this:

configs = [Calculator.ReadAtoms('al100_conf%d.nc'%n for n in range(5)]
aneb = AdaptiveNudgedElasticBand(configs,prefix='al100',level=2,subdivide=False)

level should be set to the level ANEB was working on before stopping. If subdivide=True ANEB will start up by subdividing the path and go to the next level. If subdivide=False ANEB will continue converging the current level, this is the default.

A complete list of keywords is listed below:

prefix:

Will set the name of the files. Restart files will be named:

out_<prefix>_conf<n>.nc if *keep_files=False*,
out_<prefix>_conf<n>_level<m>.nc if *keep_files=True*

Trajectory files will be called:

<prefix>_level<n>_conf<m>.nc`
keep_files:
If keep_files=False the same set of restart files will be used (default), otherwise new file names will be generated for each level.
level:
Enter the current level.
subdivide:
If subdivide=False ANEB will continue converging the current level, otherwise it will subdivide the path and move to the next level.
linear_interpolation:
If linear_interpolation=True a linear interpolation will be made between the first and the last configuration in configs.

ANEB use the quasi-Newton dynamics to minimize the band at each level. For each level the Hessian file name is set to be:

<prefix>_hessian_<level>

If this file exists then restarting a level it will be used. This file is deleted then ANEB proceed to a new level.

The tool anebfit can be used to assemble the path:

anebfit <prefix>

References:

[1] P. Maragakis, S.A.Andreev,Y. Brumer,D.R.Reichman, and E.Kaxiras, Journal of Chemical Physics,
117, page 4651 (2002)

ase2: Dynamics (last edited 2010-10-20 09:11:16 by localhost)