An improved smoothed molecular dynamics method by alternating with molecular dynamics

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294 www.elsevier.com/locate/cma

An improved smoothed molecular dynamics method by alternating with molecular dynamics Nianfeng He, Yan Liu ∗ , Xiong Zhang School of Aerospace Engineering, Tsinghua University, Beijing 100084, PR China Received 23 January 2015; received in revised form 8 August 2015; accepted 10 August 2015 Available online 17 August 2015

Highlights • • • •

An alternating smoothed molecular dynamics (AltSMD) method is proposed. AltSMD can adopt much larger time step size than MD and achieve accurate results. Influences of parameters in AltSMD method are investigated in detail. A criterion for automatic alternating is proposed for AltSMD method.

Abstract An effective method for nano-mechanics is presented by combining smoothed molecular dynamics (SMD) method and molecular dynamics (MD) method. SMD method was proposed in our previous work to sharply increase the feasible time step of MD. A set of background mesh was introduced in SMD and the equations of motion are solved on the background mesh nodes rather than on atom sites as in MD, which converts the controlling factor of critical time step size to the background mesh size. Much larger time step size, which can be one order larger than MD step size, can be adopted in SMD computation. But the application of SMD into problems with moving atom disorders and at finite temperature is still not satisfactory. An improvement scheme is proposed in this paper by alternating SMD computation and MD relaxation in the solution process. SMD is used at the beginning of the computation, then it is converted to MD relaxation whenever required. SMD computation continues after MD relaxation until next relaxation is needed. The conversion between MD and SMD is very convenient and straightforward owing to their similarities. The accuracy can be guaranteed by MD relaxation process, and the overall efficiency is better than MD efficiency since large time step size can be adopted in SMD computation. Examples of nano-indentation and nanowire tension validate the alternating SMD method. Influences of factors including the length of relaxation intervals and the number of MD relaxation steps are studied. A criterion automatically determining the alternating process is also proposed and validated. c 2015 Elsevier B.V. All rights reserved. ⃝

Keywords: Molecular dynamics; Smoothed molecular dynamics; Multiscale computation; Nano-indentation; Nanowire

∗ Corresponding author.

E-mail address: [email protected] (Y. Liu). http://dx.doi.org/10.1016/j.cma.2015.08.005 c 2015 Elsevier B.V. All rights reserved. 0045-7825/⃝

274

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

1. Introduction Molecular dynamics (MD) method is widely used in many scientific research fields such as fracture mechanics [1], nanowire mechanics [2,3], composites [4], and problems of extreme condition [5]. MD results have straightforward physical meaning, and a large amount of important information can be extracted. The applicable spatial scale and temporal scale of classical MD method, however, are still very limited even with fast-developed massive computing devices. The limitation in spatial scale lies in that one point in MD method usually represents one atom, which can result in a system of millions of degrees of freedom (DOFs) in one computation. It is very costly to simulate as little as 1 µm3 ˚ Concurrent multiscale methods, which combine material, because the typical distance between atoms is of the order A. molecular dynamics and continuum-based modeling techniques, have been focused on to expand the spatial capability of MD. A typical concurrent multiscale method includes two kinds of regions. The core atomic region surrounds the crack tips, the dislocations or the contact lines where traditional continuum theory may fail. The continuum method is used in the far-field region. The computational cost in the far-field region therefore can be saved owing to the reduction of atomic DOFs. Handshaking region, or called the bridging region or the transition region, links the atomic region and the continuum region. The idea of concurrent multiscale method can be traced back to the work of Kohlhoff et al. [6] in early 1990s. Other pioneer work includes Tan and Yang [7], Broughton et al. [1], Shilkrot et al. [8], E et al. [9], and the quasi-continuum method [10]. Rising interests in nano-mechanics and nano-electro-mechanical system encourage researchers to develop more accurate multiscale methods. Huge differences between discrete atomic system and continuum system may lead to mismatch on the interface of the two regions. Liu and his collaborators [11,12] proposed a bridging scale method, which can eliminate the wave reflection if high-frequency waves approach the interface. Xiao and Belytschko [13] proposed another method, called the bridging domain method, which bridged the two kinds of regions with Lagrangian multiplier and augmented Lagrangian multiplier methods. The mismatch between atomic region and continuum region was also recognized as “ghost force” in quasi-continuum method [10,14], which can be carefully eliminated from force terms. Fish et al. [15] blended the continuum equilibrium equation and the atomic force in the transition region, and the constraints in the interface region were implemented in the weak form. Bridging the giant gap can also be dealt with by keeping the atomic potential or forces whilst reducing atomic DOFs by introducing continuum assumptions. The coarse-grained molecular dynamics [16] is a pioneer work in this field. The Cauchy–Born rule was adopted in quasi-continuum method [10] to convert atomic potential to the function of local deformation gradient. Quasi-continuum method was further improved by sampling atomic potential in clusters [17]. Liu et al. [18] proposed atomic finite element method (AFEM), which directly computed nodal forces in finite element method based on atomic potential. A linkage between atomic finite element and conventional finite element was also developed by non-local interaction. Biyikli et al. [19] proposed a multiresolution molecular dynamics framework based on sampling the atomic energy in finite elements. They carefully classified different kinds of representative atoms, sampling atoms and non-sampling atoms to improve the accuracy. Recent investigations of concurrent multiscale methods focused more on finite-temperature problems. Dupuy et al. [20] developed finite temperature quasi-continuum method. The widely used Cauchy–Born rule in quasicontinuum method and other multiscale methods needs to be adjusted at finite temperature, and a quasi-harmonic assumption is always adopted to simplify the formulation. A finite temperature continuum theory [21] was constructed from finite temperature Cauchy–Born (FTCB) rule and Taylor expansion of atomic potential. Liu and Li [22] proposed a non-equilibrium multiscale model based on non-equilibrium MD and FTCB rule. Liu and Li [23] further combined FTCB rule and cohesive element to construct a multiscale interphase zone model, and applied the model into dynamic fracture problems. Fish et al. [24] proposed a generalized mathematical homogenization (GMH) model for finite temperature based on investigation of dynamic problem on unit cell coupled with thermo-mechanical continuum equations. Xiang et al. [25] investigated the atomic motions in a representative volume element (RVE) to construct thermo-mechanical constitutive equations. A separation of mechanical motion and thermal motion was employed, and the equations were constructed based on conservation laws. Originating from the heterogeneous multiscale method (HMM) [9], a concurrent method at finite temperature was developed and then successfully applied for brittle crack problems [26]. Anciaux et al. [27] focused on the bridging domain method, and they removed the thermal gradient in the continuum region by improving the coupling technique. Mathew et al. [28] separated the mechanical

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

