Solving the Equation of Motion

In principle, MD simulations describe the motions of particles (ions or atoms) by solving the equation of the motion. That is, the position of the particle is predicted from the previous and present ones by adding the forces acting from other particles. In the case of classical MD, it would not be an exaggeration to say that potential functions and its parameters (force fields) determine the fate of the particles, if other conditions are reasonably selected.

The equation of the motion for the /-th particle can be written as

Algorithm

In classical MD simulations, the equation of the motion is solved numerically. In other words, the next position of the /-th particle is calculated based on the position at t and that of one or several steps before. Several algorithms to solve the equation of motions have so far been proposed. Verlet algorithm [7] and Gear's method [8], which is one of the predictor-corrector methods with several steps, are also used for MD simulation [1-3].

The former is suitable for the calculation of motion including jump or hops of ions often encountered in the simulations of melts and glasses because it does not drag the memory of past motions for a long time.

Verlet method is explained as follows [7]:

where the vector r(t + At) is for positions of the /-th particle after At and r(t- At) is that before At.

One can obtain the following relation using the sum of Eqs. (4.2) and (4.3):

while a following equation was obtained, from the difference of Eqs. (4.2) and (4.3):

That is, new position of the particle i is represented by and a new velocity is represented by

A numerical error in solving the equation of motion thus depends on the time of each step, At, and hence too large step should not be used. However, if the particle is localized for a long time, often a longer time step can be used because the structures of cages and or paths can have longer relaxation times. In many of ionic systems, each time step At can be 1-4 fs. For the case of the system including light particles such as hydrogen, the smaller time step is necessary (typically, that less than 1 fs for the case of classical MD and even a shorter time is necessary for ab initio MD).

Required Conditions for MD Simulations

In conventional MD, periodic boundary conditions (PBC) are used to mimic the situation of bulk When one uses PBC, caution is required for the size of the basic cell to be used. It may modify both the motion and the structure by repeating those in the basic cells. It also determines the wavenumber, which can be addressed.

Needless to say, good quality of potential parameters suitable for the purpose of the work is also required.

When one handles nanoscale structures with different treatments of boundary conditions, additional caution will be required. For example, if the surface of the materials is treated, often the region used for summation to take the statistical treatment needs to be changed and additional normalization procedures should be used for comparison with other systems.

In this chapter, some of these problems are discussed one by

one.

Periodic Boundary Conditions and Finite Size Effects

To treat the properties of bulk, periodic boundary condition (PBC) is often used to remove the effects of the surface in MD simulations [2-4].

That is, the basic cell containing particles is surrounded by the periodically repeating image cells.

For a particle within the basic cell, interactions from particles within the sphere with a certain cutoff length (typically chosen to be L/2, where L is a length of a basic cell) are considered, including ghost particles (images of the same particle found in the basic cell) in surrounding image cells. By this treatment, properties of the bulk can be simulated by the limited number of particles. If the particle moved beyond the boundary, the ghost particle comes into the basic cell. Therefore, in the case of constant number simulations such as NVE (constant number of particles, volume, and energy) or NPT (constant number of particles, pressure, and temperature) ensembles, the number of the particles is kept constant by the PBC.

Usually, Coulombic force can be calculated by the Ewald method [9] or by other methods such as particle mesh Ewald [10]. Even in the case of liquid or glasses can be treated by using the PBC, which mimics the periodicity like crystals. In many cases, the minimum image (within 2/L) is used for the calculation of short-range forces. The same length is often used for the real part of the Ewald summation for the calculations of Coulombic forces.

In liquids or glasses, the structure shown by pair correlation function g(r), decays within a certain length. For ion-ion interactions, the correlation of the structure (deviation from the line of g(r) = 1) is not clear at distance longer than 8-12 A in ionic systems such as lithium silicate glass [3]. Therefore, about twice this length is a minimum length of the unit cell for MD simulations to avoid the ionic structure being affected by the PBC.

Sometimes, system size used in molecular dynamics simulations is not large enough and it causes several kinds of finite size effects in the results of simulations. This problem limits the usage of ab initio MD, which requires more calculation resources than classical MD. Even in the classical MD, it is not easy to treat the large system size required for the treatment of many systems with different purposes. Here we consider several origins of finite size effects to be considered.

