New Bending Algorithm for Field-Driven Molecular Dynamics

A field-driven bending method is introduced in this paper according to the coordinate transformation between straight and curved coordinates. This novel method can incorporate with the periodic boundary conditions in analysis along axial, bending, and transverse directions. For the case of small bending, the bending strain can be compatible with the beam theory. Consequently, it can be regarded as a generalized SLLOD algorithm. In this work, the bulk copper beam under bending is analyzed first by the novel bending method. The bending stress estimated here is well consistent to the results predicted by the beam theory. Moreover, a hollow nanowire is also analyzed. The zigzag traces of atomic stress and the corresponding 422 common neighbor type can be observed near the inner surface of the hollow nanowire, which values are increased with an increase of time. It can be seen that the novel bending method with periodic boundary condition along axial direction can provide a more physical significance than the traditional method with fixed boundary condition.


Introduction
The nano-scale mechanical properties become important since the size of electrical components is successively reduced for the portable convenience [1,2]. Most of studies focused on mechanical properties related to the tension and compression. The problems of bending are actually met more frequently although it is composed by tension and compression. The bending tests of nanomaterials by using atomic simulation were widely applied. Liu et al. [3] simulated the pure bending of defect-free Al single crystals to investigate dislocation nucleation from free surfaces. They found that dislocation nucleation is not well represented by a critical value of the resolved shear stress but is reasonably well represented by the critical stress-gradient criterion. On the other hand, the size effects were also discussed widely. Miller and Shenoy [4] found that the surface elastic constant is the same order as the bulk elastic constant. The surface effect was also discussed in the bending case.
Unlike the case of tensile or compression tests of nanowire (NW) or nanofilm (NF) where the periodic boundary conditions (PBCs) were applied along the axial direction to remove the size effect, almost all the bending simulations took the ends of nanowire or nanofilm as fixed boundary conditions (FxBCs) [3]. The FxBC is essentially inducing the size effects into the simulated objects. From the viewpoint of thermodynamics, the fixed atoms are viewed as zero velocities, and, thus, zero temperature at the fixed ends. In other words, all thermodynamic variables involving atom velocities are not well defined at the fixed boundary.
For the purpose of the computational efficiency, there are many methods to improve the computational speed. One of the methods is to use the synthetic system instead of real system. For example, Nose-Hoover algorithm [5][6][7], the synthetic thermostat variable generates the NVT ensemble more stably and efficiently than the rescaled velocity method [8], the latter cannot generate the NVT ensemble exactly. Moreover, the synthetic system usually combines the physical response into the equations of motion, thus can prevent the discontinuous trajectory of atoms and save the time to do the local equilibrium.
Non-equilibrium molecular dynamics (NEMD) can be described as two representations [9]. One is the boundarydriven (BD) representation, the other is the field-driven (FD) representation. The FD method is belonging to the synthetic system. The BD method was used to calculate the thermal transport coefficients while the FD method was used to calculate the mechanical ones. One can mathematically transform the non-equilibrium boundary conditions for a thermal transport process into a mechanical field. The two representations of the system are said to be ''congruent''. Almost all FD methods can combine with PBC while BD methods usually combine with FxBCs. In addition, since the non-equilibrium response is reflected in the equations of motion for FD method, it is no need to use the stepwise equilibrium-nonequilibrium cyclic driven that usually used in the BD method. Thus the FD method can be more efficient than BD method.
From the above reasons, a novel-bending algorithm is proposed and investigated in this paper. Based on the coordinate transformation from flat coordinate to curved one, the straight material is transformed to the curved one. The method belongs to FD method, and can be viewed as the generalized SLLOD algorithm [9]. It also removes the fixed atoms generally used at FxBCs so that all thermodynamic variables involving atom velocities can be defined everywhere. The method for the bending algorithm is introduced in ''Methodology'' section. ''Numerical Tests'' section shows some numerical tests for both macroscopic and microscopic systems. Finally, it is concluded in ''Conclusion'' section.