275

vibration and thermal oscillation. Only the low-frequency vibration were allowed into continuum region, and the energy of high-frequency motion was conveyed into the continuum region by thermal flux. The frequency separation can be finished by low-frequency filter. Ramisetti et al. [29] suggested an effective digital filter into the multiscale method. Even with well-developed concurrent multiscale methods, the problem scale is still limited due to expensive computational cost in MD region. Typical MD step size is at the order of femtoseconds, which is the result of stability requirement of explicit integration in time domain. The critical time step size in MD computation is restricted by highfrequency vibrations. Increasing feasible time step size, which in sequence increases the capability of MD to simulate larger time scale, is required. Methods for larger time step size include constraint method [30], subspace method [31], and multiple-time-step (MTS) method [32–34]. Constraint method focused on the simulation of large molecules. High-frequency vibrations of bonds were suppressed in constraint method [30], which released the restrictions on the critical time step size. Similar idea was adopted in subspace method [31], where an optimal cut-off frequency was determined and only the subspace containing low-frequency motions was reserved. The integrations of low-frequency motions and high-frequency motions were divided in MTS method [32], and the time step for low-frequency motions can be several times larger than that for high-frequency motions. A double MTS method [34] was proposed by further categorizing the vibrations of light and heavy atoms in fast motions. The above methods can effectively increase the available time step size, but the enlargement of time step size is still limited. Qian et al. [35] further extended MTS method by coupling enriched space–time finite element method with molecular dynamics, which has been successfully applied to dynamic fracture of lattice structures. In our previous work [36], smoothed molecular dynamics (SMD) method was proposed by introducing one set of background mesh into MD method. Much larger time step size, even one order of magnitude higher than MD critical time step size, can be adopted in SMD. The idea of SMD method originated from material point method (MPM). An Eulerian background mesh is used along with Lagrangian material points in MPM. The background mesh serves to solve equations of motion, which brings an essential advantage into MPM that the critical time step size is controlled by background mesh size instead of characteristic length between Lagrangian points. Even undergoing very large deformation, the critical time step size will not sharply decrease as the background mesh does not deform, so that MPM can be very efficient in large deformation problems. Numerical examples confirmed that MPM can be much more efficient than meshfree SPH method and even more efficient than explicit finite element method (FEM) [37]. MPM has been successfully applied in problems such as dynamic fracture [38], high-velocity impact [39,5], and fluid–structure interaction [40]. MPM was also adopted in multiscale simulation. Guo and Yang [41] and Lu et al. [42] coupled MPM to MD method to construct concurrent multiscale methods. Shen and Chen simulated delamination problem of film using both MD and MPM [43]. Liu et al. [5] proposed a multiscale framework by sequentially combining MD and MPM for hyper-velocity impact problem. Accurate parameters for equations of state under extreme condition can be extracted, and then they will be adopted in MPM computation. Similar to MPM, the atomic equations of motion are mapped onto background mesh nodes and solved there in SMD method. Then the atomic variables are updated by nodal variables. SMD successfully changes the factor controlling critical time step size to background mesh size. If an element size much larger than lattice constant is adopted, the available time step size in SMD can be much larger than that in conventional MD [36]. The enlargement of time step can be attributed to smoothing out high-frequency motions of MD results, which will not violate the overall accuracy if the low-frequency deformation is focused on. Parallel version of SMD method was developed [44]. Wang et al. [45] proposed adaptive SMD method where the background mesh can be hierarchically refined. The resolution of local atom disorder in SMD, however, is not good due to smoothing out high-frequency motions. A straightforward way is to couple SMD with MD [36]. The two methods can be coupled very easily and efficiently owing to their similarity. Coupled MD–SMD method can capture plastic yield point and local atom disorder correctly with less computational resource. SMD method is much improved in this paper based on the key idea of prediction–correction method. The solution process is accelerated with standard SMD, and the accuracy is improved by MD relaxation. So the improved method can obtain results close to MD results while keeping high efficiency. Section 2 introduces the theory of SMD method. The improved SMD method, called the alternating smoothed molecular dynamics (AltSMD) method, is proposed in Section 3. AltSMD method is also validated with nano-indentation problem in Section 3. The criterion controlling alternating between SMD computation and MD relaxation is then developed in Section 4. Section 5 further validates the new method with tension problems of nanowire. The entire paper is concluded in Section 6.

276

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

2. Smoothed Molecular Dynamics (SMD) 2.1. Brief introduction to molecular dynamics The equation of motion of an atom in conventional MD can be written as [46] m i r¨i = Fi = Fi int + Fi ext

(1)

where the subscript i is the atom index, m i is the atom mass, ri is the position vector, the superposed dots denote differentiations with respect to time, and Fi is the total force exerted on the atom i. The superscripts ext and int represent the external part and the internal part, respectively. The internal force, which is the interaction between atoms, is calculated by Fi int = −

∂ Ep ∂ri

(2)

where E p is the potential energy. Eq. (1) is numerically integrated in time domain. MD simulations are involved in large number of atoms, which hinders the use of implicit integration schemes solving a set of linear equations at every time step. Explicit integration schemes, such as various forms of central-difference algorithm, are popular to solve Eq. (1). For example, the velocityVerlet algorithm updates equations of motion from step n to n + 1 as follows, ∆tMD n r¨ (3) 2 i n+1/2 rin+1 = rin + ∆tMD r˙i (4) ∆tMD n+1 n+1/2 + r˙in+1 = r˙i r¨ (5) 2 i where ∆t is the time step size, and the subscript MD is used to distinguish with SMD time step in the context. ∆t MD is limited by atom spacing. Numerical disasters may be aroused if inappropriate time step size is chosen. The potential plays an important role to obtain realistic and accurate results. Embedded atom method (EAM) potential [47] is adopted in the following examples, n+1/2

