IMD
IMD Home
User Guide

Integrators
The integrator is the routine which actually moves the atoms,
depending on the current forces and velocities. The active integrator
is selected in the parameter file via the ensemble parameter.
Currently, each integrator to be used must be enabled at compile time.
There are integrators simulating different
thermodynamic ensembles, integrators used for
relaxation simulations, and some special
purpose integrators for the simulation of fracture, shear flow,
separate electronic temperature (TTM) or heat conduction.
Each integrator requires the use of a suitable
time step.
 nve

Constant volume, constant energy ensemble. This is the most reliable
integrator, and thus the main work horse. No particular parameters
are needed, except perhaps the initial temperature starttemp, if not
all atoms have their velocities given in the configuration file.
 nvt

Constant volume, constant temperature ensemble, with
NoséHoover thermostat (more precisely, this is an
ensemble with external temperature control  the temperature
can actually vary during the simulation). The imposed temperature
varies linearly from starttemp to endtemp (parameters).
A further required parameter is the
time constant of the temperature control, tau_eta.
Alternatively, the inverse of tau_eta, inv_tau_eta, can be specified.
For historical reasons, even the squared inverse isq_tau_eta is
accepted. inv_tau_eta or isq_tau_eta can be zero, which means
constant energy simulation (nve).
If the option UNIAX is switched on, an
additional parameter tau_eta_rot is required. It provides the
time constant for the coupling to the heat bath adjusting the
temperature for rotational motion. Like tau_eta, its inverse
inv_tau_eta, or its squared inverse isq_tau_eta, can be used
as parameter.
 and

This is actually an option which can be used together with the
nve or mik integrators. Every tempintv steps, the velocities
are reinitialized to a random maxwell distribution of the imposed
external temperature, which is linearly interpolated between
starttemp and endtemp. With the nve integrator,
this corresponds to an NVT ensemble with Anderson thermostat.
 ber

Use the Berendsen thermostat instead of NoseHoover.
 finnis

Local temperature control which resembles an electronic heat bath for the
ions according to Finns et al (PRB 44, p.567, 1999). It is basically a
simplyfied implementation of ftg.
 npt_iso

Constant pressure, constant temperature ensemble, with isotropic
volume scaling. NoséHoover type thermostat and volume control
are employed. Temperature control
is as in nvt, with the same parameters. (Again, if the option
UNIAX is switched on, an additional
parameter tau_eta_rot is required.) The pressure varies linearly
from pressure_start to pressure_end (parameters). A
time constant for pressure control,
tau_xi, or its inverse inv_tau_xi, or its squared inverse, isq_tau_xi,
is required as a parameter. If inv_tau_xi or isq_tau_xi is zero,
the volume remains constant. The cell subdivision includes a
rescaling tolerance cell_size_tol (parameter), which defaults to
0.05. If this tolerance is exceeded, the cell subdivision is
changed on the fly.
 npt_axial

Constant pressure, constant temperature ensemble,
with volume rescaling along each axis independently.
This works only for orthogonal boxes correctly.
Works the same way as npt_iso, with the same parameters, except that
the three axial pressures can be specified separately.
The parameter relax_dirs (which takes an integer, 0 or 1,
for each dimension) can be used to switch off pressure control in
those directions, for which the corresponding parameter component is
zero. This then corresponds to a mixed NVT/NPT ensemble.
 mik

Microconvergence integrator. After each step, if the velocity of an
atom goes "uphill" in the potential landscape, it is reset to zero.
No particular parameter is needed, except perhaps the initial
temperature, if not all atoms have their velocities given in the
configuration file.
 glok

This is a global version of the mik integrator: if the scalar product
of the global force and momentum vectors (containing the force and
momentum components of all atoms) is negative (momentum goes "uphill"
in the potential landscape), all momenta are reset to zero. Note that
this concerns only translational forces and momenta, not the
rotational ones in UNIAX. The momenta are
also reset to zero if the atoms are getting too fast, i.e., the
kinetic energy per atom exceeds the value of the parameter
glok_ekin_threshold. In some situations, glok relaxes
considerably faster than mik.
 adaptglok

The adaptive glok algorithm is an improvement of glok along the lines of fire.
The parameters are given in the parameter file.
 fire

The fire algorithm has been implemented by E. Bitzek, P. Koskinen,
F. Gähler, M. Moseler, and P. Gumbsch, see Phys. Rev. Lett. 97
(2006) 170201.
The parameters are given in the parameter file.
 cg