Methodology
The Coordinate Transformation For a tension simulation, one may view the stretch as a coordinate transformation described as The Jacobian and inverse Jacobian are From Fig. 1, the original straight material that has length L in the (x,y,z) coordinate transforms to (x 0 ,y 0 ,z 0 ) coordinate. In the view of (x 0 ,y 0 ,z 0 ) coordinate, the straight material has a length of aL. The factor 1/a of the transformation can be viewed as the factor of stretch. One can set the coordinate of next time (x(t),y(t),z(t)) to be equal to the coordinate (x 0 ,y 0 ,z 0 ), thus the dynamical, uniform, and field-driven stretch can be performed.
From the above idea, one can search a coordinate transformation for the bending purpose. For the simplest bending case, one may consider the curved axis x 0 as a quadratic curve of the formỹ ¼ ax 2 (see Fig. 2). The slope at x 0 can be obtained by where d is a horizontal distance between x and x 0 . In order to satisfy the assumption of ''a plane normal to the axis remains a plane normal to the curved axis after bending'' in the beam theory [10], the y 0 axis is set to be normal to the curved axis x 0 . Thus for an arbitrary point P(x,y) = P(x 0 ,y 0 ), the coordinate transformation can be obtained and given by The distance d can be evaluated by the assumption that the normal at x 0 is orthogonal to the tangent at x 0 , i.e., fððx þ dÞ;aðx þ dÞ 2 ÞÀðx;yÞg Á ð1;2aðx þ dÞÞ ¼0 It is a cubic equation having three roots. The discriminant D can be used to check the roots [11], For D \ 0, there are three different real roots; for D = 0, there are triple or double real roots; and for D [ 0, there is only one real root and two imaginary roots. The geometrical condition requests that the distance d is a real root, thus D must be greater than zero. One can see that if we set y\ 1 2a ; then D must be greater than zero. Thus the condition y\ 1 2a is selected as a limit range of the coordinate transformation.
For a very small deflection, i.e., a 2 ! 0; the Eqs. 3-7 can be reduced as Thus the coordinate transformation and inverse transformation can be obtained, And the corresponding Jacobian and inverse Jacobian are It can be seen that J j j J À1 ¼ 1: After constructing the J and J -1 , the bending simulation can be performed as Fig. 3. The coordinate (x,y,z) is transformed to (x 0 ,y 0 ,z 0 ) with the curvature-related coefficient a. Thus the simulated system is curved in the viewpoint of (x 0 ,y 0 ,z 0 ) coordinate. The (x 0 ,y 0 ,z 0 ) can then be taken as the coordinate at next time step to perform the bending dynamic simulation. Fig. 1 Tension/compression by using coordinate transformation. The coordinate (x,y,z) is transformed to (x 0 ,y 0 ,z 0 ) with scale 1/a along z 0 -direction. Thus the length L becomes aL in the (x 0 ,y 0 ,z 0 ) scale. The (x 0 ,y 0 ,z 0 ) can then be taken as the coordinate at next time step to perform the tension/ compression dynamic simulation Fig. 2 Curved coordinate transformation between coordinates (x,y,z) and (x 0 ,y 0 ,z 0 ). The point P is located at (x,y) in the flat coordinate while at (x 0 ,y 0 ) in the curved one. The y 0 axis is turned with h refers to the vertical direction and perpendicular to x 0 axis. The distance between x and x 0 in the horizontal plane is d Fig. 3 Schematic diagram of bending by using coordinate transformation. The coordinate (x,y,z) is transformed to (x 0 ,y 0 ,z 0 ) with the curvature-related coefficient a. Thus the simulated system is curved in the (x 0 ,y 0 ,z 0 ) scale. The (x 0 ,y 0 ,z 0 ) can then be taken as the coordinate at next time step to perform the bending dynamic simulation The Bending Strain For the curvature, it can be shown that Note that the curvature is independent of coordinates; it is convenient to characterize the bending status. For the displacement u x = x 0 -x = d and u y = y 0 -y = -a(x ? d) 2 , the strain components can be obtained And the volume change can be calculated as where L, b, and c are the length, width, and depth of the beam, respectively. Thus the volume of the beam can be viewed as no change after bending. At y = 0, the axial strain e xx = 0, thus the axis x 0 can be considered as centroid axis or neutral surface. The linear relation between e xx and y is also consistent with the assumption of beam theory [10]. For e xy = 0, it is also conformed to the assumption that plane sections initially normal to the beam axis remain plane and normal to that axis after bending. The transverse strain e yy % 0 also meets the assumption of beam theory. Thus the model can be used to verify the suitability of beam theory in the nanoscale model with slight bending.

The SLOOD Algorithm
Let a(t) be a function of time t, the rate of displacement gradient tensor can be written as During the motion of bending, a particle i can experience a velocity where x i is the position of particle i. Thus the equations of motion for the position of atom i can be written as By using the mean value theorem for vector-valued functions [12], the above equation can be rewritten as where r _ Uðx; tÞ ¼ R 1 0 dkr _ uðkx; tÞ: Assuming the equation of motion for the conjugate momentum induced by the bending can be written as where p i , F i , and m i are the conjugate momentum, force, and mass, respectively, of particle i. Equations 25 and 26 are one example of the general NEMD equations of motion where C i and D i are the phase variables coupling of the field F e (t) to the system. Equations 25 and 26 can be reduced to original SLLOD algorithm if the rate of displacement gradient tensor is independent of coordinate [9]. Thus Eqs. 25 and 26 can be viewed as the generalized SLLOD algorithm. Note that Eqs. 25 and 26 can also be used in the case with large deformation.

Periodic Boundary Conditions
Since the curved coordinate and flat coordinate can be transformed to each other at each time step, one can transform the curved system to rectangular one, and then the PBCs can be applied (see Fig. 4). After the PBCs are applied, the atomic positions are restored to the curved coordinates and do the simulation at next time step. The minimum image criteria (MIC) can also be implemented. As shown in Fig. 4, the curved cell can be transformed to rectangular one by J -1 . After building the image cells, the neighbor list can be built up [8], and then the rectangular cell is restored to the curved cell by J. The forces between any neighbor atoms thus are calculated according to the neighbor list and the curved positions.
In the mathematical description of the primary cell, one can consider the position of atom i that is changed due to the bending essentially, i.e., the second term of Eq. 25. It is convenient to use a shape matrix h(x,t) to describe the primary cell so that x i (t) = s i Áh(x i ,t), where s i is the scaled position of atom i [13]. Thus the position x i is fully changed by h if the scaled position s i is not a function of t, and the second term of Eq. 25 becomes Equation 29 then serves the equation of motion for the primary cell.

Bulk Copper Beam
The bulk bending simulation is first tested for verifying the new bending method. The copper atomic system is constructed as 20a Lattice 9 10a Lattice 9 10a Lattice (equals to 72.38 9 36.19 9 36.19 Å 3 ) with 8000 atoms, where a Lattice = 3.619 Å is the lattice constant for copper. The x-, y-, and z-directions are along to the [100], [110], and [111] orientations, respectively. The PBCs/MIC are imposed along the x-, y-, and z-directions, and the Morse potential is adopted [8]. Verlet neighbor list combined with cell-link method is used [14]. Gesar fifth-order predictorcorrector algorithm [8] with time step 1 fs is also applied. The atomic stress formula for the atomic system is calculated by [15] where r ab i is the atomic stress with component ab for atom i; m i , v a i ; r a ij ¼ r a i À r a j are the mass for atom i, a component velocity for atom i, a component distance between atoms i and j, respectively; and / is the potential energy. The total stress for the whole system is where the total volume V and atomic volume V i are related by The system is equilibrated first in the NrT ensemble with GGMT thermostat method (with five thermostat variables) [16] at 300 K and MTK barostat method [17] at 0 GPa during 0.1 ns. Then the system is bended with bending rate _ j ¼ 2 _ a ¼ 5 Â 10 À8 fs -1 Å -1 , which equals to the bending rate at the beam end of _ y 0 ¼ 2 _ axj x¼36:19 ¼ 1:81 rad ns -1 , and the temperature is controlled at 300 K with GGMT method.
The front and isometric views of the bending system with curvature j = 4.5 9 10 -3 Å -1 are shown in Fig. 5. According to the beam theory [10], the normal stress r xx can be estimated by where E is the Young's modulus. For Cu, E = 117.2 GPa [18], thus the maximum and minimum stresses are r xx ¼ 117:2 Â 4:5 Â 10 À3 Â ðAE18:095Þ ¼ AE9:54 GPa: The simulated values with the bending algorithm are r max = 7.81 GPa and r min = -6.79 GPa which take the average among the top and bottom atomic layers, respectively. The values of atomic simulation are close to the one estimated by continuous mechanical beam theory.
If the crystalline effect is considered in the Young's modulus which is given by [19] 1 x n 2 y þ n 2 y n 2 z þ n 2 z n 2 where C 11 , C 12 , and C 44 are the elastic moduli, (n x ,n y ,n z ) is the crystalline orientation. For Cu, the moduli are C 11 = 168.4 GPa, C 12 = 121.4 GPa, and C 44 = 75.4 GPa, respectively, at 300 K [20]. Thus the Young's modulus along [100] orientation is E = 66.69 GPa, and gives r xx ¼ 66:69 Â 4:5 Â 10 À3 Â ðAE18:095Þ ¼ AE5:43 GPa: The values of the axial stress obtained in this work are distributed just in between. The possible cause is that the lattice constant is stretched/compressed so that the elastic moduli are no longer the same as reference [20]. Thus the axial stresses at top and bottom are not symmetric, and deviate from Eq. 36. Once the model is compatible with macroscopic bulk, it is confident to use the bending method to simulate the microscopic nano-system.

Hollow Copper Nanowire
Nanowires (NWs) exhibit an interesting quantum conductance behavior even at room temperature. Electron transport properties for NWs have been investigated extensively due to their significant importance in a variety of applications [21]. Diao et al. [22] investigated the elastic properties of Au NWs by molecular statics, and found that due to the surface effects, the smaller the cross-sectional area the higher the Young's modulus in the NWs without undergoing the phase transformation. Chen and Chen [23] studied the Au NWs subjected to uniaxial tension at high strain-rate under different temperatures. They found the microstructures of NWs were transformed first from FCC to face-centred-orthorhombic-like crystalline, and then changed to the amorphous state. Moreover, it was predicted that the conductance at high strain-rate deformation may be no longer quantized. Recent research has revealed that geometry, including surface orientation and the hollowness of nanomaterials, can also greatly affect their behavior [24][25][26][27][28].
The works of Jiang and Zheng [27,28] are referred here to compare the effect induced by different boundaries. The system size studied in this paper is same as Zheng's work (outer and inner cross-section parameters are 10a Lattice and 4a lattice , respectively), except the PBC is applied here along the axial direction instead of FxBC used in Zheng's work. Other settings remain the same as previous sub-section. The system is equilibrated in NVT ensemble with GGMT method at 10 K before bending.
The axial stress distribution of hollow NW after bending with j = 2.25 9 10 -3 Å -1 at 10 K is shown in Fig. 6. It is observed that higher tensile stresses are created near the corner and dislocations. Higher stress at corner is mainly induced by the surface tension; while the dislocations are found at (100) surfaces, identical to report in the reference [28].
The technique of common neighbor analysis (CNA) is adopted to analyze the local structure distribution [23,29,30]. The CNA has three indices, j, k, and l, which denote the number of common neighbor (CN) particles, the pair number of CN particles, and the number of CN pairs that makes up a chain, respectively. A pair is constructed by two particles whose distance apart is less than a cutoff radius. The cutoff radius is chosen to be 1.2d NN [23], where d NN is the nearest-neighbor distance. For the perfect FCC structure, the probability of 421 CN type is 100%, while the perfect HCP structure contains 50% 421 CN type and 50% 422 CN type. The structures of 421 and 422 CN types are shown in Fig. 7 for reference.  Fig. 8 show the 422 CN-type distributions. The lower pictures show the axial stress r xx distributions with color range from -10 to 10 GPa. The growing atomic tensile stress can be observed, and the 422 CN-type growths are accompanied with. The stress trace and 422 CN type initially grow along 110 Â Ã direction, and reflected to [110] direction by the inner corner, then form a zigzag trace. Note that the depth of the stress trace and 422 CN type is only about one or two atomic layers from the inner surface.
Since the PBC is applied here, the atoms at the boundary are movable so that the stress trace can grow through the ends. On the contrary, if the FxBC is applied, the atoms at the ends will be fixed, and the interface between movable and fixed atoms will lead to an artificially induced crack, obviously violating the physical phenomenon.

Conclusion
In this study, the synthetic, field-driven bending method is introduced by using the coordinate transformation between straight and curved coordinates. The new method can incorporate with PBCs along axial, bending, and transverse directions. For problem with small bending effect, the bending strains evaluated by this method are well consistent with those predicted by the beam theory. Furthermore, it can be regarded as the generalized SLLOD algorithm. The accuracy and reliability of this novel bending method are verified by two examples, which are the bulk copper beam and the hollow NW under bending, respectively. The bending stress of the bulk copper beam estimated here is quite close to those predicted by the beam theory; while the atomic stress and the corresponding microstructure of 422 CN type near the inner surface of the hollow NW are increased with an increase of time. These results are well consistent with the earlier work. Moreover, the performance of this novel bending method can be significantly enhanced by using PBC along axial direction in the bending model by eliminating the artificial crack which is easily created by using traditional method with FxBC. Fig. 7 The a 421 and b 422 CN types for CNA. Red atoms indicate any two atoms within a neighbor distance and form a pair. Blue atoms (#3-#6) are the common neighbor atoms for red atom pair (#1 and #2). The black line between two atoms indicates these two atoms form a pair