r˙i

Ei =

= r˙in +

1 2

n 

Φi j (ri j ) + ψi (ρi )

(6)

j=1, j̸=i

where Φi j is the contribution of pair interaction between atom i and atom j, ψi is the embedding energy of atom i which reflects interaction with the environment of electrons, and ρi is the electron density. The readers can refer to literature [47] for details. 2.2. Theory of SMD Inspired by the idea of material point method, a background finite element mesh was introduced in SMD method. The element size of background mesh is usually chosen much larger than atomic lattice distance, but theoretically speaking any mesh size can be used. In fact, if the element size is so small that any two elements containing atoms do not share nodes, SMD will degenerate to classical MD [36]. The mesh is usually uniform cubic mesh in 3D and uniform square mesh in 2D though any type of mesh can be used without difficulty. All the physical variables are still carried and recorded on atoms. The difference from MD is that the equations of motion are constructed and solved on background mesh nodes instead of on atom sites and the atomic variables are updated by nodal increments. The background mesh nodal variables and atomic variables are linked through standard finite element (FE) approximation gi =

ne  J =1

N J (ri ) g J

(7)

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

277

where g can be any variable such as displacement, velocity, and acceleration. N J (ri ) is standard FE shape function, and n e is the number of nodes of the element atom i is in. The lowercase subscript indices i and j refer to atoms and the uppercase ones refer to nodes. Inside each time step, the atoms are bound to the background mesh and they deform together. Substituting Eq. (7) into Eq. (1) leads to mi

ne 

N J (ri )¨r J = Fi .

(8)

J =1

The functionality of atoms in SMD method is similar to that of quadrature points in finite element method. The momentum equations of background mesh nodes are assembled from atomic equations of motion, that is, the momentum equation of one node is constructed by mapping Eq. (8) of all the atoms inside the linking elements onto the node with finite element shape functions,   na ne na    N I (ri )m i N J (ri )¨r J = N I (ri )Fi , I = 1, 2, . . . , n node , (9) i=1

J =1

i=1

or in the concise form ne 

M I J r¨ J = F I ,

I = 1, 2, . . . , n node ,

(10)

J =1

where MI J = FI =

na  i=1 na 

N I (ri )N J (ri )m i I3×3

(11)

N I (ri )Fi

(12)

i=1

n node is the number of total nodes, n a is the number of atoms influenced by node I , and I3×3 is the identity matrix of order three. Deriving Eq. (10) can also be understood as a least-square solution of Eq. (8), since the number of independent variables (i.e. the number of nodal variables) is smaller than the number of equations in (8). Eq. (9) is actually the normal equation of Eq. (8). Eq. (11) is usually referred as consistent mass matrix, which requires solving linear algebraic equations in every time step. The lumped mass matrix is much more popular in explicit solver since no linear system needs to be solved. The lumped mass matrix can be obtained by accumulating all the elements in one row of the consistent mass matrix to the diagonal position and dropping the other matrix elements, that is M I I = M I I3×3 =

n node J =1

MI J =

na 

N I (ri )m i I3×3 ,

(13)

i=1

  M = diag M11 , M22 , . . . , Mn node n node (14) n node where the property of partition of unity J =1 N J (ri ) = 1 for finite element shape function is invoked in deriving the third equal sign in Eq. (13). The nodal momentum equations (10) then change to the following form by using lumped mass matrix, p˙ I = M I r¨ I = F I ,

I = 1, 2, . . . , n node .

(15)

The nodal momentum p I is calculated by accumulating atomic momenta to the background mesh nodes similar to the assembling process of nodal forces pI =

na  i=1

N I (ri )m i r˙i .

(16)

278

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

The flowchart of one SMD time step is listed below. It is assumed that the variables at and before nth step are known, and the variables at (n + 1)-th step are to be solved. The velocity-Verlet scheme is shown as an example, though any other explicit integration scheme can be adopted. 1. For each atom, determine ID of the element containing the atom. Calculate the values of nodal shape functions on the atom site. The nodal masses M In , the nodal forces FnI and the nodal momenta pnI at the beginning of step n + 1 are calculated respectively by Eq. (13), Eq. (12) and Eq. (16). The superscript n is used to denote that all the above variables are calculated by the values at step n. 2. Update nodal momenta and atomic velocities at intermediate step n + 1/2 with nodal forces n+1/2

pI

n+1/2

r˙i

∆t SMD n FI 2 ne N I (rin )FnI ∆t SMD  = r˙in + 2 M In I =1 = pnI +

(17) (18)

where ∆t SMD is SMD time step size. 3. Update atomic positions with interpolated nodal momenta rin+1

=

rin

n+1/2 ne  N I (rin )p I + ∆t SMD M In I =1

(19)

and then calculate atomic forces Fin+1 at new positions. 4. The new atomic forces are mapped to the nodes to construct new nodal forces Fn+1 by Eq. (12). The shape function I N I is not changed inside one step since the atoms are assumed to be bound and deform with the background mesh. 5. Update atomic velocities from intermediate step n + 1/2 to step n + 1 with interpolated nodal forces n+1/2

r˙in+1 = r˙i

+

ne N I (rin )Fn+1 ∆t SMD  I . n 2 M I I =1

(20)