• 1. Even when we considered the infinitive system using PBC, the wavenumber accessible by the simulation is limited by the size of the simulation box. Furthermore, when one considered the motion of particles in the system with PBC, it is repeated as well and some artificial waves or vibrations in the particle motions will be formed.
• 2. Especially in the case of crystals, a basic box of the MD is formed by several repeating basic lattices of the crystal and therefore the number of repetitions of them in each axis direction will affect the periodicity of the motions. Consistently, mobility of ions is found to be affected by the system size considerably [11]. To reduce such effects, the system size used should be large enough to the possible extent, while ensuring the practical usability.
• 3. In the case of melt and glasses, cooperative motions of network formers and ions also may be affected by PBC. If the system size is too small, the particle might be affected by its own ghost in an image cell, which is moving in the same direction.

It is useful to change the box size to check the finite system size effect mentioned above. Interestingly, the size effect in nanosystems is also found in the study of ferroelectricity, characterized by a spontaneous electric polarization. Shimada and collaborators [12] have examined the ferroelectricity in edged PbTi0₃ nanowires under axial tension by ab initio density-functional theory calculations. They have reported that "surprisingly, the smallest PbO-terminated nanowire with a cross section of only one-unit cell can possess ferroelectricity while ferroelectricity disappears in the Ti0₂-terminated nanowires with a cross section smaller than four-by-four cells diameter of about 17 A. However, ferroelectricity is recovered by axial tension, where smaller nanowires require larger critical strains.” These results are interesting to consider both the size effects and roles of the epitaxial strain.

System Size Effects Related to Crystallization

The remarkable size effect is also found in the calculation of thermal conductivity [13]. For the study of glasses, further caution is required to avoid the crystallization of the system. If only a small number of particles are contained in each basic box of the simulation, the system may easily crystallize and/or behaves like crystals because of PBC. When systems have network structures, long-ranged oscillation tends to be formed and continued by PBC. The motions of particles to be examined will be modified by these effects.

Required System Sizes and Related Conditions

Here we will discuss how large system size is required for meaningful arguments of structures and dynamics obtained by MD simulations. Such requirements are common for both classical and ab initio MD simulations. A part of arguments is found in [4] and is expanded for complex systems treated in the present book here. Of course, it will depend on the purpose of simulations, but here we consider some typical scales for the examination of structures and motions of particles such as the diffusion.

In general, one needs larger systems to access the structures and dynamics related to the smaller wavenumber. If the system size is too small, the particle might be affected by its own ghost in an image cell, which is moving in the same direction. In the case of slow dynamics near the glass transition temperatures or in viscose liquids, the size of the cooperatively rearranging region (CRR) of the systems [14, 15] should also be considered. With the changes in the size of CRR with temperatures and/or pressure of the system, the required system size will also change. Typical systems treated in recent MD works of ionic systems seem to be containing several thousand particles. However, a larger system size will be required if one examines the more complex systems with hierarchy structures and/or different domain sizes. The reason for the requirement of larger size can be easily found in the example of the structure of gels (see Chapters 2, 8 and 9.)

Length scales in ionic systems such as Debye length

More specifically, we will consider some typical systems.

In simple ionic systems such as molten NaCl, long correlation length, including the Debye length, is found by the existence of charge density wave as discussed later in Section 4.8. Even for the crystals, larger systems will be required if one considers the fluctuation of the structures or long-ranged motions. Therefore, system size must be determined to represent these length scales reasonably.

Then we will consider the example of silicate glasses. In the strong system with three-dimensional networks, a large finite- size effect is known, while the effect is relatively small in fragile liquids [16]. If attention is paid to the local structures in a fragile system with larger alkali contents, relatively small systems can be used, while if one would like to discuss the network statistics such as the distribution of rings, larger system size will be necessary. Of course, it is better to check several length scales in the systems. As discussed in Chapter 2, these length scales can be modified by nano-sizing. They will be modified by other perturbations such as the mixing of components and formation of composites.

Next, we will consider the example of simple ionic liquid (ЕМ1М⁺-ЫОз), (here EMIM⁺ means 1-methyl 3-ethyl imidazolium ion). In this system, correlation in g(r) is observed until ~18 A at 400 К and therefore the system size with ~36 A of L will be a good choice for many purposes, while larger cell size will be required if the structure has long tails (chains) and/or if one is considering the existence of large-scale domains.

Besides the system size, other conditions such as the shape of the basic cell may need to set in the MD simulations, especially for complex systems. When some alignments of molecules or ions are found, a non-cubic basic box of the simulations might be better to use.

