IMD

IMD Home

User Guide

Compiling
Usage
General Remarks
Units
Parameter File
Valid Parameters
Configurations
Integrators
Ensembles
Relaxators
Special
Interactions
Potential Files
Pair Potentials
EAM
ADP
MEAM
Covalent
Coulomb
Gay-Berne
OpenKIM
Options
Output
Utilities
Distributions
Configurations
Correlation
Potentials
Miscellaneous
Online Visualization
Limitations
Developer's Corner

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.

Thermodynamic Ensembles

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.

langevin

Another option for nve and mik that enables the Langevin thermostat. Every timestep the velocities are reduced by viscous damping and random forces are added to adjust the temperature. The required friction coefficient is the same as in the option viscous.

ber

Use the Berendsen thermostat instead of Nose-Hoover.

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.

Relaxators

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 Fletcher-Reeves variant of cg; by default, the Polak-Ribiere 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 non-zero, 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.

Special Purpose Integrators

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 Lees-Edwards 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 yx-component 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 xy-component. The components of shear_rate2 specify its zy-, the xz-, and its yx-component.

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 FD-lattice 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 non-zero 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