6. Abandon the deformed background mesh at the end of the step. The undeformed background mesh will be used again at the beginning of next step. It should be noted that atomic velocities are updated by nodal forces instead of new nodal velocities, which is similar to the updating process in material point method to avoid artificial damping [48]. Similarly, atomic displacements are updated by nodal momenta instead of nodal displacements. The above flowchart is very similar to MD flowchart except for the mapping process between atoms and background mesh nodes (i.e. atomic equations are assembled to nodal sites and nodal increments are interpolated to update atomic variables). Typical techniques in MD method, such as thermostat or velocity scaling for NVT ensemble, can be implemented directly in SMD method. 3. Alternating Smoothed Molecular Dynamics Method 3.1. Theory and procedure As was shown in our previous work [36], SMD method can use much larger time step because high frequency atomic motions are filtered out during the assembling and interpolation process. The name “smoothed” also came from the smoothing feature of the assembling and interpolation process. If local deformation gradient is small and no local atomic disorder happens, SMD can well predict the behavior of atom system and the computational cost can be greatly saved. If local deformation is very large and dominates the major feature of the whole problem, SMD results may deviate from MD solution and the deviation increases as deformation grows. Fig. 1 shows the curves of indentation force versus indentation depth of a nano-indentation problem. SMD curve is higher than MD curve because the indentation force spreads a larger area due to the smoothing process. But SMD method can obtain a nice approximate prediction of the local deformation. As SMD method can use much larger time step size, it can be adopted as a fast integrator to an approximate status close to real MD results. Based on

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

279

Fig. 1. SMD and MD results of nano-indentation problem.

this point and inspired by the idea of predictor–corrector algorithm, we propose an alternating smoothed molecular dynamics (AltSMD) method as sequentially alternating SMD computation and MD relaxation in time domain. The procedure of AltSMD method is shown schematically in Fig. 2. SMD method is used at the beginning of computation. State variables are recorded and monitored during SMD computation. If the variables meet some criteria showing that large deviation may exist between SMD results and MD results, the computation is alternated to MD relaxation. The MD relaxation proceeds until that the computational results have been corrected to real status. After MD relaxation, the computation continues with SMD method until another MD relaxation is needed. SMD computation and MD relaxation are alternated in time domain throughout the computation. The key idea of AltSMD method can be understood as a predictor–corrector algorithm with one cycle of SMD computation and MD relaxation as one huge prediction–correction step. SMD computation can be viewed as an efficient predictor with very large time steps, which can approach the real status much faster than MD computation. MD relaxation can be identified as a corrector which pulls the inaccurate status to an accurate status. AltSMD possesses both the advantages of MD and SMD. MD relaxation ensures AltSMD method much more accurate than standard SMD method. SMD stepping process brings much better efficiency since MD relaxation only occupies a very small part of computational cost. The implementation of AltSMD method is very simple and straightforward owing to the similarity between MD and SMD. In fact, nearly no additional computation cost is needed in the alternating process. If SMD computation is “on”, then the background mesh is used and all the momentum equations are constructed and solved on mesh nodes. If MD relaxation is “on”, then the external loading is kept unchanged and all the momentum equations are constructed and solved on atom sites. Both SMD method and AltSMD were implemented in the code LAMMPS [49], with which all the following examples were computed. The performance of AltSMD at fixed ratio of relaxation steps is emphasized in this section, and the specific criterion will be discussed in Section 4. 3.2. Validation of AltSMD method with nano-indentation The first validation example is an ideal copper block indented by a soft indenter, as shown in Fig. 3. The alternating process is executed by fixing SMD and MD intervals in this subsection. The dimensions of the copper block are ˚ is the lattice constant of copper. A spherical indenter is placed above the 20a0 × 20a0 × 20a0 , where a0 = 3.615 A ¯ direction. Periodical boundary conditions are applied in [1 0 0] (X ) copper block and moves downward along [0 0 1] and [0 1 0] (Y ) directions, and free boundary conditions are applied on the top face. The atoms in the bottom two layers are rigidly fixed. The amplitude of repelling force between the indenter and an atom  −K (ri − R)2 for ri ≤ R F(ri ) = (21) 0 for ri > R

280

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 2. Schematic flow chart of AltSMD method.

Fig. 3. Schematic diagram of nano-indentation example.

˚ 3 . The where ri is the distance between atom i and the center of the indenter. Constants R = 5a0 and K = 10 eV/A system is firstly relaxed for 10 ps, then the indenter moves with constant velocity 1 m/s. The indenter is kept stationary during MD relaxation process in AltSMD simulation. The background mesh size is 5a0 × 5a0 × 5a0 . ∆t MD = 0.01 ps and ∆t SMD = 0.04 ps. AltSMD results are shown in Figs. 4–7 along with MD results. 300 MD relaxation steps are inserted between every 250 SMD steps in the first calculation. The calculated indentation forces during both SMD stepping and MD relaxation are plotted, so that it can be seen clearly that MD relaxation can pull SMD results back to real force–depth

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

281

Fig. 4. AltSMD results with 300 MD relaxation steps between every 250 SMD steps.

a

b

c

Fig. 5. Comparison of configurations at various indentation depths. (a) MD results; (b) SMD results; (c) AltSMD results. Color denotes CSP value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

curve, which shows “zig-zag” segments in Fig. 4. The zig-zag segments are not shown in the following figures and time average is applied to indentation force to remove oscillations. Fig. 5 shows the configurations of the middle slice. AltSMD method can capture local atom disorder right under the indenter, and the configuration is close to MD result in both the elastic regime (indentation depth 0.4nm) and the plastic regime (indentation depths 0.6nm and 0.8nm). The deformation under the indenter (region in dashed box in Fig. 5) is magnified in Fig. 6. Disorder atoms are plotted in Fig. 7. Color in Figs. 5–7 denotes the centro-symmetry parameter (CSP) value, which is widely used for indication of atom disorders. Both CSP distribution and defect structures of AltSMD results agree well with those of MD results. The influences of the number of SMD computation steps (i.e. MD relaxation intervals) and the number of MD relaxation steps are investigated. Firstly the number of SMD computation steps are fixed at 250, and the number of MD relaxation steps are changed from 10 to 300. Typical results are shown in Fig. 8. The legend “n SMD : n MD ”

282

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 6. Comparison of deformations under the indenter. Color denotes CSP value.

