Basic simulations of molecular dynamics - Molecular Dynamics Simulation

Chapter 2 Molecular Dynamics Simulation

2.1 Basic simulations of molecular dynamics

Molecular dynamic (MD) has generally used in simulating the structure of liquids, solids, and droplets. We will introduce briefly the molecular dynamic simulation and how to define the model applied to simulate the behavior of a droplet impacting a solid substrate or another droplet.

First of all in a simulation, we should consider the location of all atoms. By the location and the time derivative, all atoms can be defined their point masses in molecular dynamic simulation. Each atom as a particle interacts with the other particle through interaction forces derived from interaction potentials in the system, and time evolution is governed by Newtonian mechanics. At each time step, the acceleration of a particle used to update the position of the particle would be

calculated by Newton’s second law of motion as follows (2.1).

Since Newton’s second law is time independent or equivalently

is invariant under time translations. Consequently, we expect there to be some function of the positions and velocities whose value is constant in time; this function

is called the Hamiltonian, H, ( ^N, ^N)

H r p =const (2.2)

Here, the momentum of molecule is defined in terms of its velocity by . For an isolated system, total energy E is conserved, where E is equal to

the sum of kinetic energy and potential energy. Thus, for an isolated system, we identify total energy as the Hamiltonian; then for N spherical molecules, H can be written as

where the potential energy results from the intermolecular interactions. First consider the total time derivative of the general Hamiltonian (2.2),

( ^N)

H has no explicit time dependence, then the last term on the RHS If, as is (2.3),

of (2.4) vanishes and we are left with

i 0

Now consider the total time derivative of the isolated-system Hamiltonian given in (2.3),

On comparing (2.5) and (2.6), we find for each molecule i.

Substituting (2.7) into (2.5) gives

Since the velocities are all independent of one another, (2.10) can be satisfied only, for each molecule i, we have

Equations (2.7) and (2.11) are Hamilton’s equation of motion. For a system of N particles, (2.7) and (2.11) represent 6N first-order differential equations that are equivalent to Newton’s 3N second-order differential equations (2.1) .

In the Newtonian view, motion is a response to an applied force. However, in the Hamiltonian view, motion occurs in such a way as to preserve the Hamiltonian function, where the force does not appear explicitly.

For an isolated system, the particles move in accordance with Newton’s second law, tracing our trajectories that can be represented by time-dependent position vectors ri(t).

Similarly, we also have time-dependent momentum pi(t). At one instant, there are positions and moment of the N particles in a 6N-dimensional hyperspace. Such a space, called phase space, is composed of two parts: a 3N-dimensional configuration space, in which the coordinates are the components of position vectors ri(t), and a

the components of the momentum vectors pi(t). As time evolves, the points defined by positions and momentum in the 6N phase space moves, describing a trajectory in phase space.

2.2 Equations of motion

Potential energy represents the functions of all atoms in the system. Because of the complexity of these functions, there is no analytical solution to the equations of motion, therefore we must use numerical theorems in the operation. There are many numerical algorithms has developed to calculate the integrating equations of motion.

One of them is Leapfrog method [Frenkel, et al., 1996]. Leapfrog method is the most general method applied in MD simulation and there are some reasons for it. When we put the leapfrog method in the integrating equations of motion, the values of calculation will be accurate to third order. From the viewpoint of energy conservation, it tends to be considerably better than higher-order methods. Furthermore, its requirement of the memory storage of the computer calculation is lower than other methods.

The introductions of Leapfrog method are as follows : The Leap-Frog method use the velocity at

t−dt and to compile the interaction

force of atoms. Use the calculated force to compute the velocity at t

t+dt after a

time step, as (2.12) equation. And then apply the velocity at 2

t+dt to get the

location of the atom at 2

t+dt , as (2.13). The process of operation above is relatively

simple because the method we use has no need of modification.

* ( )

To diminish the arithmetic error, we adopt the velocity at 2

t+dt to solve the

location at next time step rather than the velocity at . In other words, the computed solutions of the time average velocity between the velocity at and t would have better accuracy than which only applied the velocity at . The calculating processes as Fig. 2.1.

