Abstract
This work presents a detailed description of the formulation and implementation of the Atomistic Finite Element Method AFEM, exemplified in the analysis of one and twodimensional atomic domains governed by the Lennard Jones interatomic potential. The methodology to synthesize element stiffness matrices and load vectors, the potential energy modification of the atomistic finite elements (AFE) to account for boundary edge effects, the inclusion of boundary conditions is carefully described. The conceptual relation between the cutoff radius of interatomic potentials and the number of nodes in the AFE is addressed and exemplified for the 1D case. For the 1D case elements with 3, 5 and 7 nodes were addressed. The AFEM has been used to describe the mechanical behavior of onedimensional atomic arrays as well as twodimensional lattices of atoms. The examples also included the analysis of pristine domains, as well as domains with missing atoms, defects, or vacancies. Results are compared with classical molecular dynamic simulations (MD) performed using a commercial package. The results have been very encouraging in terms of accuracy and in the computational effort necessary to execute both methodologies, AFEM and MD. The methodology can be expanded to model any domain described by an interatomic energy potential.
Keywords:
AFEM; Interatomic potentials; Lennard Jones potential; molecular dynamics
1 INTRODUCTION
The last decades have witnessed an enormous increase in the development of nanotechnologies, which have brought a large impact on areas such as engineering, medicine, biology, chemistry, as well as in aerospace, construction and manufacturing industries. The scientific community has been working to analyze and develop mechanical, electronic and optoelectronic devices at nanoscale, which can be applied to many areas which range from sensors engineering to medical therapies and drug delivery systems (The Royal Society and The Royal Academy of Engineering, 2004Nanoscience and nanotechnologies: opportunities and uncertainties (2004). The Royal Society and The Royal Academy of Engineering, London.; Novoselov et al., 2012Novoselov, K.S., Fal′Ko, V.I., Colombo, L., Gellert, P.R., Schwab, M.G., Kim, K. (2012). A roadmap for graphene, Nature 490: 192.; Chen et al., 2013Chen, C., Hone, J. (2013). Graphene Nanoelectromechanical Systems. Proceedings of the IEEE 101(7).). Nanomechanics of materials is a new branch of mechanics which studies the properties and behavior of nanoscale materials and structures in response to applied loads and prescribed boundary conditions. From the dimensional perspective, a structure with at least one dimension smaller than 100 nm is considered to be a nanostructure (The Royal Society and The Royal Academy of Engineering, 2004).
The mechanical behaviour of nanoscale systems can be analyzed by using Abinitio methods (Dirac, 1981Dirac, P.A.M. (1981). The Principles of Quantum Mechanics, Clarendon Press, Oxford.), semiempirical quantum methods (Tadmor and Miller, 2011Tadmor, E.B., Miller, R.E. (2011). Modeling materials: continuum, atomistic, and multiscale techniques. Cambridge University Press, New York.) and also modified classical continuum mechanics. Abinitio calculations are most accurate compared to semiempirical quantum methods, but computationally very expensive and the analyses are still limited to a few hundred or thousand atoms. The molecular dynamics (MD) (Alder et al., 1959Alder, B.J., Wainwright, T.E.J. (1959). Studies in Molecular Dynamics. I. General Method, The Journal of Chemical Physics 31: 459466.) has been one of the most commonly method used to analyze the behaviour of nanomaterials. Although MD methods can be applied to analyze larger systems, than those that can be treated by Abinitio methods, it is still not applicable to systems consisting of many millions of atoms, which are frequently necessary to model practical devices.
For nanostructures greater than a few nanometres, multiscale continuum methods such as GurtinMurdoch theory (GMT) (Gurtin and Murdoch, 1975Gurtin, M., Ian Murdoch, A. (1975). A continuum theory of elastic material surfaces, Archive for Rational Mechanics and Analysis 57: 291323.), that bridges surface and bulk energies, is quite accurate and computationally inexpensive (Dingerville et al., 2005Dingerville, R., Qu, J., Cherakoui, M. (2005). Surface free energy and its effect on the elastic behavior of nanosized particles wires and films, Journal of the Mechanics and Physics of Solids 53: 18271854.; Miller et al., 2000Miller, R.E., (2000). Sizedependent elastic properties of nanosized structural elements, Nanotechnology, 11:139147.; Liu et al., 2010Liu, C., Rajapakse, R.K.N.D. (2010). Continuum models incorporating surface energy for static and dynamic response of nanoscale beams, IEEE Transactions on Nanotechnology 9: 422431.; Sapsathiarn et al., 2012Sapsathiarn, Y., Rajapakse, R.K.N.D. (2012). A model for large deflections of nanobeams and experimental comparison, IEEE Transactions on Nanotechnology 11, 247254.). However, GMT and other continuum methods are not strictly applicable to model discrete atomic structures. The computational cost in modelling nanomaterials is a fundamental challenge. In order to overcome this challenge, the Atomicscale Finite Element Method (AFEM) was developed and applied for multiscale analysis of carbon nanotubes (CNTs) by Liu and his coauthors (Liu et al., 2004; Liu et al., 2005Liu, B., Jiang, H., Huang, Y., Qu, S., Yu, M.F. (2005). Atomicscale finite element method in multiscale computation with applications to carbon nanotubes, Physical Review B, 72: 10980121.). Presently, the AFEM can model and analyze the mechanical behavior of many materials at the nanoscale (Malakouti, 2006Malakouti, M., Montazeri, A. (2016). Nanomechanics analysis of perfect and defected graphene sheets via a novel atomicscale finite element method, Superlattices and Microstructures 94: 1 12.; Tao et al., 2016Tao, J., Sun, Y., (2016). The Elastic Property of Bulk Silicon Nanomaterials through an Atomic Simulation Method. Journal of Nanomaterials 2016: 6.; Damasceno et al., 2015Damasceno, D.A., Rajapakse, R.K.N.D., Mesquita, E. (2015). Mechanical behavior of atoms using atomicscale finite element method. Proceedings of the XXXVI IberoLatin American Congress on Computational Methods in Engineering.).
The AFEM is based on an energy approach. It requires an interatomic energy potential, describing local or nonlocal bonding forces of an atom interacting with a chosen set of other surrounding atoms. The choice of the interacting surrounding atoms and their array will, ultimately, lead to the formulation of the specific Atomic Finite Element. In principle, the method can be applied to all atomic systems that may be described by an interatomic energy potential.
The purpose of this work is to introduce and discuss fundamental issues about the formulation and application of the AtomicScale Finite Element Method (AFEM). This methodology will be applied, exemplary, to analyze the mechanical behavior of onedimensional (1D) and twodimensional (2D) atomic structures and lattices. The LennardJones interatomic potential will be used in the formulation and analysis. In particular, the article describes how to relate the number of nodes of a 1D Atomic Finite Element with the cutoff radius of the LennardJones potential. It also addresses, in detail, the modifications at the element level that have to be performed to model bounded atomic domains and to introduce proper boundary conditions. The obtained numerical results are compared with classical MD simulations. Accuracy and computational costs involved in both methodologies are addressed.
2 ATOMIC SCALE FINITE ELEMENT METHOD FORMULATION
The AFEM, as described by Liu et al. (2004Liu, B., Huang, Y., Jiang, H., Qu, S., Hwang, K.C. (2004). The atomicscale finite element method, Computer Methods in Applied Mechanics and Engineering 193: 18491864.), is based on an energy approach. It requires an interatomic energy potential, U, describing local or nonlocal bonding forces of an atom interacting with a chosen set of other surrounding atoms. The interatomic energy potentials are described in terms of the coordinates of the individual atoms, x_{i}, and some constitutive parameters of the considered atomic elements and arrays.
For a system with N atoms, the total interatomic energy considers the contribution of all interacting atoms:
Consider the energy W_{f}, related to the sum of the work of external forces, ${\overline{)F}}_{i}^{ext}$ acting on individual atoms indicated by ‘i':
The total energy in this system with N atoms, E_{tot}, is composed of the interatomic bonding energy and the (eventual) work of external forces.
As usual, the desired equilibrium condition is obtained by determining the stationary value of the total energy in the system with respect to the position x_{i} of the individual atoms:
The total energy E_{tot} can be expanded in a Taylor series around the equilibrium position of the atoms, x^{(0)}:
Defining the displacement u as the difference between the actual, x, and the equilibrium position, x^{(0)}
And substituting the Equation (5) into Equation (4) leads to the AFEM equation system:
Where [K(u)] corresponds to the nonlinear stiffness matrix, {u}, the displacement increment vector and {P} the nonequilibrium load vector, respectively, given by:
And
3 LENNARD JONES POTENTIAL
In this section, a brief review of the LennardJones (LJ) potential proposed by Jones (1924Jones, J.E. (1924). On the determination of molecular fieldsII. From the equation of state of a gas. Proceedings of the Royal Society A 106: 463477.) is introduced. It is used to calculate the potential energy of the interaction between two nonbonding atoms based on their distance of separation, r_{ij}:
Where ${\text{r}}_{\text{ij}}=\sqrt{{\left({\text{x}}_{\text{i}}{\text{ x}}_{\text{j}}\right)}^{2}+{\left({\text{y}}_{\text{i}}{\text{ y}}_{\text{j}}\right)}^{2}}$ is the bond length, and the r_{c} is some prescribed cutoff radius, which switches off the atomic interaction when the bond length exceeds its value. Only two LennardJones parameters are sufficient to describe the intermolecular interactions, σ and ( _{LJ} , as shown in Figure 1. The parameter ( _{LJ} is the depth of the potential well and σ is the interparticle distance at which the potential is zero. The ${\left(\raisebox{1ex}{$\sigma $}\!\left/ \!\raisebox{1ex}{${\text{r}}_{\text{ij}}$}\right.\right)}^{12}$ term represents the strong repulsion and the high energy on the bond, as the distance between the pair of atoms decreases. The ${\left(\raisebox{1ex}{$\sigma $}\!\left/ \!\raisebox{1ex}{${\text{r}}_{\text{ij}}$}\right.\right)}^{6}$ term represents the attractive term, which describes the weak attraction on the bond, as the distance between the pair of atoms increase.
The internal force, F(r_{ij}), between the two atoms is obtained by differentiating the potential with respect to the intermolecular distance, r_{ij}:
Figure 1 shows the LJ potential, U(r), and the interatomic force, F(r), as a function of the distance, r_{ij}. By setting the force F(r) =0, the equilibrium distance between the atoms, r_{eq} is found to be:
If the distance between two atoms is less than r_{eq}, the LennardJones force is repulsive. If the distance is greater than r_{eq}, the force is attractive. The maximal force is obtained at the intermolecular distance,
It should be noted that for intermolecular distances larger than r_{fmax} the interatomic bonding force decreases continuously. When the distance between atoms reaches the equilibrium distance between the atoms, r_{eq}, the bonding force is assumed to vanish.
The cutoff radius r_{c} is an important parameter to be chosen. It determines with how many neighboring atoms, a given atom is interacting. The cutoff radius r_{c} also determines whether the atomic interaction has a local or nonlocal character. In this article the cutoff radius of the interatomic bonding will be related to the number of nodes within a specific Atomic Finite Element. This will be illustrated in the next paragraphs within the context of 1D elements.
Consider an atomic chain with 5 atoms in equilibrium position, as shown in Figure 2a. Every considered atom interacts with one or more of its neighbors, according to the relation between the interatomic equilibrium distance r_{eq} and the chosen cutoff radius r_{c}. Three distinct conditions will be analyzed. First, consider that the cutoff radius obeys the relation (r_{eq} < r_{c} < 2r_{eq}). This means that one atom only interacts with the first nearest neighboring atoms, as shown in Figure 2b. For the second condition, (2 r_{eq} < r_{c} < 3 r_{eq}), every atom interacts with the first and the second nearest neighboring atoms, as shown in Figure 2c. The third condition, (3 r_{eq} < r_{c} < 4 r_{eq}), is illustrated in Figure 2d.
An atomic chain with 5 atoms and distinct relations between the equilibrium distance r_{eq} and the chosen cutoff radius of the LJ potentials, r_{c}.
Two aspects should be stressed. First, when the cutoff radius is smaller than two times the equilibrium distance, r_{c} <2 r_{eq}, the atomic bondings have a local character. For larger cutoff radius, r_{c} >2 r_{eq}, the bonding extends to the second or more distant neighbors and gives a nonlocal character to the interatomic interactions. The second comment is related to the finitude of the atomic chain considered. The atoms 1 and 5, indicated in Figure 2a, represent the physical limit of the considered chain. The lack of atoms to the left of atom number 1 and to the right of atom number 5 impacts directly the energy bonds of the chain. If r_{c} >2 r_{eq}, nonlocal bonds exist and not only the bond energy of the most external atoms (1 and 5, at the present case) if affected. This must be carefully considered when modelling a bounded atomic domain with the AFEM.
4 ATOMIC FINTE ELEMENTS
In this section the formulation for 1D and 2D Atomic Finite Elements are, exemplarily, presented for the case the LennardJones interatomic potential. The idea is to focus on the fundamental issues of the AFEM formulation and not on material properties of the domains to be analyzed.
4.1 The 1D Atomic Element
The notion of the atomic finite element will be introduced based on three examples. Consider Figure 3, showing an atomic finite element with three nodes, designated as i1, i and i+1. The definition of the atomic finite element and all the derived calculations, such as internal forces and stiffness matrix, are related to the central node, in this case, node i. It is assumed that the relation between the equilibrium distance and the potential cutoff radius satisfies (r_{eq} < r_{c} < 2 r_{eq}). The arrows in Figure 3 a shows that the central atom i, is only connected to its first neighbor atom to the left (i1) and to the right (i+1). The total bonding energy for this 3 node atomic finite element i is composed of the bonding (i1, i) and (i, i+1):
The expressions for element stiffness matrix, ${\text{K}}_{\text{i}}^{\text{3e}}$ and nonequilibrium force vector, ${\text{P}}_{\text{i}}^{\text{3e}}$ and for this atomic three node finite element are obtained by substituting the expression for the bonding energy (14) and the external forces in equations (8) and (9), resulting delivering, respectively:
The LennardJones potential is a pairwise potential, describing the interaction between two atoms, that is the interaction of atom (i) over atom (i+1) and the reciprocal interaction of atom (i+1) on atom (i). In the present formulation the atomic element is represented by node (i) and only the interaction of node (i) with (i1) and (i) with (i+1) will be considered. This is implied by the directions of the arrows in Figure 3. The reciprocal interaction, of the neighboring atomic element centered in the node (i1) on the node (i), of the present element, will be added later at the mesh assembly procedure. So only half of the bonding energy between (i) and (i1), and (i) and (i+1) has to be considered. This is expressed by the 1/2 factor in the element stiffness matrix.
Figure 4a shows a 5node atomic finite element. The central node is denoted by (i) and it interacts with the two nearest nodes, (i1) and (i2) to the left and (i+1) and (i+2) to the right. For this element the condition (2r_{eq} < r_{c} < 3r_{eq}) is considered. Analogously, for the 7node element, shown in Figure 4b, the relation between equilibrium distance and cutoff radius is (3r_{eq} < r_{c} < 4r_{eq}).
The total bonding energy of these elements are, respectively, given by:
4.2 Atomic Element in Two Dimensions
In this section twodimensional atomic finite elements are presented. The starting point is the 7node 2D element shown in Figure 5. As it will be shown, it may be used to model the geometry of hexagonal Bravais type lattices (Liu et al., 2006Liu, W.K., Karpov, E.G., Park, H.S. (2006). Nano Mechanics and Materials: Theory, Multiscale Methods and Applications, John Wiley & Sons.). The central atom, (i), has connection with all of surrounding atoms, (i+1), (i+2), (i+3), (i+4), (i+5) and (i+6). The total energy for this basic 2D 7node atomic finite element is the sum of all pairwise interacting nodes:
The element stiffness matrix and the element nonequilibrium force vector for the atomic finite element with seven nodes are given by, respectively:
Figures 6a and 6b illustrate two other 2D atomic finite elements. These atomic finite elements are defined based on the types of Bravais lattices (Liu et al., 2006Liu, W.K., Karpov, E.G., Park, H.S. (2006). Nano Mechanics and Materials: Theory, Multiscale Methods and Applications, John Wiley & Sons.). shown in Figure 7a and 7b respectively. The extension of the AFEM to these elements is straight forward.
5 AFEM IMPLEMENTATION
The assemblage procedure in one and two dimensions, the concept of modified atomic finite element and the application of boundary conditions will be discussed in this section. Figure 8 illustrates the AFEM assemblage procedure for a 1D atomic chain containing 4 atoms. The atomic finite element with three nodes is considered. The atomic chain is bounded so that to the left side of atom 1 and to the right side of atom 4 there are no atoms. Because of this lack of neighboring atoms, the atomic finite elements on the edges of the mesh must be modified. So even before the imposition of the boundary conditions, the potential energy of the edge element must be modified to account for the nonexistence of one or more neighboring atoms.
Using a 3node atomic finite element, the assemblage procedure for the 4 element chain, shown in Figure 8, requires 4 atomic finite elements, called, respectively, El=1, El=2, El=3, El=4. Initially, all the atomic elements have three local nodes, designated by (1), (0) and (+1). In the sequence, the atomic elements, El=1 and El=4, must be modified. For the atomic El = 1, the contribution of energy from the central node (0) onto left node (1) does not exist and, therefore, must be eliminated from the calculations of the stiffness matrix and global nonequilibrium force vector. Also, the degree of freedom of node (1) of EI=1 must be eliminated from the assemblage procedure. This is graphically indicated by the dashed arrow between node (0) and node (1) for EI=1. Analogously, for EI=4, the energy contribution of node (0) to the node (+1) must be eliminated from the energy calculation procedure. Again, the dashed arrow in EI=4, between node (0) and node (+1) illustrates the required modification.
The expression for the total energy of the 4 atom chain, Utotal, is given by the sum of the contributions of the energy of individual elements, U^{EI}, as shown in equation (22):
The the expressions for the energy of each of the 4 atomic elements used to model the 4atom chain of Figure 8 are given in the expressions (23) to (26):
It should be noticed that the energy expressions for the elements in equations (23) to (26) are given in terms of the distance between the atoms, considering the global numbering scheme. Elements EI=1 and EI=4 only present the energy contribution from one bond. On the other hand, elements EI=2 and EI=3 present two energy contributions, because these elements present neighboring atoms to the right and to the left.
The third important stage of the assemblage produce in AFEM is the application of the boundary conditions. Figure 9 illustrates the same atomic chain of Figure 8 but including a set of boundary conditions. In this case, the displacement of atom 1 is blocked, so, the degree of freedom of this atom must be excluded from the global stiffness matrix, and from of the global nonequilibrium force vector.
The assemblage procedure of AFEM in two dimensions system is the same used for the one dimensional system. Figure 10 illustrates an example, which consists of a triangular (hexagonal) mesh with 38 atoms. The basic atomic finite element considered has 7 nodes, and was shown in Figure 5. In order to apply boundary conditions or to analyze meshes with defects and missing atoms, the basic atomic elements must also be modified to account for the missing neighboring atoms. This is illustrated in Figure 10. Element number 23 is an example of a complete atomic element. Elements 1, 6 and 25 are incomplete, or modified elements. Just like the 1D case, the bond energy of the missing neighboring atoms must be excluded from the total energy of the element.
The total energy expression for the mentioned elements 23, 1, 6, and 25 are:
The total bond energy of the complete element 23 and of the modified elements 1, 6 and 25 is given in Table 1. The element energy determined within the atomistic finite element procedure and the corresponding energy obtained by the Molecular Dynamics (MD) package LAMMPS (Plimpton, 1995Plimpton, S. (1995), Fast Parallel Algorithms for ShortRange Molecular Dynamics, Journal of Computational Physics 117: 119.) are also indicated. It can be seen that the energy for each AFE is distinct, according to the number of missing surrounding atoms. Table 1 also shows that the energy values calculated by AFEM and MD are in accordance.
The matrices and load vectors for the complete or modified elements are assembled and, in the sequence, the boundary conditions are imposed. It is stressed again, that modifying the atomic elements to account for missing bonds and the imposition of boundary conditions are two distinct solution steps within the AFEM. In the present implementation the resulting nonlinear system, equation (7) was solved iteratively by the NewtonRaphson method.
6 RESULTS
6.1 Comparison of AFEM and MD in One Dimension
In this section, the 1D AFEM is investigated. Initially, the influence of the number of element nodes in the atomic elements is addressed. As already discussed, the number of nodes in the AFE is related to the chosen cutoff radius for the Lennard Jones potential. The first example is shown in Figure 11. It consists of a 7 atom chain, fixed at atom 1 and loaded by a force F at node 7. This example will be solve using atomic elements with 3, 5 and 7 nodes. Before imposing the boundary conditions, the distance between the atoms is set to be equal to the equilibrium position, r_{eq} = 1.122 Å. After applying the boundary conditions, the load vector is increased until the atomic chain fails by bond breaking. At every load step the nonlinear equation (7) is solved iteratively.
Figure 12 shows the resulting forcestrain relations. The strain measure ( is obtained by dividing the displacement of atom 7, u_{7} by the original chain length L_{o}, ( = u_{7}/L_{o}. The curves show that there is very little difference between the results obtained considering the atomic elements with 5nodes and 7nodes. It suggests that 1D AFEM simulations using atomic finite element with 5 nodes is almost as accurate as atomic finite element with 7nodes.
In the sequence, the accuracy of the obtained AFEM results will be validated by comparison with the results obtained from MD package LAMMPS. To simplify, the values of parameters of the LJ potential were assigned unit values, σ = 1 and ε_{LJ} =1. The AFEM results were determined for the 5 node element. In the MD simulations, the Lennard Jones potential was used at a temperature of 1 K. Nonperiodic boundary conditions were used. For the AFEM and MD simulations the cutoff radius is set equal to r_{c} = 2.5σ (Smit, 1992Smit, B. (1992), Phase diagrams of LennardJones fluids, The Journal of Chemical Physics 96: 86398640.). The time integration step for the MD simulations is 0.05 fs.
Figure 13a shows forcestrain curves obtained by both methodologies. Figure 13c depicts the relative error between strains obtained by both solutions. There is a good agreement between both solutions, with a relative error smaller that 10^{5}, except at the initial equilibrium position. The MD simulation needed 1000 iterations to achieve the equilibrium position for every load step. On the other hand, the number of iterations for the AFEM to reach a nonequilibrium force vector smaller than 10^{11} for distinct values of the load force is shown in Figure 13b. The number of iterative steps for every force in the AFEM varied from 4 to 8 according to the loading stage.
1D AFEM results  a) forcestrain relation; b) number of iterations for AFEM, c) MDAFEM relative error, d) convergence measure.
Figure 13d shows a typical behavior of the quantity u•P, that is the dot product between displacement vector and the global nonequilibrium force vector. It also represents a measure of convergence of the iterative solution (Kim, 2006Kim, K. (2006). The Atomicscale Finite Element Method for Analyzing Mechanical Behavior of Carbon Nanotube and Quartz, M.Sc. Dissertation, Virginia Polytechnic Institute, Virginia.). The figure shows that convergence obtained for the present results is very fast.
The results of Figures 13a to 13d suggest that the 1D AFEM formulation and implementation is correct, that the results are accurate and the number of iterative steps, required for the solution, is significant smaller than those of the typical MD simulations.
6.2 Comparison of AFEM and MD in Two Dimensions
This section is dedicated to the study of a 2D AFE. A twodimensional hexagonal (or triangular) lattice, shown in Figure 14, presenting 655 atoms is analyzed. The 2D atomic finite element (AFE) with 7 nodes, given in Figure 5, is considered. Zero displacement boundary conditions are imposed on all atoms on the left edge, shown inside of the red dashed box. External force with the same amplitude are applied on all atoms on the right edge in x direction, as shown by the arrows in Figure 14. The Lennard Jones potential was used in MD at a temperature of 1 K. The LennardJones potential parameters are, σ = 1= and ( _{LJ} = 1. The cutoff radius for the MD simulations is set to r_{c} = 1.5Å. Nonperiodic boundary conditions were used.
(a) Forcestrain relation; (b) forcedisplacement relation; (c) convergence of AFEM; (d) number of iterations for distinc force levels.
Figure 15a shows a forcestrain diagram obtained by AFEM and MD. The results are obtained for the central atom at the left border of the mesh. A similar result, showing the change in total length of the atomic mesh, ∆L, is given if Figure 15b. The AFEM and MD solutions agree very well. The residual (u•P) is shown in Figure 15c and the number of iterations required to solve the nonlinear system for every load step can be found in Figure 15d. The MD simulation required 50000 iterative steps to achieve a final configuration with 4 significant digits, whereas the AFEM never needed more than 7 steps to achieve the same accuracy.
As a second 2D example, consider the same mesh of the previous example, shown in Figure 16. Zero displacement boundary condition imposed to the atoms within the red dashed box at the left border. An external force F_{ext} is applied in the vertical direction upon atom number 739. The black lattice in Figure 16 represents the original atomic mesh. The blue mesh shows the deformed mesh at its last loading step before bond breaking. The atom 739 is connected to atoms 728 and 729 as illustrated in the detail shown in Figure 16. It is simple to see that the bonding force between atoms 739 and 729 is in compression and that the force between atoms 739 and 728 is traction.
Figure 17 shows in dashed lines the interatomic bonding force developed by the LennardJones potential as a function of the interatomic distance r _{ij} , as anticipated by equation (11). This figure also show the bonding forces between atoms 739729 (compression) and between atoms 739728 (traction) for distinct amplitude stages of the external force, F_{ext}. Traction is marked in blue color and compression in red color.
As expected, the interatomic bonding forces start with a zero value at the equilibrium position and increase continuously in traction and compression for increasing values of the external load. The bonding forces follow exactly the path predicted in equation (11). It should be noticed that this is a highly nonlinear behavior, in which interatomic displacements for compressed bonds is much smaller than those of the traction bonds for the same external loading level. Figure 17 also shows the break bonding situation. After the condition, in which the traction bond, in blue, reaches the maximum bonding force at r_{fmax} = 1.2445 σ, given by equation (13), the bonding forces decrease for larger distances, the bond can no longer withstand the imposed external loads and breaks. The iterative system will no longer converge, unless a new mesh without the lose atom is considered.
In the third example, the effect of missing atoms or vacancies on the mechanical behavior of an atomic mesh is investigated. Figure 18a shows a pristine mesh, with no vacancies. On the other hand, Figure 18b depicts the same mesh but with 4 missing atoms. In both cases, zero displacement boundary conditions are prescribed for the atoms in the right edge within the red dashed box. A force of equal, increasing amplitude is applied to the 4 atoms at the left border. The external load is increased incrementally until one or more atomic bonds within the mesh break. Figure 18c shows the stress (force)strain behavior of the pristine and of the defect mesh. The same strain measure of the previous examples is used. It can be seen that vacancies do have a large effect on the forcestrain behavior. For validation purposes, Figure 18d depicts the deformed configuration of the defect mesh (Figure 16b), determined by AFEM and MD. There is a good visual agreement between the two methodologies.
(a) Pristine mesh; (b) Defected mesh, (c) stressstrain curves, and (d) final equilibrium configuration from AFEM and MD.
As a last example, consider the square Bravais lattice shown in Figure 19a. The square 2D AFE with 5 nodes, shown in Figure 6a, is used to model the atomic domain. The displacement of all atoms with the box at the left edge of the mesh are blocked. Upon all atoms within the dashed box at the right side, an equal and increasing force is applied. Figure 19b shows the resulting forcestrain curves obtained by both methodologies. Both curves show the same trend, but there is no perfect agreement. This should be further investigated, because all previous and even more complex meshes had presented a much better agreement between both methodologies.
7 CONCLUSIONS
This work presents a detailed description of the formulation and implementation of the Atomistic Finite Element Method (AFEM), exemplified in the analysis of one and twodimensional atomic domains governed by the LennardJones interatomic potential. The methodology to synthesize element stiffness matrices and load vectors, the potential energy modification of the atomistic finite elements (AFE) to account for boundary edge effects, the inclusion of boundary conditions were carefully described, in a way that the authors had not previously found in the AFEM literature. The conceptual relation between the cutoff radius of interatomic potentials and the number of nodes in the AFE is addressed and exemplified for the 1D case. For the 1D case elements with 3, 5 and 7 nodes were addressed. The AFEM has been used to describe the mechanical behavior or onedimensional atomic arrays as well as twodimensional lattices of atoms. The reported studies also included the analysis of pristine domains, as well as domains with missing atoms, defects, or vacancies. Almost all results were compared with classical molecular dynamic simulations (MD) performed using a commercial package. The results have been very encouraging in terms of accuracy and in the computational effort necessary to execute both methodologies, AFEM and MD. The methodology presented can be extended to other potential, like Tersoff (Tersoff, 1987Tersoff, J. (1987), New empirical approach for the structure and energy of covalent systems, Physical Review B 37: 6991.), REBO (Brenner, 1990Brenner, D.W. (1990), Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Physical Review B 42: 9458.), AIREBO (Stuart, 2000Stuart, S.J.; Tutein, A.B.; Harrison, J.A. (2000), A reactive potential for hydrocarbons with intermolecular interactions, The Journal of Chemical Physics 112: 6472 6486.), which have been applied to model nano structures composed of graphene and materials alike.
Acknowledgments
The research leading to this article has been funded by the São Paulo Research Foundation (Fapesp) through grants 2012/179484, 2013/230851, 2015/002092 and 2013/082937 (CEPID). Support from CAPES and CNPQ has been fundamental to the development of the present research. This is gratefully acknowledged.
References
 Alder, B.J., Wainwright, T.E.J. (1959). Studies in Molecular Dynamics. I. General Method, The Journal of Chemical Physics 31: 459466.
 Brenner, D.W. (1990), Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Physical Review B 42: 9458.
 Chen, C., Hone, J. (2013). Graphene Nanoelectromechanical Systems. Proceedings of the IEEE 101(7).
 Damasceno, D.A., Rajapakse, R.K.N.D., Mesquita, E. (2015). Mechanical behavior of atoms using atomicscale finite element method. Proceedings of the XXXVI IberoLatin American Congress on Computational Methods in Engineering.
 Dingerville, R., Qu, J., Cherakoui, M. (2005). Surface free energy and its effect on the elastic behavior of nanosized particles wires and films, Journal of the Mechanics and Physics of Solids 53: 18271854.
 Dirac, P.A.M. (1981). The Principles of Quantum Mechanics, Clarendon Press, Oxford.
 Gurtin, M., Ian Murdoch, A. (1975). A continuum theory of elastic material surfaces, Archive for Rational Mechanics and Analysis 57: 291323.
 Jones, J.E. (1924). On the determination of molecular fieldsII. From the equation of state of a gas. Proceedings of the Royal Society A 106: 463477.
 Kim, K. (2006). The Atomicscale Finite Element Method for Analyzing Mechanical Behavior of Carbon Nanotube and Quartz, M.Sc. Dissertation, Virginia Polytechnic Institute, Virginia.
 Liu, B., Huang, Y., Jiang, H., Qu, S., Hwang, K.C. (2004). The atomicscale finite element method, Computer Methods in Applied Mechanics and Engineering 193: 18491864.
 Liu, B., Jiang, H., Huang, Y., Qu, S., Yu, M.F. (2005). Atomicscale finite element method in multiscale computation with applications to carbon nanotubes, Physical Review B, 72: 10980121.
 Liu, C., Rajapakse, R.K.N.D. (2010). Continuum models incorporating surface energy for static and dynamic response of nanoscale beams, IEEE Transactions on Nanotechnology 9: 422431.
 Liu, W.K., Karpov, E.G., Park, H.S. (2006). Nano Mechanics and Materials: Theory, Multiscale Methods and Applications, John Wiley & Sons.
 Malakouti, M., Montazeri, A. (2016). Nanomechanics analysis of perfect and defected graphene sheets via a novel atomicscale finite element method, Superlattices and Microstructures 94: 1 12.
 Miller, R.E., (2000). Sizedependent elastic properties of nanosized structural elements, Nanotechnology, 11:139147.
 Nanoscience and nanotechnologies: opportunities and uncertainties (2004). The Royal Society and The Royal Academy of Engineering, London.
 Novoselov, K.S., Fal′Ko, V.I., Colombo, L., Gellert, P.R., Schwab, M.G., Kim, K. (2012). A roadmap for graphene, Nature 490: 192.
 Plimpton, S. (1995), Fast Parallel Algorithms for ShortRange Molecular Dynamics, Journal of Computational Physics 117: 119.
 Sapsathiarn, Y., Rajapakse, R.K.N.D. (2012). A model for large deflections of nanobeams and experimental comparison, IEEE Transactions on Nanotechnology 11, 247254.
 Smit, B. (1992), Phase diagrams of LennardJones fluids, The Journal of Chemical Physics 96: 86398640.
 Stuart, S.J.; Tutein, A.B.; Harrison, J.A. (2000), A reactive potential for hydrocarbons with intermolecular interactions, The Journal of Chemical Physics 112: 6472 6486.
 Tadmor, E.B., Miller, R.E. (2011). Modeling materials: continuum, atomistic, and multiscale techniques. Cambridge University Press, New York.
 Tao, J., Sun, Y., (2016). The Elastic Property of Bulk Silicon Nanomaterials through an Atomic Simulation Method. Journal of Nanomaterials 2016: 6.
 Tersoff, J. (1987), New empirical approach for the structure and energy of covalent systems, Physical Review B 37: 6991.
Publication Dates

Publication in this collection
Dec 2017
History

Received
21 May 2017 
Accepted
03 Aug 2017