˚ 2 are plotted. The indentation depth is 0.8 nm. Fig. 7. Comparison of defects under the indenter. Only the atoms whose CSP value larger than 0.4 A

Fig. 8. Influences of the number of MD relaxation steps (n MD ).

stands for that “n MD ” MD relaxation steps are inserted between every “n SMD ” SMD computation steps. It can be seen that even a very small number of MD relaxation steps can avoid overestimating the repelling force, but the plastic

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

283

Fig. 9. Influences of the number of SMD computation steps between MD relaxations (n SMD ).

Fig. 10. Combination of SMD and AltSMD methods.

yield point, which is indicated as a sudden drop on the curve, cannot be captured if the number of MD relaxation steps is too small. A clear yield point is obtained when n MD is 100. Increasing n MD leads to better prediction of yield point, and AltSMD results with n SMD : n MD = 250 : 300 agree well with MD results even in the plastic regime. Then n MD is fixed at 300, and n SMD is varied to investigate its influence. The results are shown in Fig. 9. Increasing n SMD does not influence the capture of yield point much. The yield point can be well predicted even when n SMD is large if enough MD relaxation is executed. But too large n SMD will decrease the accuracy in the plastic regime, which implies that timely alternating controlled by criterion may be needed for the plastic regime. An extreme case is then examined. SMD method is used in the first half of solution process, and then AltSMD method is used for the remaining solution process. The results are shown in Fig. 10. It is appealing to see that AltSMD method can obtain nice results even when the first half of the curve deviates much from MD results. This case ensures the correctness of using criteria to automatically alternate between MD relaxation and SMD computation, where SMD computation may last a long time before the criterion is satisfied.

284

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294 Table 1 Comparison of CPU time for fixed n SMD : n MD . n SMD : n MD

CPU time (s)

Ratio to MD cost

Estimated ratio (Eq. (22))

MD method 250 : 10 250 : 100 250 : 200 250 : 300 600 : 300 900 : 300

7351 2601 3351 4145 4971 3389 3062

1.0 0.354 0.456 0.564 0.676 0.461 0.417

— 0.36 0.45 0.55 0.65 0.48 0.43

Fig. 11. Indentation at fixed temperature 100 K.

Compared to typical MD step, SMD step has extra computation for the assembling and interpolation process. The ratio of AltSMD computational cost to MD computational cost for fixed n SMD : n MD can be roughly estimated as η=

n SMD · (1 + α) + n MD t SMD n SMD ∆ ∆t MD

(22)

where α is the ratio of the additional cost in one SMD step to the total cost in MD step. For the present indentation example, α is about 40%. The computational cost of AltSMD method is listed in Table 1. The total number of steps of reference MD computation is 80,000. The estimation equation (22) agrees well with the practical cost. Computational saving is from 32% to 64% depending on n SMD : n MD , which is not very high because ∆t SMD is only four times of ∆t MD . The efficiency will be increased much if a much larger background mesh size is used. Performances of AltSMD at finite temperature are then investigated. The results for the temperatures 100 K, 200 K and 300 K are shown in Figs. 11–13. The canonical (NVT) ensemble is used in these examples. Conclusions similar to the above NVE examples are obtained. Even a very small n MD can achieve good curves in elastic region, but cannot capture plastic yield point. When n MD increases, AltSMD method can capture better yield point, and its results are close to MD results in both elastic and plastic regions. 4. Alternating Criterion for AltSMD Method As shown in the previous section, SMD force–depth curve can be corrected to MD curve well in the elastic region even when large difference occurs. But the alternating process should be more frequently invoked before the plasticity mechanism occurs to correctly capture the yield point. So the criterion controlling the alternating process should indicate whether local atom disorder or local large deformation may happen.

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

285

Fig. 12. Indentation at fixed temperature 200 K.

Fig. 13. Indentation at fixed temperature 300 K.

A criterion is proposed based on the difference of SMD and MD computations. One major difference is that MD updates atomic velocities with real atomic interaction forces, while SMD updates atomic velocities with nodal forces. The assembling and interpolation process smoothes out high-frequency part of atomic forces, which can produce an indicator of local atom disorder or local large deformation. This is because when the deformation is uniform and gradual, the difference between forces before and after the assembling and interpolation process should be small; otherwise the difference can be large due to smoothing out local atom disorders. The criterion is in the following form εE =

|Fi − Fi∗ |ave−E ≥ε |Fi − Fi∗ |ave−all

(23)

where ε is the threshold value, Fi is MD atomic force, and Fi∗ is SMD force, which is defined as the product of the interpolation of nodal accelerations and the atom mass   ne  ∗ Fi = N I (ri )¨r I · m i . (24) I =1

286

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 14. Results of indentation example with different criterion threshold values.

Fig. 15. Dimensionless MD relaxation steps for different criterion threshold values.

The subscript “ave-E” means averaging over atoms in one element E, and “ave-all” denotes averaging over all the atoms. Spatial average of force differences is used instead of difference of single atom in Eq. (23) to avoid too frequent alternating between SMD computation and MD relaxation. It should be noted that the elements to calculate criterion can be different from the background mesh in SMD computation. The threshold value ε is usually obviously larger than unity. If the criterion Eq. (23) is satisfied in any element, which means that local atom disorder may happen, SMD computation should be converted to MD relaxation. The nano-indentation example is examined with the above force criterion. Minimum and maximum values of n MD and n SMD are set to avoid missing correct yield point or too frequent alternating. For the current example, at least 5 MD steps must be executed if the relaxation is invoked, and the number of steps per relaxation is not allowed more than 300. The results are shown in Fig. 14. It is not surprising to see that the computation with criterion controlling can obtain better results since more frequent alternating and larger number of relaxation steps can be automatically invoked around the yield point. Fig. 15 shows dimensionless MD relaxation steps (n ∗MD ) for different criterion values. n ∗MD is calculated by dividing total number of relaxation steps by total computation steps in pure MD simulation. Larger criterion value (ε = 4.0) can significantly reduce the number of relaxation steps and computational cost for

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

287

Table 2 Comparison of CPU time for criterion-controlled AltSMD computations. Method