In the study of nanostructured systems such as a nanoparticle and/or each isolated system, the system without PBC might be used. In that case, a direct sum or by multipolar expansion can be used for calculations of the Coulombic Forces. If one uses the system without PBC, additional care is necessary for the calculation of the structural properties such as pair correlation functions g(r). Usually, the function is calculated using forces within a radial distance r and normalized by the number of the particles within the sphere. It needs to be modified by the structures of boundary, dimensionality, and shapes in the case of nanostructured systems.

As shown in Chapters 8 and 9 for the study of aggregates and gels, the contribution of large length scales is found. For example, to treat nanocolloidal-silica-water-NaCl systems in Chapters 8 and 9, the system containing tens of thousands of particles was used. Even a larger system is required in some objectives of works. A similar situation is also observed in other systems, such as network-forming systems and/or composites.

With an increase in systems size, the time required for equilibration of the system also becomes longer. For the simulation of equilibrated properties of large systems, one should be careful to equilibrate systems sufficiently, although, for some special purposes related to the problem of the glass transition, rapidly quenched systems are also examined. The required equilibration time is much affected by the dynamics of the particles in the system. In the system with slower dynamics, longer equilibration time is required.

Time Scales to Be Covered

In non-equilibrium systems such as supercooled liquids and glasses, non-equilibrium relaxation, that is aging, overlaps. In such a case, temperature increases gradually during the runs in the NVE condition and this fact is useful to distinguish the status of the system. Due to this reason, we need to distinguish the term such as "NVE condition” from "NVE ensemble” in some works [17]. For example, in the case of the porous system formed in NVE conditions, relaxation of the system to form the larger voids is found as shown in Chapter 6. This formation of larger voids caused by the stress relaxation is closely related to the mechanism of the fracture. If one examines the mechanical response of the system, the role of stress relaxation time and its relation to the conditions of MD should be considered. It should be noted that the fracture in NPT conditions in MD simulations is not necessarily directly comparable to the experiments under atmospheric conditions. Crystallization also overlaps the equilibration process.

In the case of ab initio MD simulations, available equilibration time may be restricted, and hence, using the equilibrated structure obtained by classical MD will be a good choice for the initial configurations, if available.

Due to the complexity of structures, dynamics of ions as well as atoms, molecules in nanostructured materials may contain several different time scales of dynamics. Even simple glassforming materials and/or simple ionic system, it shows several length-scale and time-scale regions in the mean squared displacement (MSD), which corresponds to the caging to diffusive regimes through the power law region (see Chapter 5). Sometimes, a time region of local motions seems to be erroneously regarded as the diffusive regime, and hence, a careful judgment of the regions to be examined is required.

Relationship between Surfaces of Nanostructured Materials and PeriodicBoundary Conditions in MD Simulations

In the field of molecular dynamics (MD) simulations, PBC is usually used to reduce the effects of the surfaces when one treats the bulk systems as already mentioned.

However, even though such conditions to mimic the infinite structure are applied, it is known that properties of systems strongly depend on the size of a basic box used in the simulations and the number of particles included there. Such finite size effects are caused by restriction by introduced repeating structures and motions imposed by PBC as already discussed. That is, the length scale, which can be treated in the simulation, is limited by the system size. This fact is interesting to consider the situation of nanostructured materials.

The situation of some of the nanostructured materials such as nanoparticles is comparable to the MD cell for which PBC is removed. Some systems may not be stable under such conditions. Obviously, in this case, long-ranged structure and motions are terminated at the surface. In other words, nanoparticles or films, materials trapped in pores, properties are affected by the absence of the short and long length-scaled structures and cooperative motions at the surface. Of course, rearrangement of atoms occurred on the surface.

This situation will be enhanced in further nano-sizing.

It is possible to simulate a surface of nanoparticles without PBC. Actually, simulations of small clusters have been done without PBC in several works. For example, in the simulation of several sizes of the cluster of silica [18], it was found that smaller clusters have a higher density, but the values are in close agreement with that calculated for bulk silica in the coesite phase. If PBC is removed to mimic the surface of the nanoparticles, one should not forget that the long scale forces and dynamics are also cut by it. That is, properties of nano-sized particles are not necessarily caused by the effects of the large surface area alone but the inclusion of the changes in structures and their length scales.

As shown in Chapter 8, nanocolloidal silica was prepared as the cluster of Si0₄ units without PBC [19] in our MD simulations. Then such model colloids are used as a part of the colloidal solutions.