- Email: [email protected]
A pseudo-atom method in potential energy minimization for crystals
Computational Materials Science 44 (2008) 547–551
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
A pseudo-atom method in potential energy minimization for crystals Oleg Vinogradov * Department of Mechanical and Manufacturing Engineering, University of Calgary, 2500 University Drive NW, Calgary, Alberta, Canada T2N 1N4
a r t i c l e
i n f o
Article history: Received 31 March 2008 Received in revised form 17 April 2008 Accepted 22 April 2008 Available online 20 June 2008 PACS: 81.07.Bc 62.25.-g Keywords: Atomistic simulation Molecular statics Energy minimization Computational efficiency
a b s t r a c t Three energy minimization algorithms analogous to molecular dynamics, when each atom is displaced independently in the direction of the unbalanced force, are presented and compared. The displacements are governed by the local potential field, defined by the atom’s neighbours, and it is assumed that its collective effect corresponds to that of a single pseudo-atom. In one of the algorithms the atom at each iteration step is moved directly into the position of the local minimum, and in the other two it moves toward this position according to either a linear law or an exponential one. The examples for 2D hexagonal Lennard–Jones lattices with different numbers of atoms are used in this study. It is shown that the computational efficiency is of the order O(N2) for all three approaches but the rate of convergence for the algorithm with an exponential law of displacements gives the best results for the examples considered. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Traditionally, molecular statics (MS) methods are based on various iterative energy minimization algorithms (like steepest descent, conjugate gradient, and Newton’s method) for the entire Natoms system and the equilibrium state of the lattice corresponding to the global minimum of its potential energy is found by solving a coupled system of N equations. Since the number of minima on the potential energy surface increases exponentially with the size of the system [1], finding the global minimum for a large N represents a computationally challenging problem, and in large scale simulations it is usually solved with some approximation to the ‘‘real” minimum or a saddle point on the above surface [2]. For a coupled system of N equations the computational efficiency of convergence is usually of the order of either O(N3) or O(N2) for the classical energy minimization algorithms [2], and of O(N2) for an alternative approach developed in [3]. Even more important, in all gradient-type methods the initial conditions chosen might prevent finding a minimum associated with dislocation formation since the iterative algorithm may not be able to overcome the corresponding potential barrier (the solution may be trapped in a local minimum). For example, in [4] the MS technique was used up to the onset of dislocation nucleation
* Tel.: +1 403 220 7187; fax: +1 403 282 8406. E-mail address: [email protected] 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.04.028
and then the molecular dynamics (MD) to find the equilibrium state during dislocation nucleation. To uncouple the system of equations, an algorithm was suggested in [5], which was based on an ordered sequence of minimization steps such that at each step the unbalanced resultant forces acting on atoms were found, and then each atom was moved (similar to MD method) into the direction of the unbalanced force vector by the amount governed by an exponential function, which depended on two parameters: the magnitude of the unbalanced force and the displacement of this atom from the local equilibrium. It was assumed in [5] that the unbalanced force was in effect the force between two atoms: the given atom and the imaginary (or a pseudo) atom, the latter representing the collective effect of all neighbours within the cut-off distance. The method was applied to study dislocations formation in single [5] and polycrystal systems [6]. It was shown that it was capable of finding the equilibrium state with any prescribed accuracy [7]. In addition to its robustness, the above method, since the atoms are moved independently, allows parallel processing, similar to MD algorithms. For brevity, we will call this approach exponentially controlled atom displacement (ECAD) method. The algorithm suggested in [5] was, in fact, based on local energy minimization for a given atom among its fixed neighbours, and at each step this atom was moved into the direction of the local minimum without reaching this minimum. The question then arises whether an algorithm in which the given atom is moved into the local minimum of the potential energy exactly will be more
548
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
computationally efficient. This local minimum can be found by using conventional minimizations methods. Correspondingly, the global minimum then is found by an iterative local energy minimization (LEM) scheme. Another alternative to the algorithm described in [5] is to define the atom displacement as a fraction of the distance to the local equilibrium corresponding to the given unbalanced force. Since in this case the function describing atom displacements is linearly, instead of exponentially [5], depends on the distance to the equilibrium, the convergence to equilibrium might be faster. We will call it below a linearly controlled atom displacement (LCAD) method. The objective of this paper is to explore these two alternative approaches in terms of their computational efficiencies and compare the results against the one suggested earlier [5]. The numerical experiments are done for various 2D Lennard–Jones (LJ) hexagonal lattices.
The main idea of defining displacements by Eq. (4) was to move each atom into the direction of minimum local potential in steps smaller then the maximum possible step corresponding to the above minimum. As the first alternative, the displacement of the atom can be found by minimizing its local energy potential, assuming that its neighbours are fixed. The potential function for an atom i among its fixed neighbours is
2. The algorithms
aij ¼ c rij r0
In the following the LJ pair-wise interaction potentials between the atoms are used in a normalized form
where 0 < c < 1. Similar to the algorithm given in [5], at each iterative step, for each approach, a subset of atoms nk 2 N is identified and then each atom is moved into a new position based on the specific definition of aij.
Uðxi ; yi ; zi Þ ¼
ð1Þ
where / ¼ U=4e, r ij ¼ Rij =r, e defines the strength of the interaction, r is the distance scaling factor, Rij is the distance between the atoms i and j, 0 < r ij 6 rc , and rc is a cut-off distance, which is defined here by the second-nearest neighbours. The corresponding interatomic normalized force between the atoms i and j is given by
fij ðr ij Þ ¼ 12r 13 þ 6r 7 ij ij
ð2Þ
It follows from the above equation that the maximum normalized force for a pair of atoms is achieved at the normalized distance rmax = (26/7)0.5 = 1.224, and then the corresponding maximum normalized force is fmax = 0.599. The normalized distance between the atoms at equilibrium is equal to r0 = 1.122. We consider a system of N-atoms with coordinates ui, (i = 1,2,. . .,N) making up a vector u. The iterative processes is described by the following relationship
ukþ1 ¼ uk þ ak pk
/ij ðrij Þ
ð5Þ
j¼1
where m is the number of atoms within the cut-off distance of the atom i. As the second alternative, a fraction of the equilibrium distance between the atom i and the pseudo-atom can be used to define aij. Namely, aij in this case is given by
h i /ij ðr ij Þ ¼ r 12 r 6 ij ij
m X
ð6Þ
3. Results The program was written in Mathematica 5.2, using DELL Precision 670 computer. The numerical experiments were carried out on various 2D samples of hexagonal closed-packed lattices with one deficiency. The length-to-width ratio of samples was within 2–2.1 and the total numbers of atoms in samples were as follows: 275, 940, 1414, 1984, 2979, 3785, 4966, 6307, and 7808. An example of a 275-atom lattice is shown in Fig. 1. In applying the LEM scheme, the system of two equations (for a 2D case) Eq. (5) was solved by using FindRoot command, where (xi, yi) were the unknowns, and the starting point (xio, yio) were the coordinates of the moving atom at the previous iteration (or initial
ð3Þ
where pk is the unit vector defining the directions of displacement from the point uk, and ak is the vector of scalar factors defining the magnitudes of these displacements. In any intermediate state of relaxation the resultant force acting on each atom is non-zero. Thus if aij ak , where i denotes a real atom number and j denotes a pseudo-atom, is positive then the atom i will be moving towards the equilibrium position with the atom j. The process of convergence and its efficiency depend on how the vector ak is defined. In [5] the coefficients aij were defined by the following exponential function
aij ¼ a 1 ebjfij jjrij r0 j
ð4Þ
where jfij j fk is the magnitude of the unbalanced force acting on the atom i, fk is the vector of unbalanced forces acting on a subset of atoms nk N at the k th iteration step, r ij is the distance between two atoms corresponding to the magnitude of the interatomic force |fij| between the given atom i and a pseudo-atom, representing the effect of all given atom neighbours within the cut-off distance, r0 is the normalized distance between the two atoms in equilibrium, and a and b are constants. Note that the subset of atoms nk is defined by all those atoms for which the unbalanced forces exceed the minimum force characterising the convergence criterion.
Fig. 1. A basal plane of 275-atoms hexagonal close-packed lattice with deficiency.
549
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
position). The FindRoot command in our case is based on using the Newton’s method to find a numerical solution to the coupled system of equations (Eq. (5)) characterising the local potential of atom i among its fixed neighbours within the cut-off distance. We used a default, 10 digits, setting for the accuracy. Each crystal was subjected to tensile boundary displacements. Accordingly, its upper boundary, more specifically, a layer of atoms corresponding to a cut-off distance was fixed, while a similar layer at the bottom was subjected to incremental displacements. Thus the set of all atoms in a crystal, S0, was divided into two sets: Sf (atoms belonging to a fixed or a displaced boundary) and a complementary set, Sm = S0 Sf (all moving atoms). The (x,y)-coordinate system is placed at the lower left corner of the crystal before relaxation, where x is horizontal and y is vertical. For any atom i 2 ðSf Þ its trajectory ui(x,y) satisfies the boundary condition oui ðx; yÞ=ox ¼ 0, whereas for all atoms i 2 Sm , including left and right boundaries, it will be that oui ðx; yÞ=ox–0. Note that the atoms in the Sm set are moving according to Eqs. (3) and (4). The lower boundary displacements were done in equal steps of 0.4(rmax r0). The iterative process at each loading step was stopped when the resultant force acting on any atom was smaller than 0.00083fmax. In Fig. 2 the normalized stress vs strain data are shown for the 275-atoms lattice for the three approaches described above. The normalized stress and strain correspond to the engineering definitions, namely, in our case, for stress it is force divided by the original width of the crystal, and for strain it is displacement divided by the original length of the crystal. In Fig. 2 the dots indicate the results for each loading step. As one can see, the three approaches produce identical results. The sudden drop in stress corresponds to nano-instabilities in the system associated with dislocations formation. Again, it was found that for the three approaches the lattices before and after dislocation formations were identical. The results are shown in Fig. 3, in which the network of scaled forces (tensile and compressive) is shown before and after dislocation formations. As one can see, the nanocavities produced by dislocations result in local inhomogeneous
Fig. 3. The network of interatomic forces (scaled) just before (a) and during (b) dislocations formation (red – compressive and blue – tensile). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
atomic force field. As it has been shown previously [5], the process of convergence to equilibrium at the loading step just before and during dislocation formations is qualitatively different since
0.12
0.1
Stress
0.08
0.06
0.04
0.02
0.02
0.04
0.06
0.08
Strain Fig. 2. Stress vs. strain data for a 275-atoms lattice (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
550
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
1.2
1
Error, %
0.8
0.6
0.4
0.2
0 0
5000
10000
15000
20000
25000
30000
Iteration number Fig. 4. Convergence to equilibrium just before (red) and during (black) dislocations formation for a 275-atoms lattice. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
4.5 4.25 4
Time, h
during dislocation formation it is necessary to overcome potential barriers. This is shown in Fig. 4 for a 275-atoms lattice. We should note that all three approaches were able to overcome the above barriers during dislocation formations and had similar convergence patterns. In order to compare the three different minimization schemes, the empirical coefficients in ECAD and LEM methods should be optimised in terms of computational efficiency. The coefficient b in Eq. (4) was taken to be the same, b = 500, since its effect on the rate of convergence was relatively small. The effect of the coefficient a, on the other hand, is more pronounced, as can be seen in Fig. 5 for a 275-atoms lattice at the first step of loading. The larger this coefficient is the smaller the convergence time becomes. However, large coefficients a lead to larger atoms displacements and since these displacements are independent from each other at each iteration step the convergence can fail at other steps during the loading, especially during dislocations formation. It was found that for this particular system this coefficient should be around 1.2 in order to simulate the entire loading process shown in Fig. 2. In the following it was taken that a = 1.2. Similar tests were performed for the c coefficient in Eq. (6) using the 4966-atoms lattice. The results are shown in Fig. 6. In
3.75 3.5 3.25 3 0.3
0.4
0.5
0.6
γ
0.7
0.8
0.9
1
Fig. 6. The effect of coefficient c on a convergence time for a 4966-atoms lattice during the first step of loading.
8
Time, h
6
1.75 1.5
4
2
Time, h
1.25 0
1
0
2000
4000
6000
8000
Number of atoms
0.75
Fig. 7. Convergence times vs number of atoms during the first loading step (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method); broken lines indicate quadratic approximations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.5 0.25 0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Coefficient a Fig. 5. The effect of coefficient a on a convergence time for a 275-atoms lattice during the first step of loading.
the following c = 0.85, resulting in the shortest computer time, was chosen for the comparative analysis.
551
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
2500
Subset Length
2000
1500
1000
500
0 0
200
400
600
800
Iteration number Fig. 8. Variability of the length of the subset nk during convergence to equilibrium of a 6307-atoms lattice during the first step of loading (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
In Fig. 7 the convergence times at the first step of loading vs the number of atoms is shown for the three simulation algorithms. The broken line curves are the quadratic polynomials fits for each model. As one can see, for the number of atoms tested the convergence is of the order O(N2) but the rate of convergence is different for different models. The explanation is rooted in the variability of the subsystem of atoms nk, which defines all atoms with resultant forces exceeding the minimum allowable. The smaller this subsystem is the faster the convergence. Since atoms are moved into new positions independently, the most efficient process of self-adjustment is the one producing at each iteration the smallest subsets nk. In Fig. 8 the variability of the subsets nk 2 Sm during the convergence to equilibrium for the 6307-atoms lattice during the first step of loading is shown. Note that the total number of atoms in the moving set was Sm = 6045. As one can see, the LCAD and LEM methods produce qualitatively different results as those for the ECAD one. The first two cause rapid growth and then decay of the subsystem of atoms with resultant forces exceeding convergence criterion, whereas the ECAD method leads to more steady change and much smaller subsystems. In fact, for the ECAD method the number of atoms in nk is 10-fold smaller than the total number of atoms in the relaxation subset Sm. This explains the better efficiency of the ECAD method in Fig. 7 compared to the two other alternatives. This property of the ECAD method, in addition to its ability to utilize parallel processing, makes it potentially more advantageous compared to the classical minimization algorithms. 4. Conclusion The idea of a pseudo-atom, representing the collective effect of all neighbour atoms to a given one, was introduced and
investigated for three different numerical algorithms. The investigations were carried out for a system of LJ 2D crystals with different number of atoms and were limited to a tensile extension of these crystals. The results for a 275-atoms lattice show that all three algorithms produce identical stress–strain data and interatomic force fields, and capable of detecting dislocations formation. However, the algorithm introduced previously in [5] gives the best computational efficiency. This is explained by examining the variability of the subsystem of atoms with resultant forces exceeding the minimum allowable one during the system relaxation to equilibrium. It is shown that in a self-organized system the convergence to equilibrium by independent displacement of atoms is efficient but strongly depends on the magnitudes of atoms displacements. All three methods allow parallel processing, which can be utilized to further increase computational efficiency. The above numerical investigation was limited to a specific range of system sizes, to planar lattices, and a specific loading. Further applications of the method suggested in [5] to larger 3D systems will allow to clarify its potential for simulations of more general molecular statics problems. References [1] D.J. Wales, H.A. Scheraga, Science 285 (5432) (1999) 1368–1373. [2] A.R. Leach, Molecular Modelling: Principles and Applications, Prentice Hall, 2001. [3] I. Ye. Telitchev, O. Vinogradov, Int. J. Comp. Methods 3 (1) (2006) 71–81. [4] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Phil. Mag. 64 (4) (1990) 851–878. [5] O. Vinogradov, Int. J. Comp. Methods 3 (2) (2006) 153–161. [6] O. Vinogradov, Comp. Mater. Sci 39 (3) (2007) 611–615. [7] O. Vinogradov, Comp. Mater. Sci. 42 (2008) 478–482.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
A pseudo-atom method in potential energy minimization for crystals Oleg Vinogradov * Department of Mechanical and Manufacturing Engineering, University of Calgary, 2500 University Drive NW, Calgary, Alberta, Canada T2N 1N4
a r t i c l e
i n f o
Article history: Received 31 March 2008 Received in revised form 17 April 2008 Accepted 22 April 2008 Available online 20 June 2008 PACS: 81.07.Bc 62.25.-g Keywords: Atomistic simulation Molecular statics Energy minimization Computational efficiency
a b s t r a c t Three energy minimization algorithms analogous to molecular dynamics, when each atom is displaced independently in the direction of the unbalanced force, are presented and compared. The displacements are governed by the local potential field, defined by the atom’s neighbours, and it is assumed that its collective effect corresponds to that of a single pseudo-atom. In one of the algorithms the atom at each iteration step is moved directly into the position of the local minimum, and in the other two it moves toward this position according to either a linear law or an exponential one. The examples for 2D hexagonal Lennard–Jones lattices with different numbers of atoms are used in this study. It is shown that the computational efficiency is of the order O(N2) for all three approaches but the rate of convergence for the algorithm with an exponential law of displacements gives the best results for the examples considered. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Traditionally, molecular statics (MS) methods are based on various iterative energy minimization algorithms (like steepest descent, conjugate gradient, and Newton’s method) for the entire Natoms system and the equilibrium state of the lattice corresponding to the global minimum of its potential energy is found by solving a coupled system of N equations. Since the number of minima on the potential energy surface increases exponentially with the size of the system [1], finding the global minimum for a large N represents a computationally challenging problem, and in large scale simulations it is usually solved with some approximation to the ‘‘real” minimum or a saddle point on the above surface [2]. For a coupled system of N equations the computational efficiency of convergence is usually of the order of either O(N3) or O(N2) for the classical energy minimization algorithms [2], and of O(N2) for an alternative approach developed in [3]. Even more important, in all gradient-type methods the initial conditions chosen might prevent finding a minimum associated with dislocation formation since the iterative algorithm may not be able to overcome the corresponding potential barrier (the solution may be trapped in a local minimum). For example, in [4] the MS technique was used up to the onset of dislocation nucleation
* Tel.: +1 403 220 7187; fax: +1 403 282 8406. E-mail address: [email protected] 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.04.028
and then the molecular dynamics (MD) to find the equilibrium state during dislocation nucleation. To uncouple the system of equations, an algorithm was suggested in [5], which was based on an ordered sequence of minimization steps such that at each step the unbalanced resultant forces acting on atoms were found, and then each atom was moved (similar to MD method) into the direction of the unbalanced force vector by the amount governed by an exponential function, which depended on two parameters: the magnitude of the unbalanced force and the displacement of this atom from the local equilibrium. It was assumed in [5] that the unbalanced force was in effect the force between two atoms: the given atom and the imaginary (or a pseudo) atom, the latter representing the collective effect of all neighbours within the cut-off distance. The method was applied to study dislocations formation in single [5] and polycrystal systems [6]. It was shown that it was capable of finding the equilibrium state with any prescribed accuracy [7]. In addition to its robustness, the above method, since the atoms are moved independently, allows parallel processing, similar to MD algorithms. For brevity, we will call this approach exponentially controlled atom displacement (ECAD) method. The algorithm suggested in [5] was, in fact, based on local energy minimization for a given atom among its fixed neighbours, and at each step this atom was moved into the direction of the local minimum without reaching this minimum. The question then arises whether an algorithm in which the given atom is moved into the local minimum of the potential energy exactly will be more
548
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
computationally efficient. This local minimum can be found by using conventional minimizations methods. Correspondingly, the global minimum then is found by an iterative local energy minimization (LEM) scheme. Another alternative to the algorithm described in [5] is to define the atom displacement as a fraction of the distance to the local equilibrium corresponding to the given unbalanced force. Since in this case the function describing atom displacements is linearly, instead of exponentially [5], depends on the distance to the equilibrium, the convergence to equilibrium might be faster. We will call it below a linearly controlled atom displacement (LCAD) method. The objective of this paper is to explore these two alternative approaches in terms of their computational efficiencies and compare the results against the one suggested earlier [5]. The numerical experiments are done for various 2D Lennard–Jones (LJ) hexagonal lattices.
The main idea of defining displacements by Eq. (4) was to move each atom into the direction of minimum local potential in steps smaller then the maximum possible step corresponding to the above minimum. As the first alternative, the displacement of the atom can be found by minimizing its local energy potential, assuming that its neighbours are fixed. The potential function for an atom i among its fixed neighbours is
2. The algorithms
aij ¼ c rij r0
In the following the LJ pair-wise interaction potentials between the atoms are used in a normalized form
where 0 < c < 1. Similar to the algorithm given in [5], at each iterative step, for each approach, a subset of atoms nk 2 N is identified and then each atom is moved into a new position based on the specific definition of aij.
Uðxi ; yi ; zi Þ ¼
ð1Þ
where / ¼ U=4e, r ij ¼ Rij =r, e defines the strength of the interaction, r is the distance scaling factor, Rij is the distance between the atoms i and j, 0 < r ij 6 rc , and rc is a cut-off distance, which is defined here by the second-nearest neighbours. The corresponding interatomic normalized force between the atoms i and j is given by
fij ðr ij Þ ¼ 12r 13 þ 6r 7 ij ij
ð2Þ
It follows from the above equation that the maximum normalized force for a pair of atoms is achieved at the normalized distance rmax = (26/7)0.5 = 1.224, and then the corresponding maximum normalized force is fmax = 0.599. The normalized distance between the atoms at equilibrium is equal to r0 = 1.122. We consider a system of N-atoms with coordinates ui, (i = 1,2,. . .,N) making up a vector u. The iterative processes is described by the following relationship
ukþ1 ¼ uk þ ak pk
/ij ðrij Þ
ð5Þ
j¼1
where m is the number of atoms within the cut-off distance of the atom i. As the second alternative, a fraction of the equilibrium distance between the atom i and the pseudo-atom can be used to define aij. Namely, aij in this case is given by
h i /ij ðr ij Þ ¼ r 12 r 6 ij ij
m X
ð6Þ
3. Results The program was written in Mathematica 5.2, using DELL Precision 670 computer. The numerical experiments were carried out on various 2D samples of hexagonal closed-packed lattices with one deficiency. The length-to-width ratio of samples was within 2–2.1 and the total numbers of atoms in samples were as follows: 275, 940, 1414, 1984, 2979, 3785, 4966, 6307, and 7808. An example of a 275-atom lattice is shown in Fig. 1. In applying the LEM scheme, the system of two equations (for a 2D case) Eq. (5) was solved by using FindRoot command, where (xi, yi) were the unknowns, and the starting point (xio, yio) were the coordinates of the moving atom at the previous iteration (or initial
ð3Þ
where pk is the unit vector defining the directions of displacement from the point uk, and ak is the vector of scalar factors defining the magnitudes of these displacements. In any intermediate state of relaxation the resultant force acting on each atom is non-zero. Thus if aij ak , where i denotes a real atom number and j denotes a pseudo-atom, is positive then the atom i will be moving towards the equilibrium position with the atom j. The process of convergence and its efficiency depend on how the vector ak is defined. In [5] the coefficients aij were defined by the following exponential function
aij ¼ a 1 ebjfij jjrij r0 j
ð4Þ
where jfij j fk is the magnitude of the unbalanced force acting on the atom i, fk is the vector of unbalanced forces acting on a subset of atoms nk N at the k th iteration step, r ij is the distance between two atoms corresponding to the magnitude of the interatomic force |fij| between the given atom i and a pseudo-atom, representing the effect of all given atom neighbours within the cut-off distance, r0 is the normalized distance between the two atoms in equilibrium, and a and b are constants. Note that the subset of atoms nk is defined by all those atoms for which the unbalanced forces exceed the minimum force characterising the convergence criterion.
Fig. 1. A basal plane of 275-atoms hexagonal close-packed lattice with deficiency.
549
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
position). The FindRoot command in our case is based on using the Newton’s method to find a numerical solution to the coupled system of equations (Eq. (5)) characterising the local potential of atom i among its fixed neighbours within the cut-off distance. We used a default, 10 digits, setting for the accuracy. Each crystal was subjected to tensile boundary displacements. Accordingly, its upper boundary, more specifically, a layer of atoms corresponding to a cut-off distance was fixed, while a similar layer at the bottom was subjected to incremental displacements. Thus the set of all atoms in a crystal, S0, was divided into two sets: Sf (atoms belonging to a fixed or a displaced boundary) and a complementary set, Sm = S0 Sf (all moving atoms). The (x,y)-coordinate system is placed at the lower left corner of the crystal before relaxation, where x is horizontal and y is vertical. For any atom i 2 ðSf Þ its trajectory ui(x,y) satisfies the boundary condition oui ðx; yÞ=ox ¼ 0, whereas for all atoms i 2 Sm , including left and right boundaries, it will be that oui ðx; yÞ=ox–0. Note that the atoms in the Sm set are moving according to Eqs. (3) and (4). The lower boundary displacements were done in equal steps of 0.4(rmax r0). The iterative process at each loading step was stopped when the resultant force acting on any atom was smaller than 0.00083fmax. In Fig. 2 the normalized stress vs strain data are shown for the 275-atoms lattice for the three approaches described above. The normalized stress and strain correspond to the engineering definitions, namely, in our case, for stress it is force divided by the original width of the crystal, and for strain it is displacement divided by the original length of the crystal. In Fig. 2 the dots indicate the results for each loading step. As one can see, the three approaches produce identical results. The sudden drop in stress corresponds to nano-instabilities in the system associated with dislocations formation. Again, it was found that for the three approaches the lattices before and after dislocation formations were identical. The results are shown in Fig. 3, in which the network of scaled forces (tensile and compressive) is shown before and after dislocation formations. As one can see, the nanocavities produced by dislocations result in local inhomogeneous
Fig. 3. The network of interatomic forces (scaled) just before (a) and during (b) dislocations formation (red – compressive and blue – tensile). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
atomic force field. As it has been shown previously [5], the process of convergence to equilibrium at the loading step just before and during dislocation formations is qualitatively different since
0.12
0.1
Stress
0.08
0.06
0.04
0.02
0.02
0.04
0.06
0.08
Strain Fig. 2. Stress vs. strain data for a 275-atoms lattice (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
550
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
1.2
1
Error, %
0.8
0.6
0.4
0.2
0 0
5000
10000
15000
20000
25000
30000
Iteration number Fig. 4. Convergence to equilibrium just before (red) and during (black) dislocations formation for a 275-atoms lattice. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
4.5 4.25 4
Time, h
during dislocation formation it is necessary to overcome potential barriers. This is shown in Fig. 4 for a 275-atoms lattice. We should note that all three approaches were able to overcome the above barriers during dislocation formations and had similar convergence patterns. In order to compare the three different minimization schemes, the empirical coefficients in ECAD and LEM methods should be optimised in terms of computational efficiency. The coefficient b in Eq. (4) was taken to be the same, b = 500, since its effect on the rate of convergence was relatively small. The effect of the coefficient a, on the other hand, is more pronounced, as can be seen in Fig. 5 for a 275-atoms lattice at the first step of loading. The larger this coefficient is the smaller the convergence time becomes. However, large coefficients a lead to larger atoms displacements and since these displacements are independent from each other at each iteration step the convergence can fail at other steps during the loading, especially during dislocations formation. It was found that for this particular system this coefficient should be around 1.2 in order to simulate the entire loading process shown in Fig. 2. In the following it was taken that a = 1.2. Similar tests were performed for the c coefficient in Eq. (6) using the 4966-atoms lattice. The results are shown in Fig. 6. In
3.75 3.5 3.25 3 0.3
0.4
0.5
0.6
γ
0.7
0.8
0.9
1
Fig. 6. The effect of coefficient c on a convergence time for a 4966-atoms lattice during the first step of loading.
8
Time, h
6
1.75 1.5
4
2
Time, h
1.25 0
1
0
2000
4000
6000
8000
Number of atoms
0.75
Fig. 7. Convergence times vs number of atoms during the first loading step (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method); broken lines indicate quadratic approximations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.5 0.25 0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Coefficient a Fig. 5. The effect of coefficient a on a convergence time for a 275-atoms lattice during the first step of loading.
the following c = 0.85, resulting in the shortest computer time, was chosen for the comparative analysis.
551
O. Vinogradov / Computational Materials Science 44 (2008) 547–551
2500
Subset Length
2000
1500
1000
500
0 0
200
400
600
800
Iteration number Fig. 8. Variability of the length of the subset nk during convergence to equilibrium of a 6307-atoms lattice during the first step of loading (black dots – ECAD method, red dots – LCAD method, and blue dots – LEM method). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
In Fig. 7 the convergence times at the first step of loading vs the number of atoms is shown for the three simulation algorithms. The broken line curves are the quadratic polynomials fits for each model. As one can see, for the number of atoms tested the convergence is of the order O(N2) but the rate of convergence is different for different models. The explanation is rooted in the variability of the subsystem of atoms nk, which defines all atoms with resultant forces exceeding the minimum allowable. The smaller this subsystem is the faster the convergence. Since atoms are moved into new positions independently, the most efficient process of self-adjustment is the one producing at each iteration the smallest subsets nk. In Fig. 8 the variability of the subsets nk 2 Sm during the convergence to equilibrium for the 6307-atoms lattice during the first step of loading is shown. Note that the total number of atoms in the moving set was Sm = 6045. As one can see, the LCAD and LEM methods produce qualitatively different results as those for the ECAD one. The first two cause rapid growth and then decay of the subsystem of atoms with resultant forces exceeding convergence criterion, whereas the ECAD method leads to more steady change and much smaller subsystems. In fact, for the ECAD method the number of atoms in nk is 10-fold smaller than the total number of atoms in the relaxation subset Sm. This explains the better efficiency of the ECAD method in Fig. 7 compared to the two other alternatives. This property of the ECAD method, in addition to its ability to utilize parallel processing, makes it potentially more advantageous compared to the classical minimization algorithms. 4. Conclusion The idea of a pseudo-atom, representing the collective effect of all neighbour atoms to a given one, was introduced and
investigated for three different numerical algorithms. The investigations were carried out for a system of LJ 2D crystals with different number of atoms and were limited to a tensile extension of these crystals. The results for a 275-atoms lattice show that all three algorithms produce identical stress–strain data and interatomic force fields, and capable of detecting dislocations formation. However, the algorithm introduced previously in [5] gives the best computational efficiency. This is explained by examining the variability of the subsystem of atoms with resultant forces exceeding the minimum allowable one during the system relaxation to equilibrium. It is shown that in a self-organized system the convergence to equilibrium by independent displacement of atoms is efficient but strongly depends on the magnitudes of atoms displacements. All three methods allow parallel processing, which can be utilized to further increase computational efficiency. The above numerical investigation was limited to a specific range of system sizes, to planar lattices, and a specific loading. Further applications of the method suggested in [5] to larger 3D systems will allow to clarify its potential for simulations of more general molecular statics problems. References [1] D.J. Wales, H.A. Scheraga, Science 285 (5432) (1999) 1368–1373. [2] A.R. Leach, Molecular Modelling: Principles and Applications, Prentice Hall, 2001. [3] I. Ye. Telitchev, O. Vinogradov, Int. J. Comp. Methods 3 (1) (2006) 71–81. [4] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Phil. Mag. 64 (4) (1990) 851–878. [5] O. Vinogradov, Int. J. Comp. Methods 3 (2) (2006) 153–161. [6] O. Vinogradov, Comp. Mater. Sci 39 (3) (2007) 611–615. [7] O. Vinogradov, Comp. Mater. Sci. 42 (2008) 478–482.