CPU time (s)

Ratio to MD cost

MD AltSMD, n SMD : n MD = 250 : 300 AltSMD, n SMD : n MD = 900 : 300 AltSMD, ε = 3.0 AltSMD, ε = 3.5 AltSMD, ε = 4.0

7351 4971 3062 5200 3113 2727

1.0 0.676 0.417 0.707 0.423 0.371

Fig. 16. Tension of an initially indented nanowire.

MD relaxation, but the yield point delays a little and the curve in the plastic regime deviates more from MD curve. Results of ε = 3.0 and 3.5 accurately capture the yield point, and they are close to MD results in the plastic regime. We also tested criterion values 2.0 and 2.5, which can also give accurate results but result in too many relaxation steps. The efficiency is provided in Table 2, and the efficiency for fixed n SMD : n MD is also listed for comparison. The computation with criterion ε = 3.5 is as efficient as that with fixed n SMD : n MD = 900 : 300, but its accuracy is much better. The choice of criterion value is a balance between computation accuracy and computation burden, and the optimal value may depend on problems. For the present problem of nano-indentation the value around 3.5 is a good choice. 5. Examples of Nanowire Tension 5.1. Nanowire with initial indentation Two kinds of nanowire tension examples are then calculated to further validate AltSMD method. The first tension specimen is a copper cuboid with two initial cylindrical indentations as shown in Fig. 16. The dimensions of the cuboid are 100a0 × 20a0 × 8a0 . The two cylindrical indentations are of the radius 10a0 and the cylinder centers are located 7a0 away from the cuboid edge. Periodic boundary conditions are applied in the thickness direction, and free boundary conditions are applied in the transverse direction. Two layers of atoms at each end of the specimen are specified with constant velocity v = 0.90375 m/s. The model is first relaxed for 10 ps. The size of background mesh is 5a0 · ∆t MD = 0.01 ps and ∆t SMD = 0.04 ps. AltSMD results with different n MD are demonstrated in Fig. 17. Before dislocation appears, AltSMD results are identical to MD results no matter how many MD relaxation steps are used. When n MD is large enough, the plastic yield point, denoted as a sudden drop in the tensile force–engineering strain curve, can be accurately captured. We also tested algorithms with fixed n MD but different n SMD . The observation, as shown in Fig. 18, is similar to that in the nano-indentation example. Enough MD relaxation steps can lead to good yield point even when n SMD is large, but too large n SMD will decrease the accuracy for plastic deformation. The results for finite temperature 300 K are shown in Fig. 19, the results are also close to MD results. 5.2. Sharply-notched nanowire The second model of nanowire tension is a notched specimen as shown in Fig. 20, which is a 100a0 × 20a0 × 5a0 copper cuboid with an initial crack of the length 5a0 . The interaction forces between atoms on the two sides of the crack are excluded to implement discontinuity. The boundary conditions are the same with the previous tension

288

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 17. Comparison of MD and AltSMD results with different relaxation steps.

Fig. 18. Comparison of MD and AltSMD results with different relaxation intervals.

example. The specimen is firstly relaxed for 10ps and then stretched to broken. The background mesh size is 5a0 · ∆t MD = 0.01 ps and ∆t SMD = 0.04 ps. AltSMD results with different n MD and n SMD are shown in Figs. 21 and 22, respectively. Within elastic stage, AltSMD results are almost the same with MD results even with very few MD relaxation steps. When n MD ≥ 100, AltSMD method can capture well the development of crack. Much larger n SMD still leads to good results with enough relaxation steps. Snapshots of crack propagation of MD results and AltSMD results with n SMD : n MD = 250 : 200 are shown in Fig. 23. The region around crack tip (the dashed box in Fig. 23) is magnified in Fig. 24. AltSMD configurations are similar to MD results. AltSMD results with different criterion thresholds are shown in Fig. 25. The broken point can be well captured with ε = 3.0 and 3.5. As shown in Fig. 26, the total number of MD relaxation steps needed in this example is only less than 1.5% of the total steps in MD simulation. The proposed criterion can reasonably invoke relaxation. More relaxation work is executed around broken point, but the relaxation work prior to and after the broken point is not much. The computational cost of AltSMD method with different parameters is listed in Table 3. The total number of steps of reference MD computation is 200,000. Computations with alternating criterion are more efficient than those with

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

289

Fig. 19. AltSMD results at fixed temperature 300 K.

Fig. 20. Tension of an initially sharply-notched nanowire.

Fig. 21. Comparison of MD and AltSMD results with different relaxation steps.

fixed n SMD : n MD . Only one-third MD computational cost is required for AltSMD method with criterion, but MD results and AltSMD results are very close. The background mesh size in all the above examples is only 5a0 due to limited specimen size. Larger mesh size will result in larger critical time step size and then higher efficiency. A larger pre-cracked specimen of the dimension 300a0 × 60a0 × 5a0 is calculated. The mesh size is 15a0 , and ∆t SMD = 10 ∆t MD = 0.05 ps. The prescribed velocity

290

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 22. Comparison of MD and AltSMD results with different relaxation intervals.

Fig. 23. Comparison of snapshots during crack propagation. n SMD : n MD = 250 : 200. (a1), (b1) and (c1) are MD results; and (a2), (b2) and (c2) are AltSMD results. Color denotes CSP value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of two ends is v = 1.8075 m/s. Again results coinciding with MD results are obtained as shown in Fig. 27. The efficiency for the larger example is listed in Table 4, which demonstrates six times acceleration to MD computation. 6. Conclusion SMD approach can adopt much large time step size compared to MD approach because the factor controlling critical time step size is converted to background mesh size. An improvement of SMD method, called AltSMD method, is proposed by alternating with MD relaxation process. SMD computation is used throughout AltSMD solution, and the computation is converted to MD relaxation when the difference between SMD results and MD results exceeds threshold value.

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

291

Fig. 24. Comparison of deformations around crack tip. Color denotes CSP value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 25. AltSMD results with different criterion threshold values. Table 3 Computational cost of nanowire tension example. Method