The conjugate gradient relaxator cg is controlled by the following
parameters. linmin_tol is the tolerance in the line
minimization. This is usually the only parameter that needs to
be set. The number of interations in line minimization is limited
by linmin_maxsteps (default 100. The step width
attempted in the first line minimization is given by
linmin_dmax (default 0.01). In the following
line minimizations, the optimal step width from the previous line
minimization is used. cg_infolevel = 1 enables some
diagnostic messages. cg_fr = 1 enables the FletcherReeves
variant of cg; by default, the PolakRibiere variant is used.
If cg_reset_int (default 0) is positive, the cg history
is reset every cg_reset_int steps.
All parameters are given in the parameter file.
 acg

Adaptive conjugate gradient
The parameters are given in the parameter file.
 relaxation control

Relaxation simulations with the relaxators mik, glok, and cg, are
controlled by the following common parameters. A sample is regarded
as sufficiently relaxed, if one of the following conditions
is satisfied: the kinetic energy per atom is below ekin_threshold,
the force norm is below
fnorm_threshold, or the absolute value of the potential
energy change since the last step is below delta_epot_threshold.
The computation of these quantities is automatically enabled.
The threshold values are zero by default, but one or several of them
can be set to a positive value in the parameter file.
If the sample is sufficiently relaxed, a snapshot of the configuration
is written (snapshot files, with an ending .ss, have the same format
as checkpoints, but a separate numbering). Moreover, if the
deform option is active
(max_deform_int > 0), the next deformation step is performed,
and the relaxation starts again. Further, if
force boundary conditions are active, and the force increment
extra_dforce is nonzero, it is added to the boundary forces,
and the relaxation starts again. If neither the deform option or force
boundaries are active, the current simulation phase is terminated.
Note that for a relaxation simulation with mik or glok, the timestep
can often be chosen much larger than for a proper molecular dynamics
simulation. A factor of 10 larger is often non unreasonable. A little
experimentation may be needed to find the largest timestep, for which
the simulation still converges. Since time has little meaning in a
relaxation context, the number of force computations, nfc, is written
to the .eng file instead of the time, if some relaxator is compiled
into IMD (this is so even for an annealing phase with nve or nvt
before the relaxation). For mik or glok relaxators, nfc agrees with
the step number, but for the cg integrator a single step may comprise
a large number of force computations. Since the force computation
takes the main part of the simulation time, nfc allows for a fair
comparison of the speed of convergence between different integrators.
 pressure relaxation

In combination with relaxation integrators, pressure relaxation can
be switched on by setting the parameter relax_rate to a
positive (sufficiently small) value. At each step, the pressure tensor
of the sample is then calculated, and the box size and shape is
changed by a small amount in order to lower that pressure. The
parameters bulk_module and shear_module, given in
units of pressure, can be used to specify
the sensitivity of the pressure tensor to box changes. With correctly
set bulk and shear modules, a value of 0.1  0.5 for the
relax_rate seems to work well. For higher values, the
system might become unstable due to pressure fluctuations, which
can be quite large initially. If large volume changes are needed
to relax the pressure, or if the estimate for the elastic constants
is not accurate enough, a smaller value may be necessary. If the
bulk and shear modules are not given, they default to 1. Some
experimentation with relax_rate should then quickly reveal
a usable value.
The parameter relax_mode, which can take
values full, axial, and iso, determines wich
part of the pressure tensor is relaxed. full means full
relaxation, axial means relaxation of pressure along
main axes of a rectangular box, and iso relaxes the scalar
pressure by isotropic rescaling of the box. All pressure relaxation
modes require the compile option homdef,
and modes full and axial also the compile option
stress. Default mode is
full if stress
is compiled, and iso otherwise.
If option stress is compiled,
it is also possible to relax to a given external pressure tensor read
from the parameter file, in the following form:
presstens_ext xx yy xy # 2D
presstens_ext xx yy zz yz zx xy # 3D
The default for presstens_ext is zero. All three relaxation
modes are supported. If relax_mode = axial, only the diagonal
pressure components are relaxed, and if relax_mode = iso,
only the isotropic pressure is relaxed to the average of the
diagonal components. In any case, the full presstens_ext
must be given, even if some of its components are ignored.
If relax_mode is axial, the parameter relax_dirs
(which takes an integer, 0 or 1, for each dimension) can be used to switch
off relaxation in those directions, for which the corresponding parameter
component is zero.
 frac

Fracture simulation with stadium damping according to Gumbsch et al.
(PRB 55, p.3445, 1997). The crack (x: propagation direction, y: crack
plane, z: 2D or pbc) can propagate within the nve ensemble in an
elliptical "stadium" (characterized by center and stadium (= half axes)).
Outside this inner ellipse viscous damping or nvt temperature control is
applied (depending on dampingmode). The damping is characterized by the
damping coefficient gamma_bar which should be around 2*omega_E.
The ramping of the damping depends on the ratio stadium.x/box.x, i.e full
damping at the border of the box. With stadium2 one can introduce an
outer ellipse outside which full damping is applied, the slope of the
damping thus depends on the ratio stadium.x/stadium2.x. Additionally the
sample can be strained at strainrate in y direction.
 sllod

SLLOD ensemble for simulation of shear flow
with the help of a deformation rate tensor in combination with
LeesEdwards boundary conditions. It corresponds to an additional
velocity component applied to every particle in every step. It appears
also in the equation of the momenta so that the particles perform a
motion relatively to the shear flow. The SLLOD algorithm
must be used together with a thermostat and with
periodic boundary conditions in all directions. In 2D a vector
shear_rate must be specified whose two components correspond to the
xy resp. to the yxcomponent of the applied deformation rate
tensor. In 3D two vectors with three components, shear_rate and
shear_rate2 must be specified. The components of shear_rate specify
its yz, the zx, and its xycomponent. The components of shear_rate2
specify its zy, the xz, and its yxcomponent.
 ttm

Basic nve integrator that has been enhanced for metallic samples far from
thermal equilibrium, as they occur for example in simulations of laser ablation
with ultrashort pulses. This has been achieved by implementing a hybrid method
for the Two Temperature Model (TTM), which assigns separate temperatures to the
subsystem of the conduction electrons (T_e) and to the lattice (T_p).
The differential equation for the electronic temperature is solved in a
continuum ansatz on a discretised lattice with a finite difference (FD)
method. Every node/cell of this FDlattice corresponds to a cuboid in
simulation space comprised of a number of MD cells. The local lattice
temperature is derived by averaging over the atoms in said space, and the
energetic coupling between the two systems is made possible by a damping term
in the integrator proportional to the local temperature difference (see
D. S. Ivanov, L. V. Zhigilei, Phys. Rev. B 68, 064114 (2003) for
details).
This integrator is still a little awkward to use, and subject to limitations:
All materials parameters are constants (except for the electronic heat capacity
c_e, which can be made proportional to the electronic temperature by use of the
parameter fd_gamma). Aborted simulations cannot be restarted from snapshots but
only from the beginning.
Necessary Parameters: fd_k, fd_c or fd_gamma, fd_g, fd_n_timesteps, fd_ext.
Optional Parameters: fd_update_steps, ttm_int, init_t_el, fix_t_el.
The subdivision of the simulation volume into FD cells has to be supplied by
the user and depends on the dimensions of the automatically generated MD cell
arrays. This must be done with the parameter fd_ext, which takes a vector of
integers that give the number of MD cells per FD cell in the coordinate
directions of the simulation space. If one of these values is selected too
large, so the resulting MD cells would be smaller than the minimal cell size
(cutoff radius of the potential), IMD will issue an error and exit. Otherwise,
the MD cell size is automatically set to a value that allows the simulation
space of a single process to be completely filled by tiling FD cells of size
fd_ext, i.e. global_cell_dim.x%(cpu_dim.x*fd_ext.x)==0.
To obtain the MD cell array dimensions, one short test run with the input
checkpoint file and the desired cpu_dim can be done. Example: A test run with
two processes (without ttm code) yields the following output:
Reading atoms.
Minimal cell size:
( 5.960000 0.000000 0.000000 )
( 0.000000 5.960000 0.000000 )
( 0.000000 0.000000 5.960000 )
Actual cell size:
( 5.969499 0.000000 0.000000 )
( 0.000000 6.752167 0.000000 )
( 0.000000 0.000000 6.752167 )
Global cell array dimensions: 1018 6 6
Local cell array dimensions (incl buffer): 511 8 8
MPI process array dimensions: 2 1 1
Read structure with 600000 atoms.
...
maximal cell occupancy: 27
Done reading atoms.
Here, the sample has a size of 5.969499*1018 x 6.752167*6 x 6.752167*6 or
6076.95 x 40.513 x 40.513 Angstroms. The minimal cell size is 5.96 Angstroms
(the cutoff radius of the used potential). Because there are two processes,
which are aligned in x direction, the maximum values for fd_ext are 509 6 6
(which happens to be the local cell array dimensions, without buffer
cells). Note that choosing large (albeit allowed) values without considering
the minimal cell size may result in MD cells that are larger than necessary,
which may slow down the force computation loops, as they scale quadratically
with the number of atoms per MD cell.
The ttm code adds a column to the .eng file with the average electronic energy
per atom E_el. Note that this is _not_ the electronic temperature and will
likely differ from the average lattice temperature in the same file even for a
relaxed sample. In order to look at electron temperatures, assign a nonzero
value to the parameter ttm_int and use the information from the *.ttm
files. When compiled with the debug option, there will be some more columns in
the .eng file that will be useless if the code runs correctly.
For reasons of numerical stability, the timestep of the MD part and of the FD
part of the simulation differ. In general, the FD timestep is smaller by a
factor fd_n_timesteps, which must be given as a parameter. It has to fulfill
the stability criterion
n>2*timestep*fd_k/(h*min(c_e)),
with h the minimum thickness of the FD cells (1*5.969499 in the example,
because it is smaller than 6*6.752167) and min(c_e) the minimum electronic heat
capacity to expect in the whole course of the simulation (i.e. fd_c, or
alternatively fd_gamma multiplied by the lowest expected electron
temperature).
If you choose a wrong (i.e., too small) value for fd_n_timesteps, bad things
may happen, including periodic temperature fluctuations, growing exponentially
in time, and negative electron temperatures.
 nvx

Used for heat conductivity measurements. Forces a temperature
gradient into the sample.
 ftg

This ensemble can be used for crack simulations at finite temperatures, where a
crack moves from a cold part of the sample to the hot part and stops
there. Implemented by Cristoph Rudhart.
 stm

Who knows