It is very important to choose an appropriate potential energy model in MD simulation when we want to approach some different realistic materials. The force derived from the potential on an atom is due to the interaction with surroundings.

Lennard-Jones potential, one of potential methods derived from complicated computation, will be used in the two droplets upon a solid wall model. Why we choose Lennard-Jones potential as our main method will be explained and discussed in the next section.

2.3.1 Lennard-Jones potential

Lennard-Jones potential is proposed by J.E. [Lennard Jones, 1924]. The potential energy of a pair of atoms, and i j, locates at r_i and r_j as follows eq.(2.14).

The functional forms of the LJ and soft-sphere potentials in MD units are shown in Fig. 2.2 (a). Where r_ij = −r_i r_j, and the parameter ε govern the strength of the

interaction and σ defines a length scale. Fig. 2.2 (b) is a relative curve of the xenon and helium in argon Fig. 2.2 (a). The interaction repels at close range, then attract, and is eventually cut off at some limiting separation . While the strongly

repulsive core arising from the nonbonded overlap between the electron clouds has a rather arbitrary form, the attractive tail actually represents the van der Waals interaction. For the interactions involve individual pairs of atoms, each pair is treated independently, with other atoms in the neighborhood having no effect on the force between them.

to zero, so that there is no discontinuity at r_c ; f and higher derivatives are

discontinuous, thought this has no real impact on the numerical solution.

In addiction, these units of equation are usually expressed in dimensionless format in MD simulation. The advantage of using dimensionless format is that the equation of motion becomes simpler because some parameters defined in the model will be absorbed by dimensionless units. Another advantage from the viewpoint of a particle is that transferring these units into dimensionless format can avoid some calculating errors in computer hardware while the atoms are over-ranged.

2.4 Force computations 2.4.1 All pairs

All pairs is the simplest method to implement, but extremely inefficient when the interaction range is relatively small compared to linear size of simulation region.

In the simulation domain, all pairs of atoms must be examined in computing process as Fig. 2.3 (a). In fact that the amount of computation needed grows as rules out the method for all but the smallest value of . is the number of particles.

There are two techniques which can reduce the growth rate and they will be discussed in next two section.

Cell-link which can avoid most of unnecessary compute and improve the computational effects provides a means of organizing the information about atom positions into a form. This method is to divide the simulation region into a lattice of small cells, and the cells should exceed in width. Then if particles are assigned to

cells on the basis of their current positions, interactions are only possible between particles which are either in the same cells or in adjacent cells. Because of symmetry only half the neighboring cells need to be considered. For example, a total of 14 neighboring cells must be examined in three dimensions (include the cell itself).

The wraparound effects are readily incorporated into the scheme. Distinctly the region size must be at least 4 longfor this method to be useful. There are several

ways for implementing this cell-link list method to connect the relation between particles and cells. In the current demonstration code, it utilizes the concept of the pointers for particles and cells. Each cell stores a particle number, which may be zero or nonzero. Nonzero value represents a true particle number, while the zero value represents either the last atom in the cell or an empty cell. In addition, only one array cell List is used to represent the particles and the cells. An obvious advantage of doing this is that we could know exactly the size of this array in advance if periodic boundary conditions are used. Of course, there are several other methods to implement this idea of cell-link list technique. Ideas depicted above can be clearly

illustrated by Fig. 2.3 (b).

2.4.3 Neighbor Lists

During MD simulation, in order to consider the distance between each atom is all in the range or not, we spend most of time in calculating the force between each

atom. Adopting Neighbor Lists can cut down the time we spend in calculating [Lambrakos, et al., 1986 and Verlet, et al., 1967]. Neighbor List theory builds the relationship between each atom and its surrounding atoms in a particular interval of time steps as Fig. 2.3 (c), likes ten to twenty steps. Picking the atom as the center,

is the radius of Neighbor List like Fig. 2.4 and rc

rL r_L = + . In the particular r_c σ

interval of time steps, when we calculate the force on the atom , only use the atoms within range to determine which atom exists in the range. And then calculate

the force of atom without calculating the distance from whole atoms between each other at each time step in the system. For example, in the system containing atoms, build a Neighbor List every ten time step in a simulation to establish the relationship of distance between each two atom. The number of calculating time is