CPU time (s)

Ratio to MD cost

MD AltSMD, n SMD : n MD AltSMD, n SMD : n MD AltSMD, n SMD : n MD AltSMD, n SMD : n MD AltSMD, n SMD : n MD AltSMD, n SMD : n MD AltSMD, ε = 3.0 AltSMD, ε = 3.5 AltSMD, ε = 4.0

22986 7831 10081 12558 14023 11077 9815 7699 7652 7595

1.0 0.341 0.439 0.546 0.610 0.482 0.427 0.335 0.333 0.330

= 250 : 10 = 250 : 100 = 250 : 200 = 300 : 300 = 600 : 300 = 900 : 300

292

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

Fig. 26. Dimensionless MD relaxation steps for different criterion threshold values.

Fig. 27. AltSMD results with different criterion threshold values for the larger example. Table 4 Computational cost of the larger nanowire example. Method

CPU time (s)

Ratio to MD cost

MD AltSMD, ε = 3.0 AltSMD, ε = 3.5 AltSMD, ε = 4.0

498823 83758 81469 80077

1.0 0.168 0.163 0.161

AltSMD method can reduce the error in SMD simulation when dealing with local atom disorder. Numerical examples of nano-indentation and nanowire tension validate AltSMD method. The results show that AltSMD method can capture the plastic yield point or the broken point correctly and solve the entire problem efficiently. Efficiency comes from the fact that SMD computation can use much larger time step size, and MD relaxation ensures the accuracy. A criterion based on the difference of MD forces and SMD forces is proposed to further improve the efficiency. When the deformation is in the elastic range, the criterion value is small and the computation is mainly executed with

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

293

SMD algorithm. Much more MD relaxation steps can be invoked automatically by the procedure if the local atom disorder appears and the criterion value sharply increases. AltSMD method can be understood as a generalized prediction–correction algorithm with SMD as the predictor and MD as the corrector. In addition to satisfactory accuracy and high efficiency, the alternating process between SMD computation and MD relaxation is natural and straightforward in implementation owing to the similarity between MD and SMD methods. Future work includes efficiently and accurately analyzing problems under higher temperature and constructing an adaptive strategy to determine criterion threshold value. The strategy [45] to adaptively determine background mesh size can also be introduced to further improve AltSMD method. AltSMD method is more suitable for problems with prescribed loading, such as indentation problems and tensile problems. Modifications to MD relaxation process are needed to extend AltSMD method into problems of pure wave propagation. Adaptive coupling of MD and SMD in both spatial and time domains will be another focus in our future work. AltSMD method and its spatial coupling with MD method may be promising in analyzing nano-mechanics problems with defects. Acknowledgments The support from National Natural Science Foundation of China (Grant No. 11102097) and the Beijing Higher Education Young Elite Teacher Project (No. YETP0111) is gratefully acknowledged. References [1] J.Q. Broughton, F.F. Abraham, N. Bernstein, et al., Concurrent coupling of length scales: Methodology and application, Phys. Rev. B 60 (4) (1999) 2391–2403. [2] Y.G. Zheng, Y.T. Zhao, H.F. Ye, H.W. Zhang, Size-dependent elastic moduli and vibrational properties of fivefold twinned copper nanowires, Nanotechnology 25 (2014) 315701. [3] Y.G. Zheng, Y.T. Zhao, H.F. Ye, H.W. Zhang, Y.F. Fu, Atomistic simulation of torsional vibration and plastic deformation of five-fold twinned copper nanowires, Int. J. Comput. Methods 11 (2014) 1344010. [4] Y.L. Chen, B. Liu, X.Q. He, Y. Huang, K.C. Huang, Failure analysis and the optimal toughness design of carbon nanotube-reinforced composites, Compos. Sci. Technol. 70 (2010) 1360–1367. [5] Y. Liu, H.K. Wang, X. Zhang, A multiscale framework for high-velocity impact process with combined material point method and molecular dynamics, Int. J. Mech. Mater. Des. 9 (2013) 127–139. [6] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Crack propagation in bcc crystals studied with a combined finite-element and atomistic model, Phil. Mag. A 64 (1991) 851–878. [7] H.L. Tan, W. Yang, Atomistic/continuum simulation of interfacial fracture, part II: Atomistic/dislocation/continuum simulation, Acta Mech. Sin. 10 (3) (1994) 237–249. [8] L.E. Shilkrot, R.E. Miller, W.A. Curtin, Coupled atomistic and discrete dislocation plasticity, Phys. Rev. Lett. 89 (2002) 025501. [9] W. E, B. Enquist, The heterogeneous multiscale methods, Commun. Math. Sci. 1 (2003) 87–132. [10] V.B. Shenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics the quasicontinuum method, J. Mech. Phys. Solids 47 (1999) 611–642. [11] G.J. Wagner, W.K. Liu, Coupling of atomistic and continuum simulations using a bridging scale decomposition, J. Comput. Phys. 190 (2003) 249–274. [12] D. Qian, G.J. Wagner, W.K. Liu, A multiscale projection method for the analysis of carbon nanotubes, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1603–1632. [13] S.P. Xiao, T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, Comput. Methods Appl. Mech. Eng. 193 (2004) 1645–1669. [14] T. Shimokawa, J. Mortensen, J. Schiotz, K.W. Jacobsen, Matching conditions in the quasicontinuum method: Removal of the error introduced at the interface between the coarse-grained and fully atomistic region, Phys. Rev. B 69 (2004) 214104. [15] J. Fish, M.A. Nuggehally, M.S. Shephard, C.R. Picu, S. Badia, M.L. Parks, M. Gunzburger, Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4548–4560. [16] R.E. Rudd, J.Q. Broughton, Coarse-grained molecular dynamics and the atomic limit of finite elements, Phys. Rev. B 58 (1998) R5893. [17] B. Eidel, A. Stukowski, A variational formulation of the quasicontinuum method based on energy sampling in clusters, J. Mech. Phys. Solids 57 (2009) 87–108. [18] B. Liu, Y. Huang, H. Jiang, S. Qu, K. Hwang, The atomic-scale finite element method, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1849–1864. [19] E. Biyikli, Q. Yang, A.C. To, Multiresolution molecular mechanics: Dynamics, Comput. Methods Appl. Mech. Engrg. 274 (2014) 42–55. [20] L.M. Dupuy, E.B. Tadmor, R.E. Miller, R. Phillips, Finite-temperature quasicontinuum: Molecular dynamics without all the atoms, Phys. Rev. Lett. 95 (2005) 060202. [21] A.C. To, W.K. Liu, A. Kopacz, A finite temperature continuum thoery based on interatomic potential in crystalline solids, Comput. Mech. 42 (2008) 531–541.

