Comparison of time- and length-scales, accessible to different types of simulation techniques (quantum simulations (QM), molecular dynamics (MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Computing every force on each atom is a demanding part of a MD simulation. In order to limit the number of force calculations, the smallest possible system – in terms of number of atoms – should be simulated. The minimum system size necessary to compute a given desired property is not always obvious and a common practice is to analyze how increasing system size influences obtained results. A number of algorithmic approaches exist that help minimize system size such as applying periodic boundary conditions.
In addition to how large a system must be, one must determine how long a simulation must be run to obtain reliable answers (i.e. with sufficient accuracy and precision). Related to this is that the MD simulation time step should be small compared to the period associated with the highest frequency atomic event. A typical time step for MD simulations of metals is of order 1 fs (i.e. 1 E-15 s). The very small time step is the main drawback of MD method compared to other commonly used numerical simulations. Processes with time scale of order microseconds and larger are outside the reasonable realm for even highly parallel, efficient MD codes. Conversely, MC techniques cannot probe dynamics but they can probe long time equilibrium properties.

Figure 2.1 Schematic comparison of time- and length-scales, accessible to different types of simulation techniques (quantum simulations (QM), molecular dynamics (MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Accessible length- and time- scales of microscopic simulation should be considered. Figure 2.1 shows a schematic representation for different types of simulations in a length-time-diagram. It is clear that the more detailed a simulation technique addresses, the smaller is the accessibility of times and length scales. Hence, quantum simulations (QS), where fast motion of electrons are considered, are located in the lower left corner of the diagram and typical length time scales are of order of Å (Å) and ps. MD method puts fixed partial charges on interaction sites or adds a dynamics model for polarization effects to simulate electronic distribution.  Hence, the time scale of the system is determined by the time of collision between atoms and atomic vibration frequencies. Consequently, the accessible time scale is of order ns and corresponding accessible length scale becomes of order 10-10000 Å. Note Figure 2.1 was originally presented in 2002; accessible length scales for MD have increased dramatically. Brownian Dynamics (BD) can be applied to deal with particles in a solvent medium, where one is not interested in a detailed description of the solvent. In BD the effect of the solvent is hidden in average quantities. If one is not interested in a molecular picture but in macroscopic quantities, the concepts of hydrodynamics (HD) may be applied, where the system properties are hidden in effective numbers, like density and sound velocity.

introduction to Molecular dynamics, part 2: The MD Algorithm; why ∆t equal to 1 femtosecond?

2.2 The MD Algorithm

2.2.1 Verlet’s Algorithm and Velocity Verlet’s Algorithm

The governing equations of motions are all ordinary differential equations and a
standard method for the solutions is the finite difference method. Of a large number of
finite-difference approaches devised, the most commonly used one is the Runge-Kutta
(RK) method. [33,34] The beauty of the RK method is that the order of a finite
difference method can often be increased by using positions and velocities from
several points in time instead of only the current time. However, if a fourth-order RK
method is applied to MD simulation, four interatomic force evaluations per atom
would be required per step, which will make the simulation too slow. Thus, reduction
of order is desired; especially, the time expense of calculating intermolecular force
more than once per atom per step should be also avoided. However, the algorithm
must be stable.
The most widely used method of integrating the equations of motions in molecular
dynamics is Verlet‘s method or algorithm. It uses positions and velocities from
previously evaluated steps to provide a direct solution of equation (2.4). In this study,
the simulation code LAMMPS is used and its default algorithm is the velocity Verlet.
[35] Hence, the basic idea of Verlet‘s algorithm, and specifically the velocity Verlet,
will be presented here. Information on other integration schemes can be found in texts
on atomic scale simulation. [32]
This standard, or basic, Verlet algorithm is a combination of two Taylor‘s
expansions as following. The Taylor series for position from time t forward to 𝑡 + ∆𝑡

 And the series form t backward to 𝑡 − ∆𝑡 is

 Adding equations (2.6) and (2.7) with some simplifying produces

 (2.8)

This provides Verlet‘s algorithm for positions. The truncation error of the algorithm
when evolving the system by ∆t is of the order of ∆t^4. The acceleration needed in equation (2.8) is obtained from the intermolecular forces and Newton‘s second law. Velocities are not explicitly solved in equation (2.8), and they are calculated from first-order central difference.

A more commonly used, related algorithm is the velocity Verlet algorithm. The
approach is similar with the standard one but starts with position and velocity
expansions explicitly

The implement scheme of this algorithm is following:
1. Calculate velocities at mid-step,
Calculate the positions at next step 𝑡 + ∆t
2. Calculate the accelerations at 𝑡 + ∆𝑡, 𝑎(𝑡 + ∆𝑡) from the interaction potential using 𝑥(𝑡 + ∆𝑡).
3. Calculate the velocities at 𝑡 + ∆t

The advantage of this algorithm compared to the standard one is that we never need to simultaneously store the values at two different times for any of these
quantities. The velocity Verlet algorithm gives good feature of time reversibility,
simplicity and stability for even larger time steps compared to standard algorithm.



Ref: http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University
2011


---------------------------------------------------
Verlet is an efficient algorithm, and its benefits are stability, time-reversibility and
energy conservation.
The time step, ∆t, should be chosen in such manner that the simulation system
is stable and energy is conserved, that is to say that the provided ∆t should be
short enough. Usually, in classical MD simulations time steps of order 1–2 fs are
being used (especially when using constraints), but in some cases time steps can
be even 4–5 fs (when using virtual site algorithms). The length of the time
step is dominated by the fastest degrees of motion in the system, and usually theses
are bond vibrations. Most often constraints are used in simulations to reduce the
degrees of freedom, removing high-frequency vibrations by restricting the movement
of covalent bonds. This is implemented by fixing the covalent bonds to constant
angles and lengths. In addition to constraints, also restraints can be used to fix

atomic positions, bond lengths, bond angles and dihedrals.