. After every atom in the region being confirmed among all atoms in

the simulation, the calculation in forces or energy of the following ten time steps only need to estimate which atom is in the region by the atom within region. And

then, update the distributing information of itself in each atom after passing through rL

rc r_L

the next ten time steps. When we use Neighbor Lists, we should especially notice the choose of σ whose value should not be too less, avoiding the calculation of forces or energy being affected by the atoms which is not within range originally but gets into range during the ten time steps above. On the contrary, if we choose a large value of

σ , the number of atoms in range will increase. Estimating which

atom increases its calculating frequency in the range simultaneouslyincreases the

requirement of time in calculation.

2.4.4 Cell link + Neighbor Lists

We could combine Cell link method and Neighbor List method likes Fig. 2.5 so as to promote the performance of the simulation.

2.5 Boundary conditions 2.5.1 Periodic boundary conditions

Unless the purpose of the MD simulation is to imitate a real walls having physical meaning, walls in simulation had better be eliminated by using Periodic boundary conditions (PBC). The introducing of PBC is equivalent to considering an infinite space-filling array of identical copies of simulation region. Physical meaning of periodic boundary conditions is shown in Fig. 2.6.

2.5.2 Wall boundary conditions

Simulating in the MD system, we would like to keep the wall isothermal, so we define a corrected layer on wall. In the study, we use the Rescaling method to modify corrected layer. Rescaling method keeps the wall isothermal by modifying total kinetic energy. In microcosmic point of view, the temperature is related to kinetic energy. When we set a temperature for corrected layer, it means to set an average kinetic energy of atoms on the corrected layer. In a word, we must keep the kinetic energy fixed (2.16), therefore we should have a reference value. Continually, with the use of (2.17), we compute the total kinetic energy of atoms. Finally, we start rescaling by using (2.18) to make the total kinetic energy in the corrected layer be the same as

the reference value which we use in (2.17).

2.6 Parallel molecular dynamics method

There is no doubt about that MD simulation is a useful and valuable tool. But MD simulation is very time-consuming due to large number of time steps and possibly large number of atoms required to complete a meaningful simulation. In liquids and

the time step to be on the order of fentosecond. Many hundreds of thousand or even millions of time steps are needed to simulate a nanosecond in “real” time scale. In addition, up to hundreds of thousand or millions of atoms are needed in the MD simulation, even for a system size in the nanometer scale. In the past, there have been considerable effort [Plimpton, 1995] that concentrated on parallelizing MD simulation on the memory-distributed machine by taking the inherently parallelism e.g., [Boghosian, et al., 1990] and [Fox, et al., 1998]. Generally, parallel implementation of the MD method can be divided into three categories, including the atom decomposition, the force decomposition and the spatial decomposition among processors.

2.6.1 Atomic-decomposition algorithm

In the atom decomposition method, each processor, which owns nearly the same number of atoms as other processors and in which atoms are not necessarily geometrically nearby, integrates the Newton’s equation for all atoms and moves the atoms of their owns. However, this method requires global communication at each time step, which becomes unacceptably expensive as compared with the “useful” MD computation when the number of atoms increases to a certain amount, since each processor has to know all information (position and velocities) of all atoms at each time step. Or equivalently, the communication is O(N), where N is the number of

atoms in the system that is independent of the number of the processors, P. Thus, the atom decomposition method is generally suitable for small-scale problem

2.6.2 Force-decomposition algorithm

In the force decomposition method, it is based on a block-decomposition of the force matrix rather than a row-wise decomposition in the atom-decomposition method.

It improves the O(N) scaling g to be O(N / P) . It generally performs much better that the atom decomposition method; however, there exists some disadvantages. First, the number of processors has to be square of an integer. Second, load imbalance may become an issue. From previous experience [Plimpton, 1995], it is suitable for small- and intermediate-size problems.

2.6.3 Spatial-decomposition algorithm

In the spatially static domain decomposition method, simulation domains are physically divided and distributed among processors. This method so far represents the best parallel algorithm for large-scale problem in MD simulation for short-ranged interaction [Karypi et. al., 1998]; however, it only works well for a system, in which the atoms move only a very short distance during simulation or possibly distribute uniformly in space. MD simulation of solids represents one of the typical examples.