294

N. He et al. / Comput. Methods Appl. Mech. Engrg. 296 (2015) 273–294

[22] X. Liu, S. Li, Nonequilibrium multiscale computational model, J. Chem. Phys. 126 (2007) 124105. [23] L. Liu, S. Li, A finite temperature multiscale interphase zone model and simulation of fracture, J. Eng. Mater. Tech. 134 (2012) 031014. [24] J. Fish, W. Chen, R. Li, Generalized mathematical homogenization of atomistic media at finite temperatures in three dimensions, Comput. Methods Appl. Mech. Engrg. 196 (2007) 908–922. [25] M.Z. Xiang, J.Z. Cui, B.W. Li, X. Tian, Atom-continuum coupled model for thermo-mechanical behavior of materials in micro-nano scales, Sci. China: Phy, Mech, Astronom. 55 (2012) 1125–1137. [26] X. Li, J.Z. Yang, W. E, A multiscale coupling method for the modeling of dynamics of solids with application to brittle cracks, J. Comput. Phys. 229 (2010) 3970–3987. [27] G. Anciaux, S.B. Ramisetti, J.F. Molinari, A finite temperature bridging domain method for md-fe coupling with application to a contact problem, Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 204–212. [28] N. Mathew, R.C. Picu, M. Bloomfield, Concurrent coupling of atomistic and continuum models at finite temperature, Comput. Methods Appl. Mech. Engrg. 200 (2011) 765–773. [29] S.B. Ramisetti, G. Anciaux, J.F. Molinari, Spatial filters for bridging molecular dynamics with finite elements at finite temperature, Comput. Methods Appl. Mech. Engrg. 253 (2013) 28–38. [30] S. Reich, Smoothed dynamics of highly oscillatory hamiltonian systems, Physica D 89 (1995) 28–42. [31] B.K. Dey, H. Rabitz, A. Askar, Optimal reduced dimensional representation of classical molecular dynamics, J. Chem. Phys. 119 (2003) 5379–5388. [32] J.A. Izaguirre, Q. Ma, T. Matthey, J. Willcock, T. Slabach, B. Moore, G. Viamontes, Overcoming instabilities in verlet-i/r-respa with the mollified impulse method, Lect. Notes Comput. Sci. Eng. 24 (2002) 146–174. [33] X. Li, W. Yang, Multiple time step molecular dynamics simulation for interaction between dislocations and grain boundaries, Acta Mech. Sin. 21 (2005) 371–379. [34] Q. Pu, Y. Leng, X. Zhao, P.T. Cummings, Molecular simulation studies on the elongation of gold nanowires in benzenedithiol, J. Chem. Phys. C 114 (2010) 10365–10372. [35] D. Qian, S. Chirputkar, Bridging scale simulation of lattice fracture using enriched space–time finite element method, Int. J. Numer. Methods Eng. 97 (2014) 819–850. [36] Y. Liu, X. Zhang, K.Y. Sze, M. Wang, Smoothed molecular dynamics for large step time integration, CMES Comput. Model. Eng. Sci. 20 (3) (2007) 177–192. [37] S. Ma, X. Zhang, X. Qiu, Comparison study of mpm and sph in modeling hypervelocity impact problems, Int. J. Impact Eng. 36 (2009) 272–282. [38] L.M. Shen, A rate-dependent damage/decohesion model for simulating glass fragmentation under impact using the material point method, CMES Comput. Model. Eng. Sci. 49 (2009) 23–45. [39] W.W. Gong, Y. Liu, X. Zhang, H. Ma, Numerical investigation on dynamical response of aluminum foam subject to hypervelocity impact with material point method, CMES Comput. Model. Eng. Sci. 83 (5) (2012) 527–545. [40] J.G. Li, Y. Hamamoto, Y. Liu, X. Zhang, Sloshing impact simulation with material point method and its experimental validations, Comput. & Fluids 103 (2014) 86–99. [41] Z. Guo, W. Yang, MPM/MD handshaking method for multiscale simulation and its application to high energy cluster impacts, Int. J. Mech. Sci. 48 (2) (2006) 145–159. [42] H. Lu, N.P. Daphalapurkar, B. Wang, S. Roy, R. Komanduri, Multiscale simulation from atomisitic to continuum — coupling molecular dynamics (MD) with the material point method (MPM), Phil. Mag. 86 (2006) 2971–2994. [43] L. Shen, Z. Chen, A multi-scale simulation of tungsten film delamination from silicon substrate, Int. J. Solids Struct. 42 (2005) 5036–5056. [44] H.K. Wang, X. Zhang, Y. Liu, Parallel smoothed molecular dynamics method and its coupling with molecular dynamics, Chin. J. Comput. Phys. 25 (2008) 718–724. [45] H.K. Wang, X. Zhang, X.M. Qiu, Adaptive smoothed molecular dynamics for multiscale modeling, Comput. Mater. Sci. 46 (2009) 713–715. [46] D.C. Rapaport, The Art of Molecular Dynamics Simulation, 2nd ed., Cambridge University Press, 2004. [47] S. Foiles, M. Baskes, M. Daw, Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys, Phys. Rev. B 33 (12) (1986) 7983–7991. [48] X. Zhang, Y.P. Lian, Y. Liu, X. Zhou, Material Point Method, Tsinghua University Press, 2013. [49] S.J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19.