Molecular dynamics simulations of homogeneous solids using multi-layered structures

July 2008 EPL, 83 (2008) 10003 www.epljournal.org doi: 10.1209/0295-5075/83/10003 Molecular dynamics simulations of homogeneous solids using multi-layered structures F. A. Gilabert(a) and A. Castellanos Department of Electronics and Electromagnetism, Faculty of Physics, University of Seville Avda. Reina Mercedes s/n, 41012, Seville, Spain, EU received 19 February 2008; accepted in final form 16 May 2008 published online 11 June 2008 PACS 02.70.Ns – Molecular dynamics and particle methods PACS 34.20.Cf – Interatomic potentials and forces PACS 62.20.-x – Mechanical properties of solids Abstract – The main goal of this work is to model a homogeneous computer material with welldefined mechanical properties. To carry out the model material, an internal structure arranged in layers with different atom sizes is implemented using a simple interatomic law of Lennard- Jones type (LJ). We show that imposing an appropriate scaling law between the interatomic potentials from different layers, we obtain the same mechanical properties as if the material was homogeneous. Employing this scheme, given a fixed space volume to be occupied by the solid, this structural arrangement allows to decrease drastically (∼30–80%) the required number of atoms as compared with the case of a homogeneous solid, decreasing the computational effort and speeding up calculations. In that respect, this procedure is an analogous to mesh refinements methodologies usually applied in the continuum approaches. Copyright c EPLA, 2008 Introduction. – Computational designs of solids and soft matter built from simplified microscopic laws of interaction have become an extremely useful tool to model materials with prescribed macroscopic properties. Many of them are devoted to understand phenomena such as plasticity, phase transitions, creeps, large deformations, fracture propagation and so on. Among the extensive literature on this area, the research made in [1,2] is worth mentioning. Over the last years, all these numerical designs are commonly used within the framework of statistical mechanics, since a direct comparison with the theoretical models of liquids and solids is possible [3,4]. Nevertheless, to get a reliable statistical representation, all these realizations still demand high-performance calculations in view of the large number of particles required. For such reason, hybrid models are now very popular in material sciences thanks to the enormous quantity of saved calculations. This enables researchers to introduce much more detailed physical aspects in the model, and therefore it is no longer necessary to sacrifice the physical description by virtue of getting computational improvements. In a hybrid numerical design, a continuum description of the material coexists together with an atomistic scale level to represent the interest zone of the problem [5]. In few words, we may say that when the zone described by MD finishes, (a)E-mail: [email protected] it is immediately continued by a FEM mesh coupled in the interface. This mixed strategy is now commonly used in computational tests of nano-indentation performed on models of metallic materials (among many others, see for example [6]), where the region suffering large deformations is described by the atomic representation, and this is in turn sustained by the continuum region where the deformations are small and predominantly elastic. In this study, we present an alternative scenario, without abandoning the discrete nature of the internal representation of a solid, whose origins, as far as we know, started back with the bidimensional research made in [7,8]. In our deterministic description, corresponding to the mesoscopic level, we assume that particles can be considered as “super-atoms”, which are fully characterized by their positions, velocities, interparticle distances and forces. Using this description, we design a macroscopic solid internally structured in layers. Each layer is an amorphous assembly of super-atoms that possess a conveniently selected “size” or interaction range obeying a simple scaling law. This generates a solid bulk locally rearranged by amorphous atomic layers with different thicknesses. The interaction between atoms has been modeled by means of LJ potential, which has been used to model fluids (see pioneer research in refs. [9,10]) and solids such as metallic glasses, alloys and van der Waals solids [11–13]. This potential is simple, well known 10003-p1 F. A. Gilabert and A. Castellanos and widely used to study the microscopic to macroscopic crossover in condensed matter physics. This paper is organized as follows: Firstly, we present some basic preliminary definitions, nomenclature and the fundamental physical meaning of the variables employed in this study. Subsequently, the motivation of this alternative description is commented, together with our proposal for the interaction and scaling laws. Later on, we shall describe the central aspects of our computational implementation, that is, a detailed description of the layering procedure as well as the advantages obtained from it. Finally, the results obtained from the mechanical tests by means of uniaxial and isotropic compressions are exposed, and compared with those results obtained from the “classical” or homogenous solid bulk. Preliminary issues. – Let us consider a pair of atoms separated by a separation distance rand subjected to the action of a central potential u(r) and interaction force f(r)=u(r). Note that the following definitions are independent from the potential form except for the non-directional and pairs-additivity features. Let us use designations σ,r0,andr1for the distances at which the potential and its first and second derivatives become zero, i.e.:u(σ)=0,u(r0)=0andu (r1) = 0. For all usual potentials, σ<r 0<r 1. The physical meaning of these distances are: σ—hard sphere radius, r0—equilibrium distance and r1—break distance. We shall not use the quantity σ. The three most important dimensional characteristics of the interaction are energy, strength and stiffness of the bond, given by D=−u(r0), fmax = u(r1)andC=u(r0), respectively. All these quantities are positive. In this study we use Dand fmax as the units of energy and force, respectively. For the pressure we define the magnitude σ0=fmax/r2 0, whose meaning is the stress needed to break a single interatomic bond (changing the sign to negative). The unit of time is the period of oscillation of a single atom of mass mgiven by T0=2πm/C. Finally, the critical dissipation coefficient given by γc=√mC represents the quantity of mass that can be slowed down per unit time. The critical value corresponds to the viscosity of an effective medium where a particle is immersed and whose value prevents the particle to complete an oscillation due to the action of the potential field. In that respect, we have considered an athermal problem. Therefore the simplest way to represent a global dissipative mechanism is to fix a constant dissipation coefficient, in our case γ=0.03γc. We are aware of this deliberate simplification, since this choice does not allow us to control the thermal properties of the system. To control the temperature of the material an additional equation for γshould be included, evolving in time and readjusting the energy of the system [14]. Nevertheless, in this work we are interested in calculating the mechanical properties under very low temperatures. To perform the mechanical tests, the barostatic procedure devised by [15] has been implemented. It consists of a mathematical procedure that enables us to apply any kind of external stressed/compressed state on the solid, so that the internal stress of the system, given by the expression of the Cauchy stress tensor, σαβ =1 V i<j ∂u(rij) ∂rα ij rβ ij,(1) can match the externally applied pressure. In (1), the central character of the potential has been assumed, rα ij denotes the α-component of the interatomic distance between atoms iand j,andVis the volume occupied by the solid. Pressure changes are accomplished by changing the coordinates of the particles and the size of the computational cell under periodic boundary conditions. Using this procedure, we have carried out uniaxial and isotropic tests. From these tests, the elastic moduli can be accurately calculated using the internal state of stresses (with eq. (1)) and deformations (relative changes of the periodic cell given by αβ ), i.e. E=2σ2 xx +σxx(σyy +σzz)−(σyy +σzz)2 xx(2σxx +σyy +σzz)−(yy +zz)(σyy +σzz),(2) ν=σxx(yy +zz)−xx(σyy +σzz) (yy +zz)(σyy +σzz)−xx(2σxx +σyy +σzz),(3) B=1 3 σxx +σyy +σzz xx +yy +zz ,(4) where E,νand Bare the Young modulus, the Poisson coefficient and the bulk modulus, respectively. Formulae (2)–(4) can be obtained assuming an isotropic and homogeneous material (or a material with cubic symmetry) together with the constitutive equations of elasticity [16]. Motivation, interaction and scaling laws. – The motivation is to create a computer solid bulk material without compromising neither a good knowledge of their mechanical properties nor the computational effort associated to an elevated number of particles. This improvement is aimed at designing a computer bulk material: i) with well-defined mechanical behavior, ii) with simple but well-determined physics, iii) with an affordable number of particles and, iv) using an efficient and simple algorithm. We combine two ideas: first, an inhomogeneous internal structure and secondly, an appropriate scaling law of the interatomic potentials. We chose to build heterogeneous arrangements of atoms with ordering by layers and with different atom sizes. This structural representation of the material is specially suitable to perform uniaxial tests. Computational efficiency is notoriously improved, since this ordering allows us to classify atoms using techniques inspired by parallel computing using particle systems [17]. The simple idea behind our proposal is that using different atom sizes in a given volume, the number of atoms per volume decreases, and therefore the number of pair interactions decreases as well. The result is an obvious speed-up of the computations: the final purpose of this research is to obtain the same mechanical behavior as in the 10003-p2 MD simulations using multi-layered structures Fig. 1: (Colour on-line) Left: detail of the building material in layers. Black boxes are virtually bounding the region of the sample corresponding with a layer, and the red ones are enclosing the interfaces between a couple of layers. Right: snapshot of the interface between two layers with different atom sizes. An appropriate cut-off distance must be set. For an efficient calculation (see text), this cut-off will depend on the scaling factor used in the interatomic potential for both adjacent layers. homogeneous case but reducing the computing time significatively. A problem for which this methodology will be particulary suitable is the contact between a deformable sphere and a rigid plane. As described in [18], the sphere can be obtained from the cut of a block made of this heterogeneous material. Since the most important events will occur in the contact zone, this contact is more precisely described by the layer that contains the smallest atoms, whilst the rest of the regions far from the contact zone are represented by larger atom sizes. The idea is to create a block of the bulk material arranged by layers, as it is shown in fig. 1, where the layering process is indicated by virtual boxes. Each layer has a different atom size. When we say “atom sizes” we are refereing to equilibrium distances inside the amorphous lattice. In our simple election, the relation between atom sizes belonging to a couple of adjacent layers scales as rl+1 =kr l(l⩾0), where rlis the atom size of the layer l. The first layer of atoms contains the smallest atoms with size r0defined previously. Thus, the atom size in the n-th layer is rn=knr0. Taking into account this scaling, the potential energy between atoms from two contiguous layers land mis ulm(rij)=Dlm rlm rij 12 −2rlm rij 6 ,(5) where rij is the distance between atoms and rlm is the equilibrium distance equal to (klr0+kmr0)/2. This is a long-range potential but in practice one usually defines a cut-off distance from which the potential action becomes negligible. In the case of a homogeneous structure, we set this cut-off distance to rcut-off =2.1r0, that corresponds to a typical distance beyond which forces between pairs are essentially zero. For the scaled potential written in eq. (5), we have established an analogous cut-off distance given by rcut-off,lm =2.1rlm. The elastic moduli of a solid composed by two layers indexed by letters land m, respectively should be invariable if the potential parameters satisfy the following scaling formulae rlm =γlmr0,D lm =γ3 lmD, γlm =1 2(kl+km).(6) To obtain these relations, we have reasoned as follows. Let us suppose a hypothetic elastic modulus of our homogeneous material denoted by M. This modulus has a dependency to the potential parameters Dand r0given by M=λD/r3 0. The quantity λis a number depending i) on the spatial distribution of the lattice nodes (i.e.,the internal structure: amorphous, FCC, BCC, or whatever) and ii) on the interaction range between atoms (i.e., first, second or farther neighbor shells). For a two-layered computer material, the same elastic modulus should have a dependency given by M=λDlm/r3 lm. Then, we admit the following assumption: both types of solids, homogeneous and two-layered, will have the same effective internal structure, i.e.,λ=λ. After a straightforward operation, we obtain the relations written in (6). Therefore, if this assumption is true, a heterogeneous material composed by more than 2 layers will keep their elastic moduli equal to the homogenous ones if the scaling given by (6) is fulfilled between each interface formed by two contiguous layers. General aspects of the implementation. – An efficient MD realization requires, at least, the following features: i) to admit a suitable cut-off distance to avoid the unnecessary counting of particles that do not contribute, ii) an efficient algorithm for searching and sorting the nearest neighbors and iii) a convenient division of the problem using techniques based on the domain decomposition procedure [17]. In a polydispersive size system, as the one we are dealing with here, one can tackle the interparticle interaction in two ways: 1) to consider a mean cut-off distance for all layers using it for any pair with non-equal atom sizes (see the right picture in fig. 1), or 2) to take advantage of the layered structure to classify atoms according to their sizes and recalculating the appropriate cut-off distance as a function of the scaling factor. The first strategy is very easy to implement but produces an inefficient handling of interacting pairs and an inaccurate calculation of forces and internal stresses. In our situation, to compute the total force acting on one atom we need to consider the adjacent atoms possessing different sizes. To consider an average cut-off distance may lead to miss atoms beyond this range and whose contributions to the interaction force are not negligible. On the other hand and at the same time, since each atom has a list of potential neighbors to interact with, to assume this average cut-off distance may lead to include atoms in the list whose contributions are negligible. This implies additional and expensive operations of insertion, calculation and extraction of memory allocations. In that respect, the second strategy, although it is a bit more elaborated to implement, leads to a more efficient 10003-p3 F. A. Gilabert and A. Castellanos Fig. 2: (Colour on-line) Snapshot showing how a continuum bulk material is obtained once a layered sample with nl=5, N= 21311, k=4/3 reaches the equilibrium. In this case, the periodic boundary conditions are only imposed in the directions of the plane perpendicular to the column axis of the sample. handling of the interacting pairs. This reduces the truly necessary interacting number of atoms compared with the previous strategy and with the homogeneous case. Also, the interatomic forces and the stress tensor components are properly calculated, since any interacting pair is not missed. The problem is straightforwardly divided into two stages: first, we perform the interactions between atoms with the same size, i.e. inside a layer, which is computationally enclosed within a virtual bounding box, as the black ones depicted in the left picture of fig. 1, and secondly, we perform the interactions in the zone created by the interface between two adjacent layers, as the red virtual box shown in the same picture. This simple classification of interactions saves many unnecessary memory accesses (atoms inside virtual boxes are located in contiguous blocks of memory) and helps to control all forces calculated on each part of the material. Layering details. – This procedure has been implemented in a C++ code as follows. An amorphous and mechanically equilibrated block of monosized atoms is created at low temperature (<0.1 K). This block is replicated nl−1 times, where nlis the number of layers. Each copy is placed on the previous block, but scaling the atom size according to rl+1 =kr l,i.e., the copied block is “inflated”. Due to the expansion of the scaled block, the boundary box has to be readjusted to keep the same cross-section as the original block. Thus, the atoms falling outside the readjusted box must be removed. In each new layer the mass density remains approximately constant Fig. 3: Atom number as a function of the number of layers in a heterogeneous-layered system in comparison with a homogeneous system (solid line). while the numerical density decreases drastically. Using this procedure, it is easy to calculate the number of atoms in the layer l, given by Nl=N0/k2l, where N0is the atom number in the original block. Next, the program detects nlbounding boxes corresponding to the layers containing monosized atoms and nl−1 bounding boxes corresponding to interfaces containing two-sized atoms (with a negligible numerical cost). Then, the key leading to an efficient performance is that the interaction treatment in each layer becomes identical as in a homogeneous case, i.e., Dlm,rlm with l=m, except for the interface where l=m. The result is a complexity linear with the total particle number, therefore the computation time of this scheme is reduced according to this number. Once a configuration is obtained, the last step consists of getting a continuum bulk. Layers must stick to each other, as fig. 2 shows, whilst forces, stresses and energy are monitored. After the equilibrium, the new multi-layer system is achieved with zero stresses and low temperature. In fig. 3 we show the advantage of using this methodology for two growth rates k. These results can be easily checked, since the total number of atoms in a configuration, as fig. 2 shows, can be accurately calculated just making the following sum: N= nl−1  l=0 N0 k2l=N0 1−k−2nl 1−k−2.(7) This expression is nothing more than the predicted atom number in a multilayered specimen, corresponding to a simple geometrical series with rate kand nlterms. In fig. 3 we have represented with dashed lines those values according to eq. (7), while symbols are the results obtained after the construction of the solid for different layer numbers. The solid line marks the limit established by homogeneous-like constructions. Figure 4 shows the huge advantages obtained from the layering in terms of the saved number of atoms in a given volume: the number of atoms needed to perform the model is drastically reduced, between 30%–80%, depending on nland k. Elastic coefficients and mechanical response. – A variety of specimens has been created, varying the atom number (n), the layer number (nl) and the scaling size 10003-p4 MD simulations using multi-layered structures Table 1: Young modulus E, Poisson coefficient νand bulk modulus Bcalculated for different configurations of the computer bulk layered material. Their expressions are given by formulae (2)–(4). nlkn E/σ 0νB/σ 0 1 1.0 1000 24.5±0.80.370 ±0.131.5±0.2 1 1.0 8000 22.5±0.80.377 ±0.004 31.05 ±0.07 1 1.0 10648 22.1±0.40.379 ±0.002 31.2±0.02 2 1.40 1626 22.55 ±0.07 0.379 ±0.001 31.8±0.3 2 1.50 4375 22.45 ±0.07 0.3849 ±0.001 32.5±0.3 2 1.33 4737 22.48 ±0.02 0.380 ±0.002 32.4±0.4 2 1.50 5113 23.105 ±0.014 0.3765 ±0.0002 31.32 ±0.06 3 1.40 1432 24.4±0.70.372 ±0.008 32.2±0.5 3 1.33 2079 24.3±0.40.371 ±0.005 31.8±0.3 3 1.33 7924 23.5±0.10.375 ±0.001 31.78 ±0.12 4 1.33 8728 23.0±0.40.379 ±0.003 31.9±0.2 5 1.33 21311 24.0±0.50.373 ±0.005 31.4±0.3 Fig. 4: Improvements obtained using the layered model. The reduction in the number of particles is plotted as a function of thenumberoflayers,nl. factor (k). Table 1 shows the results of mechanical tests for these layered materials. Young modulus and Poisson coefficient are obtained from uniaxial tests, while the bulk modulus was calculated from isotropic compressions. Note that the case nl= 1 corresponds to the homogeneous material. The results from table 1 suggest that the elastic moduli seem to be constant, despite small differences probably caused by a slight size effect of the sample. These values were obtained applying stresses below 10−3σ0, producing deformations much smaller than 0.2%: we made sure that the samples never abandoned the elastic regime. Uniaxial stress was applied in the direction perpendicular to the layering, and the isotropic test in the three spatial directions. In both cases periodic boundary conditions were considered. When an uniaxial compression imposing a very slow strain rate is performed (˙T 0=10 −5), the mechanical response of the computer material is clearly linear, as fig. 5 shows. This diagram illustrates the stressstrain curve where the pressure and strain are calculated according to the original dimensions of the sample and not the instantaneous values. We also checked that the same samples behaved in a similar way if a test of uniaxial tension was carried out. In fig. 5, the continuous line was obtained from averaging over five different specimens. The Fig. 5: Stress-strain diagram for homogeneous and heterogeneous material. The material is deformed imposing a quasistatic–like but continuous cell shrinking on the periodic cell that contains the material. The layered specimen has 3 layers with a growth rate k=4/3. bars represent the standard deviation from the average for the heterogeneous material (similar ones are obtained for the homogeneous case). These heterogeneous specimens were created using an identical procedure but changing the initial distributions of positions and velocities of the atoms during the stage of amorphous preparation. Each specimen consisted of three layers with k=4/3. In order to compare with the homogeneous case, we have also plotted with dashed line the stress-strain curve obtained from compressing a homogeneous computer material occupying an equivalent volume. Both types of structures responded with the same mechanical behavior. Even in the domain of large deformations, average values were always within the numerical error bars. This situation is also present in experimental results of physical materials. For the sake of completeness, we performed (on the heterogeneous material used in fig. 5) uniaxial tests imposing the external pressure by means of a stepwise mechanism using the barostatic procedure described in ref. [15] The leaps of pressure were set to 0.05σ0. Figure 6 shows the resultant behavior. Two cycles consisting of loading and unloading were carried out. In the first cycle the specimen was compressed from 0 up to 0.25σ0(black 10003-p5 F. A. Gilabert and A. Castellanos Fig. 6: (Colour on-line) Stress-strain diagram for the heterogeneous material imposing small leaps of pressure. Each symbol represents a mechanical equilibrated state at low temperature. The material definitely breaks once the pressure exceeds the value 0.5σ0. Comparison with the values obtained from homogeneous structures gave an excellent agreement, proving that effects of the scaling law are exactly the expected ones. circles). Once the equilibrium was reached, the specimen was immediately unloaded up to zero pressure (red circles). The second cycle started reloading the previous final state at zero stress, incrementing step by step up to 0.5σ0(black crosses). Decompression process was from 0.5σ0returning to 0 (red crosses). To finish this test, the last state of the second reload was newly compressed up to 0.65σ0(black squares). In this case we observed that the material broke once the value 0.5σ0was exceeded (deformation jumped from 4% up to 12%). Concerning this breakage, there are two specific quantities that complete the characterization of the mechanical behavior together with the elastic moduli: the onset of plasticity σY, or yield strength (beyond which the material flows plastically) and the maximum stress that the material can sustain σU, or ultimate strength. From calculations averaged from over 20 different samples of homogeneous bulk materials consisting of 104atoms each, we have obtained that these quantities are σY/σ0=0.29 ±0.12 and σu/σ0=0.51 ±0.08. In figs. 5 and 6 it can be seen that the heterogenous material behaves identically within this range of values. In fact we obtained similar values for these two quantities: σY/σ0=0.3±0.1andσu/σ0=0.50 ±0.1. Conclusion. – Within the framework of the MD approach, a simple computational procedure based on the scaling of the interatomic potential was used to model heterogenous amorphous solids that possess the same mechanical properties as if they had an internal homogeneous structure. The layered model was subjected to various mechanical tests. From the tests, Young and bulk moduli, the Poisson coefficient and the stress-strain diagrams (controlling the strain rate or the applied pressure) were obtained. The mechanical response showed by layered structures exhibited the same behavior as if they were homogeneous. The combination between simple scaling laws on the interatomic potential and an appropriate layering makes it possible to reduce up to 80% the numerical density of required atoms to describe the solid. This set of ideas presented in this paper are suggested as the beginning of a methodology similar to the mesh refinements approach in FEM, although applied to discrete simulation using particles. In that sense, an immediate application that we have already started is a detailed study of the elastic-plastic and fully plastic regime in adhesive nanocontacts [18]. This method will allow us to build a much larger spherical body to carry out the contact with a rigid plane. We shall be in a position to obtain the JKR solution, commonly used to describe the contact between a spherical macroscopic soft body and a rigid flat surface, in the elastic regime [19]. ∗∗∗ The authors wish to thank Prof. A. M. Krivtsov for his very useful contacts and discussions. This work has been supported by the Ministerio de Educaci´on y Ciencia of the Spanish Government under project FIS2006-03645 and by the Junta de Andaluc´ıa under project FQM-421. REFERENCES [1] Abraham F. F. and Gao H.,Phys. Rev. Lett.,84 (2000) 3113. [2] Gao H., Huang Y. and Abraham F. F.,J. Mech. Phys. Solids,49 (2001) 2113. [3] Allen M. P. and Tildesley D. J.,Computer Simulation of Liquids (Oxford University Press, Oxford) 2002. [4] Hoover Wm. G.,Computational Statistical Mechanics (Elsevier, Amsterdam) 1991. [5] Rafii-Tabar H., Hua L. and Cross M.,J. Phys.: Condens. Matter,10 (1998) 2375. [6] Shiari B., Miller R. E. and Curtin W. A.,J. Eng. Mater. Technol.,127 (2005) 358. [7] Falk M. L. and Langer J. S.,Phys. Rev. E,57 (1998) 7192. [8] Falk M. L.,Phys. Rev. B,10 (1999) 7062. [9] Rahman A.,Phys. Rev.,136 (1964) A405. [10] Verlet V.,Phys. Rev.,117 (1967) 98. [11] Kobayashi S., Maeda K. and Takeuchi S.,Acta Metall.,28 (1980) 1641. [12] Srolovitz D., Vitek V. and Egami T.,Acta Metall., 31 (1983) 335. [13] Quesnel D. J., Rimai D. S. and DeMejo L. P.,Phys. Rev. B,48 (1993) 6795. [14] Nose S.,J. Chem. Phys.,81 (1984) 511. [15] Parrinello M. and Rahman A.,J. Appl. Phys.,52 (1981) 7182. [16] Landau L. D. and Lifshitz E. M.,Theory of Elasticity (Pergamon Press, Headington Hill Hall, Oxford) 1986. [17] Plimpton S.,J. Comput. Phys.,117 (1995) 1. [18] Gilabert F. A., Quintanilla M. A. S., Castellanos A. and Valverde J. M.,Z. Angew. Math. Mech.,87 (2007) 128. [19] Johnson K. L., Kendall K. and Roberts A. D.,Proc. R. Soc. London, Ser. A,324 (1971) 301. 10003-p6