In contrast, if the distribution of the atoms tends to vary very often in the configuration space, then the load imbalance among processors develops very fast

during simulation, which detriments the parallel performance. Thus, a parallel MD method capable of adaptive domain decomposition may represent a better solution for resolving this difficulty.

2.6.4 PCMD (Parallel Cellular Molecular Dynamics) algorithm

A new parallel algorithm for MD simulation, named parallel cellular molecular dynamics (PCMD), is developed by MuST (Multiscale Science & Technology) laboratory in NCTU in Taiwan, employing dynamic domain decomposition to address the issue of load imbalance among processors in the spatially static domain-decomposition method. We focus on developing a parallel MD method using dynamic domain decomposition by taking advantage of the existing link-cells as mentioned earlier. In this proposed method, not only are the cells used to reduce the cost for building up neighbor list, but also are used to serve as the basic partitioning units. Similar idea has been applied in the parallel implementation of direct simulation Monte Carlo (DSMC) method [Nicol et. al., 1988], which is a particle simulation technique often used in rarefied gas dynamics. Note that in the following IPB stands for interprocessor boundary. General procedures as Fig. 2.7 in sequence include:

1. initialize the positions and velocities of all atoms and equally distribute the atoms among processors;

2. Check if load balancing is required. If required, then first repartition the domain,

followed by communicate cell/atom data between processors, renumber the local cell and atom numbers, and update the neighbor list for each atom due to the data migration;

3. Receive positions and velocities of other atoms in the neighbor list for all cells near the IPB;

4. Compute force for all atoms;

5. Send force data to other atoms in the neighbor list for all cells near the IPB;

6. Integrate the acceleration to update positions and velocities for all atoms;

7. Apply boundary conditions to correct the particle positions if necessary;

8. Check if preset total runtime is exceeded. If exceeded, then output the data and stop the simulation. If not, check if it is necessary to rebuild the neighbor list of all atoms using the most update atom information.

9. If it is necessary to rebuild the neighbor list (N=8 in the current study), then communicate atom data near the IPB and repeated the steps 2-8. If not necessary, then repeat steps 3-8.

In the above, in addition to the necessary data communication when atoms cross the IPB and particle/cell data near the IPB, there are two more important steps in the proposed parallel MD method as compared with the serial MD implementation. One is how to repartition the domain effectively and the other is the decision policy for

repartitioning. These two steps are described next, respectively.

Chapter 3 Simulation of xnone and helium droplet-droplet collision dynamics

3.1 Simulation conditions

This content mainly concerned the simulation of the behavior of xenon or helium droplet pair collision. In order to compare the influence of crucial parameter, impact parameter and the relative velocity between two droplets in the collision, the diameter of single droplet was set to be similar about 10.5nm in both xenon and helium. In detail, in the composition of xenon droplet, we arranged the L-J potential atom in FCC crystal structure at first, which density was 0.2 and maintained equilibrium at 166.73K, and then used the PCMD to simulate the variation in a hundred thousand time steps, and eventually obtained a xenon droplet which contained 2800 atoms. By the same approach, the density of helium crystal structure was 1.0 and it maintained equilibrium at 0.55K. After using PCMD in a hundred thousand time steps, we would obtain a droplet containing 2700 atoms. To establish the simulation model of droplet pair collision, we copied the relative coordinates of droplets and adjusted the relative velocities between them. In addition, if the test case was head-on collision, the system dimension should be 200σ *500σ *500σ (Fig. 3.1). For detailed observation of the behavior of disruption or fragmentation, otherwise the system dimension was usually

600σ *300σ *300σ (Fig. 3.2).

3.1.1 Test conditions

When simulating a collision of xenon droplet, the impact parameter was 0 and the relative velocity ranged from 1250m/s to 2250m/s in a head-on condition, but in a non-head-on condition, the impact parameter was 1.25nm to 8.75nm and the relative velocity ranged from 250 m/s to 2250m/s. In simulating a collision of helium droplets, the impact parameter was 0, but in a non-head-on condition, the impact parameter was

在文檔中應用分子動力學模擬探討不同分子量對奈米尺寸液滴碰撞之影響 (頁 26-0)