A deformation–fluctuation hybrid method for fast evaluation of elastic constants with many-body potentials

Computer Physics Communications 183 (2012) 261–265

Contents lists available at SciVerse ScienceDirect

Computer Physics Communications www.elsevier.com/locate/cpc

A deformation–fluctuation hybrid method for fast evaluation of elastic constants with many-body potentials Yubao Zhen ∗ , Chengbiao Chu School of Astronautics, Harbin Institute of Technology, Harbin 150001, People’s Republic of China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 26 July 2011 Received in revised form 28 September 2011 Accepted 29 September 2011 Available online 5 October 2011 Keywords: Elastic constants Fluctuation method Many-body potentials Group perturbation

A universal, fast deformation–fluctuation hybrid approach with the merits of both the direct and fluctuation methods for evaluation of elastic constants is proposed. Deformation is shown to provide an inherent connection between the latter two methods. Based on this finding, the hybrid method utilizes group perturbation of atoms by means of six homogeneous deformations for fast evaluation of the Born term in the stress fluctuation method. The hybrid method is tested with three typical and widely used many-body potentials, i.e., for crystalline silicon with a Stillinger–Weber potential, for crystalline copper with second nearest-neighbor modified embedded-atom method and for diamond with a secondgeneration reactive empirical bond order potential. The calculated elastic constants agree well with either the theoretical or experimental results for all three cases. The hybrid method is especially valuable for complicated many-body potentials for which the analytical second derivative of the energy is challenging or impractical to obtain. © 2011 Elsevier B.V. All rights reserved.

1. Introduction In the calculation of elastic constants in atomistic models, the direct method [1] and the fluctuation method [1–4] are commonly regarded as two fundamentally different approaches. The current work proves in theory that they are inherently connected, and proposes a universal and fast algorithm for elastic constant calculations. The direct method differs from the fluctuation method in every aspect. The theory behind the direct method is the linear relation between stress and strain, while the fluctuation method is based on the statistical relation of elastic constants to fluctuation of certain thermodynamic quantities. Operationally, the direct method generates deformations via strain-control or stress-control and requires multiple deformation modes to obtain all of the elastic constants. The fluctuation method, on the other hand, needs only one round of calculation for all constants and has several variants, namely strain fluctuation [2], stress–strain fluctuation [3] and the so-called stress fluctuation [5,6]. Some other recently proposed methods [7,8] are all based on fluctuation methods. With respect to advantages and disadvantages, the direct method is simple both in concept and implementation but suffers from fluctuations in stress or strains at non-zero temperatures. Fluctuation methods are more efficient but also suffer in some degree from convergence difficulties; see, for example, the strain fluctuation papers [1,4].

*

Corresponding author. E-mail address: [email protected] (Y. Zhen).

0010-4655/$ – see front matter doi:10.1016/j.cpc.2011.09.006

©

2011 Elsevier B.V. All rights reserved.

Among all fluctuation methods, the stress fluctuation in EhN or ThN ensembles proposed by Ray and Rahman [5,6] features much more rapid convergence than others and provides unique insight into the contributors to the elastic constants. Here E is the system energy, h = (a, b, c) is a 3 × 3 matrix constructed from the three vectors a, b and c forming the molecular dynamics cell, N is the particle number, and T is the temperature. The distinct significance of the stress fluctuation method lies in the fact that the major contributor to the elastic constant is the Born term [4] (term with the second energy derivative), which is not a fluctuation expression and converges to an accurate value very quickly. The kinetic and fluctuation terms give only minor contributions and even the latter generally take longer to converge. However, the stress fluctuation approach requires the second energy derivative to be evaluated, which is challenging for most modern many-body potentials. The current work aims to eliminate this obstacle by linking the deformation to fluctuation and proposing a hybrid, universal and fast numerical approach for evaluation of the Born term. This work is inspired by the work of Lutsko [9] and Yoshimoto et al. [10]. This paper is organized as follows. In Section 2, the theory of the deformation–fluctuation hybrid method is proposed and general information about simulation is given. In Section 3, results of the hybrid method applied to Stillinger–Weber potential, second nearest-neighbor modified embedded-atom method, and secondgeneration reactive empirical bond order potential are presented. A discussion then follows in Section 4. Throughout this paper, Latin indices refer to Cartesian components and Greek indices refer to particle labels. Also, the

262

Y. Zhen, C. Chu / Computer Physics Communications 183 (2012) 261–265

summation convention of repeated Latin indices applies unless otherwise stated. 2. Methodology 2.1. Theory Ray and Rahman [5] derived a common formalism for adiabatic (in EhN) and isothermal (in ThN) elastic constants C pqmn of a general anisotropic material under applied stress 1 −1 −1 −1 V 0 h− 0ip h 0 jq h 0km h 0ln C pqmn

4 

1 k(α , β, η, ξ ) = α β ηξ r r



 ∂ U δαη δβξ ∂2U . − ∂ r α β ∂ r ηξ ∂ r αβ r αβ

αβ

(1)

(2)

α and β are denoted as riα β = riα − riβ ,

s i = sα − si and distance as r α β = (r i r i )1/2 . The potential U i depends only on interatomic distances. Efficient Born term calculation in stress fluctuation expressions from many-body potentials for which analytical second derivatives are unavailable or intractable is a topic of long-standing interest. Efficiency here means minimal computational demand and universal applicability to different potentials. Yoshimoto et al. [10] proposed an approximation method for Born term calculation based on the derivative of force with respect to atomic positions, where each atom is perturbed consecutively along three coordinate directions followed by a force recomputation. For system size of N, their method requires 3N system-wide force evaluations for one round of elastic constant calculation. In this work, a deformation– fluctuation hybrid method is proposed which reduces this value to six, resulting in a fast algorithm independent of system size in the above sense. The theory of the hybrid method starts with a relation derived by Yoshimoto et al. [10] for the Born term

1  V

α β α β ηξ ηξ

k(α , β, η, ξ )r i r j rk rl

α <β η<ξ

  v = δil σ jk + δik σ jlv + D i jkl ,

(3)

where

D i jkl = −

 α

σivj = −

r αj

∂ σklv , ∂ r iα

1  α ∂U 1  α α ri = r f . α V α ∂r j V α i j

α β α β ηξ ηξ

k(α , β, η, ξ )r p rq rm rn

α <β η<ξ

αβ αβ

β



   1 −1 −1 −1 v v + D pqmn . = V h− h jq hkm hln δ pn σqm + δ pm σqn ip

The three terms on the right hand side of Eq. (1) correspond to the fluctuation, kinetic and Born term, respectively. Here h0 is the reference state with zero strain of the tensor h, k B is the Boltzmann constant, T is the temperature, V is the volume and the bracket   represents the ensemble average. M = − V h−1 σ h− T /2 is the thermodynamic stress tensor, σ is the microscopic stress tensor and G = h T h is the metric tensor. The superscript ‘T’ indicates the matrix transpose and δβξ is the Kronecker delta. For atom α , its position in real space and scaled space are denoted −1 α as r iα and sα i = h ip r p , respectively. The interatomic coordinate differences between atoms

α β α β ηξ ηξ

k(α , β, η, ξ )si s j sk sl

1 −1 −1 −1 = h− h jq hkm hln ip

α <β,η<ξ

where



α <β,η<ξ



 M i j Mkl  −  M i j  Mkl  kB T  1 −1  1 −1 + 2Nk B T G − G jk + G − G jl il ik   α β α β ηξ ηξ  + k(α , β, η, ξ )si s j sk sl ,

=−

Note σivj is the symmetric virial stress tensor, and the aforementioned microscopic stress tensor σi j = σivj + Nk B T δi j / V . Eq. (3) links two differentiations, where the left hand side is about pair length and the right hand side is about atom position. The latter is computationally more favorable since for system size of N, the number of unique pairs is N ( N − 1)/2, ( N − 1)/2 times the number of atoms. αβ 1 αβ To comply with the general form of Eq. (1), with si = h− ip r p and Eq. (3), the Born term is rewritten in terms of D i jkl and σivj as

(4) (5)

(6)

To simplify notation and also for later use, here we define the tenv v sor B pqmn = δ pn σqm + δ pm σqn + D pqmn . Eq. (6) shows that the major computational demand is from D i jkl . Unlike the atom-by-atom perturbation applied in Ref. [10], our strategy here is to perturb all the atoms simultaneously in a total of six special modes. These modes are derived as follows. Suppose at an instant, for two given directions i , j, all atoms undergo an arbitrary system-wide perturbation dr iα along i. To first order of approximation, the virial stress change due to such a perturbation is

dσklv |i -ptb =

 ∂σ v kl

∂ r iα α

dr iα

(no sum over i ),

(7)

where subindex i-ptb denotes perturbation along i. If such a system perturbation is taken as a special mode

dr iα = e i j r αj

(no sum over j ),

(8)

where e i j is a small constant common for all atoms specific to i j, Eq. (7) reads

dσklv |i -ptb = e i j

 ∂σ v kl

α

∂ r iα

r αj

(no sum over i , j ).

(9)

By definition in Eq. (4), this leads to

D i jkl = −

dσklv |i -ptb ei j

= − lim

e i j →0

σklv |i-ptb − σklv0 ei j

,

(10)

where σklv0 is the initial value of the virial stress tensor. Note Eq. (10) should be evaluated at each time step and the value of the tensor D so obtained is used to calculate the tensor B previously defined. The underlying physics of Eq. (8) becomes obvious upon noticing that dr iα is in fact the atomic displacement u α . It follows then i α , which is the conventional atomic deformation that e i j = ∂ u α /∂ r i j gradient at atom α . Since for each i j, e i j is assumed to be constant for all atoms, the six perturbation modes, i.e., i j as 11, 22, 33, 23, 13, 12, correspond to three uniform uniaxial deformations along coordinate directions and three two-dimensional in-plane uniform shear deformations. These deformations are similar to those applied in direct method and are connected to the Born term in EhN/ThN fluctuation method via the virial part of the deformation stress as shown by Eqs. (10) and (6). The deformation and fluctuation are thus naturally integrated, resulting in a hybrid method. The hybrid method is computationally efficient and universally applicable to all potentials with only mutual atomic distance dependence. This can be seen via Eqs. (8) and (10), where for each of the six deformations, σivj , hence D i jkl and the Born term, all

Y. Zhen, C. Chu / Computer Physics Communications 183 (2012) 261–265

can be evaluated with only forces on atoms and no other energy derivatives are required. Also, the form of Eq. (10) implies that any numerical finite difference scheme for the first derivative can be applied, though with different precision and corresponding computational demands. Eqs. (6), (8) and (10) are the major theoretical results in this work. To end this section, the way to calculate the elastic constant tensor C from Eq. (1) can be described briefly as follows. 1 −1 −1 −1 h jq and H 0i jpq ≡ h0ip h0 jq . Due to permutaDefine H i jpq ≡ h− ip tion symmetry of matrices C , B, M, G and H (note that H jiqp = H i jpq , and similarly for H 0 ), the 9 × 9 matrix manipulation in Eq. (1) can be converted to the conventional 6 × 6 form

˜ 0 C H˜ 0T =  F  +  K  + V H˜  B  H˜ T , V0H

(11)

where F , K and B stand for fluctuation, kinetic and Born term, ˜ and H˜ 0 share the same respectively. The 6 × 6 components of H ˜ as an example, for Voigt indices ρ = (i j ) and μ = form. Taking H ( pq), this is

˜ ρμ = H



H i jpq , H i jpq + H i jqp ,

μ  3, μ > 3.

(12)

Eq. (11) can be inverted to evaluate the usual 6 × 6 elastic constant matrix C i j . 2.2. Simulations In this work, the hybrid method is applied to three typical and widely used many-body potentials that are very different in nature, namely the Stillinger–Weber [11] (SW) potential, second nearestneighbor modified embedded-atom method [12,13] (2NN MEAM) and the second-generation reactive empirical bond order potential (REBO) [14]. We refer the reader to the original literature for details of the respective formalisms. Each potential is chosen for a reason. For SW, the results of Kluge et al. [15] for a periodic system with 216 silicon atoms using EhN dynamics based on analytical potential derivatives has become a classic standard. The 2NN MEAM is widely used for metallic system, but lacks a published second energy derivative (see Ref. [16] for the sophisticated analytical first derivative). The second-generation REBO allows for covalent bond breaking and forming, and depends explicitly on bond and dihedral angles. The analytical potential derivatives of REBO are thus very complicated [17]. Even so, REBO still depends on interatomic distances, hence the stress-fluctuation and the proposed hybrid method are applicable. For each potential, a specific crystalline periodic system is first equilibrated to zero pressure and desired temperatures with a multi-stage baro- and temperature-control process. The relative standard deviations for h0 are all less than 1 × 10−4 and the average microscopic pressure tensor in the system is on the order of tens of MPa, which is a negligible stress for solid crystalline systems considered here. For all calculations, Born terms and corresponding adiabatic elastic constants are calculated in the EhN ensemble with molecular dynamics. It should be noted though the hybrid method proposed here applies equally to the ThN ensemble for isothermal elastic constants. For each case, only the first 5000 steps of the EhN dynamics are turned on for Born term calculation because of its rapid convergence. The total run, however, is much longer simply in order to obtain steady fluctuation terms. Note that at zero temperature, the static Born terms are calculated in just one step. The symmetry averaged individual contributions to elastic constants from Eq. (1) are also calculated. Standard deviations of symmetry equivalent values are taken as the error estimates.

263

Table 1 The three contributions to the adiabatic elastic constants of silicon from Eq. (1) in units of 1011 dyne/cm2 for three different temperatures. Also shown are the static Born values at 0 K and the two- and three-body contributions to total Born terms. For 1477 K, individual contributions from benchmark tests are listed as well. C 11 T = 0 K static Born terms This work 15.14 Theoretical [15] 15.10

C 12

C 44

7.64 7.60

10.97 10.90

T = 888 K 2-body Born 3-body Born Born total

9.49 4.80 14.29

9.39 −1.87 7.52

9.39 1.03 10.42

Kinetic Fluctuation

0.24 −0.34

0.00 0.01

0.12 −4.92

This work Theoretical [15] Experiment [19]

14.19 ± 0.02 14.14 ± 0.01 15.75

7.53 ± 0.01 7.52 ± 0.00 6.05

T = 1164 K 2-body Born 3-body Born Born total

8.84 5.09 13.93

8.89 −1.42 7.47

8.89 1.29 10.18

Kinetic Fluctuation

0.33 −0.53

0.00 −0.02

0.16 −5.66

This work Theoretical [15] Experiment [19]

13.73 ± 0.02 13.73 ± 0.03 15.25

T = 1477 K 2-body Born 3-body Born Born total

8.09 5.57 13.66

8.38 −0.90 7.48

8.38 1.60 9.98

Kinetic Fluctuation

0.40 −0.71

0.00 −0.05

0.20 −5.87

This work

13.35 ± 0.02

2-body Born 3-body Born Born total

8.00 5.65 13.65

8.34 −0.89 7.45

8.34 1.64 9.98

Kinetic Fluctuation

0.40 −0.73

0.00 −0.06

0.20 −5.98

Theoretical [15] Experiment [19]

13.32 ± 0.01 14.80

7.45 ± 0.01 7.43 ± 0.01 5.90

7.43 ± 0.01

7.39 ± 0.04 5.75

5.62 ± 0.54 5.24 ± 0.84 7.53

4.68 ± 0.21 4.57 ± 1.14 7.32

4.31 ± 0.36

4.20 ± 0.83 7.00

All simulations are conducted using LAMMPS [18], with our new hybrid algorithms incorporated. For the group perturbations of all atoms in the system, the usual second order central difference scheme is applied. And for simplicity, all e i j are taken to have a common value 1 × 10−5 . Temperature fluctuations are within 3% about the desired temperatures. Other specific parameters for each potential are given in Section 3. 3. Results 3.1. SW The hybrid method is tested first against the benchmark results in Ref. [15] with the SW potential for silicon [11]. All parameters from the literature are kept except for the number of total steps for EhN dynamics. In Ref. [15], a total of 150 000 steps was performed. This proved not to be long enough to have the shear elastic constants C 44 , C 55 and C 66 nearly converged to a common value. In this work the total number of steps is increased to 1 million in order to accomplish the goal. The elastic constants calculated by the hybrid method along with the theoretical results [15] and from experiment [19] are presented in Table 1. Note the two- and three-body Born terms are obtained by replacing the virial stress tensor in Eq. (10) with its two- and three-body contributions, respectively. One sees that the

264

Y. Zhen, C. Chu / Computer Physics Communications 183 (2012) 261–265

Table 2 The three contributions to the adiabatic elastic constants of copper from Eq. (1) in units of GPa at zero pressure with the 2NN MEAM. Also shown are the static Born terms at 0 K. T (K)

B C 11

K C 11

F C 11

C 11 total

B C 12

K C 12

F C 12

C 12 total

B C 44

K C 44

F C 44

C 44 total

0 300 350 400 450 500 550 600 650 700 750 800

176.16 172.07 171.26 171.08 170.67 170.44 169.78 169.15 168.21 167.55 166.91 166.80

0.0 1.36 1.54 1.82 2.04 2.32 2.53 2.73 2.90 3.07 3.30 3.62

0 .0 −4.18 −4.64 −5.35 −5.69 −6.55 −6.88 −7.43 −7.90 −8.22 −8.97 −9.34

176.16 169.25 ± 0.17 168.16 ± 0.67 167.55 ± 0.09 167.02 ± 0.38 166.21 ± 0.38 165.43 ± 0.29 164.45 ± 0.07 163.21 ± 0.50 162.40 ± 0.14 161.24 ± 0.32 161.08 ± 0.44

124.92 118.28 117.09 116.06 114.93 113.93 112.79 111.79 110.68 109.87 108.66 107.80

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0 .0 −1.13 −1.27 −1.44 −1.40 −1.59 −1.67 −1.83 −1.88 −1.97 −1.80 −1.84

124.92 117.15 ± 0.11 115.82 ± 0.37 114.62 ± 0.21 113.53 ± 0.14 112.34 ± 0.03 111.12 ± 0.17 109.96 ± 0.08 108.80 ± 0.16 107.90 ± 0.10 106.86 ± 0.33 105.95 ± 0.11

81.80 81.42 81.43 81.93 82.23 82.74 82.89 82.95 82.78 82.83 82.97 83.51

0.0 0.68 0.77 0.91 1.02 1.16 1.27 1.36 1.45 1.54 1.65 1.81

0 .0 −4.40 −4.94 −5.61 −6.11 −6.93 −7.27 −7.73 −8.16 −8.30 −9.08 −9.94

81.80 77.70 ± 0.08 77.26 ± 0.27 77.23 ± 0.09 77.14 ± 0.12 76.97 ± 0.20 76.89 ± 0.22 76.58 ± 0.14 76.07 ± 0.14 76.07 ± 0.09 75.54 ± 0.70 75.38 ± 0.20

results from the hybrid method match those from the benchmark tests on C 11 and C 12 very well. Also, from error estimates, the hybrid results are shown to be more nearly converged to a common value for the shear elastic constants C 44 , C 55 , and C 66 . The significant deviations between the calculated and experimental results are due to the quality of the SW potential for silicon. As for the large fluctuations for C 44 presented in Table 1, this behavior was already discussed in Ref. [15] and was identified to be caused by the internal strain in silicon. 3.2. 2NN MEAM The hybrid method next was tested with the 2NN MEAM for a periodic crystalline copper system of 500 particles in a 5 × 5 × 5 FCC lattice. The 2NN MEAM parameters for Cu are taken from Lee et al.’s work [13] without any modification. The radial cutoff distance was 4.0 Å and time step was 1.0 fs. Eleven different temperatures from 0 K to 800 K are considered. For each temperature, a total of 500 000 EhN molecular dynamics steps was carried out. The results are given in Table 2, with fluctuation, kinetic and Born terms all listed. Also, see Fig. 1 for a comparison against experimental results for all temperatures. The static Born terms at 0 K coincide with the elastic constants at 0 K. From Fig. 1, the calculated elastic constants at 0 K match the experimental results. However, the deviations between calculated and experimental results increase gradually with temperature, especially for C 44 . At 800 K, the relative discrepancies for C 11 and C 12 are less than 7%, while for C 44 it is around 20%. The discrepancies originate from the fitting process for the potential parameters, in which the elastic constants at 0 K are used. The overall quality of 2NN MEAM parameters for Cu on elastic constants at finite temperatures is reasonably good. For non-zero applied stress, a case at 300 K with σxx = 5 GPa is tested. The resultant configuration at 300 K and zero stress is taken as the reference state. The sample then undergoes a uniaxial tension along x direction till σxx = 5 GPa, which corresponds approximately to a strain of 6.05%. Again, a total of 500 000 steps is performed for EhN dynamics. The results are shown in Table 3. Lattice parameters before and after loading are a = 3.6325 Å, b = 3.6328 Å, c = 3.6329 Å and a = 3.8530 Å, b = 3.5521 Å, c = 3.5522 Å, respectively. Obviously, under 5 GPa tensile stress, the Table 3 Adiabatic elastic constants and their three components of copper at 300 K and pressure

Fig. 1. Adiabatic elastic constants of copper at different temperatures from 0 K to 800 K with the 2NN MEAM. Hybrid results (dotted lines with crosses) are compared against experimental results (dashed lines with squares). The latter are taken from Ref. [20] for 0 K and Ref. [21] for finite temperatures.

lattice has deformed into a tetragonal one and the elastic constants show the tetragonal symmetry. 3.3. Second-generation REBO For the second-generation REBO potential, a periodic diamond system with 1000 particles is tested. Gao et al. [22] have studied the elastic constants of such a system at different temperatures using the analytical energy derivatives developed in Ref. [17]. To reproduce the theoretical literature results, temperatures and relevant parameters all were taken from Ref. [22]. The time step was 0.25 fs. Results for REBO by hybrid methods are listed in Table 4 and depicted as well in Fig. 2 along with the theoretical ones. The two results match very well. 4. Discussions The deformation–fluctuation hybrid method inherits the merits of both the direct method and the fluctuation method. The deforσxx = 5 GPa in units of GPa with the 2NN MEAM.

Term

C 11

[C 12 , C 13 ]

[C 22 , C 33 ]

C 23

C 44

[C 55 , C 66 ]

Born total Kinetic Fluctuation

129.40 1.35 −3.92

95.31 0.00 −0.77

183.24 1.35 −4.15

139.69 0.00 −1.86

125.22 0.68 −2.99

56.90 0.68 −4.99

Total

126.83

94.54 ± 0.19

180.44 ± 0.29

137.83

122.91

52.59 ± 3.01

Y. Zhen, C. Chu / Computer Physics Communications 183 (2012) 261–265

265

Table 4 The three contributions to the adiabatic elastic constants of diamond from Eq. (1) in units of GPa at zero pressure with the second-generation REBO. Also shown are the static Born terms at 0 K. T (K)

B C 11

K C 11

F C 11

C 11 total

B C 12

K C 12

F C 12

C 12 total

B C 44

K C 44

F C 44

C 44 total

0 300 500 700 900 1100

1075.60 973.34 939.86 909.84 886.68 867.33

0 .0 2.94 4.61 6.60 8.60 10.65

0.0 −4.42 −7.07 −10.44 −14.73 −19.72

1075.60 971.86 ± 0.17 937.40 ± 0.11 906.00 ± 0.45 880.55 ± 0.29 858.26 ± 0.15

125.44 165.39 174.01 183.87 196.25 212.98

0.0 0.0 0.0 0.0 0.0 0.0

0 .0 −0.24 −0.46 −0.80 −1.29 −3.51

125.44 165.15 ± 0.28 173.55 ± 0.25 183.07 ± 0.59 194.96 ± 0.03 209.47 ± 0.62

738.02 710.87 698.02 687.26 680.74 677.77

0 .0 1.47 2.30 3.30 4.30 5.32

0 .0 −33.05 −45.54 −51.07 −58.32 −70.79

738.02 679.29 ± 6.55 654.78 ± 2.63 639.49 ± 5.15 626.72 ± 5.65 612.30 ± 1.25

tentials. The hybrid method requires only force information and computationally is highly efficient, featuring only six system-wide force calculation for each round of evaluation. Applications of the hybrid method to SW, 2NN MEAM and second-generation REBO potentials show very good accuracy. The proposed hybrid method is especially useful for many-body potentials for which the explicit analytical second derivative is impractical to obtain or hard to implement. Acknowledgements The authors would like to thank the anonymous reviewer for constructive comments and Specialist Editor S.B. Trickey for kind help on language polishing. One of the authors (Y.Z.) is grateful for helpful discussions with Dr. Zhiwei Cui. This material is based upon work supported by the National Natural Science Foundation of China under Grant No. 10802025. Fig. 2. Born terms and adiabatic elastic constants of diamond at six temperatures from 0 K to 1100 K with the second-generation REBO. Hybrid results (dotted line with cross) are compared against the theoretical results (dashed line with square, taken from Ref. [22]).

mation is interpreted as group perturbations of atom positions in only six special modes. The key barrier in the stress fluctuation method on evaluation of Born term is removed by a simple finite difference approximation of the first derivative of virial stress to deformation gradient. Only forces on atoms are involved in the calculations. This character enables both computational efficiency and universal applicability of the hybrid method to all potentials depending only on interatomic distances. The hybrid method is rooted in the stress fluctuation method, therefore it converges rapidly. It is also as easy to implement as the direct method since atomic displacements are explicitly specified in the theory. Note that for some special cases, simple explicit analytical expressions for the second derivative of the potential exist. Applying the hybrid method with such potentials might introduce some extra cost which can be avoided by utilizing the analytical derivative. 5. Conclusions This work proposes a universal and fast deformation–fluctuation hybrid method to calculate the elastic constants in EhN/ThN ensemble of molecular dynamics with sophisticated many-body po-

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

M. Sprik, R.W. Impey, M.L. Klein, Phys. Rev. B 29 (1984) 4368–4374. M. Parrinello, A. Rahman, J. Chem. Phys. 76 (1982) 2662–2666. A.A. Gusev, M.M. Zehnder, U.W. Suter, Phys. Rev. B 54 (1996) 1–4. J.R. Ray, Comput. Phys. Reports 8 (1988) 109–152. J.R. Ray, A. Rahman, J. Chem. Phys. 80 (1984) 4423–4428. J.R. Ray, A. Rahman, J. Chem. Phys. 82 (1985) 4243–4247. Z. Cui, Y. Sun, J. Li, J. Qu, Phys. Rev. B 75 (2007) 214101. J. Li, Y. Sun, Z. Cui, F. Zeng, Comput. Phys. Comm. 182 (2011) 1447–1451. J.F. Lutsko, J. Appl. Phys. 65 (1988) 2991–2997. K. Yoshimoto, G.J. Papakonstantopoulos, J.F. Lutsko, J.J. de Pablo, Phys. Rev. B 71 (2005) 184108. F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262–5271. B.J. Lee, M.I. Baskes, Phys. Rev. B 62 (2000) 8564–8567. B.J. Lee, J.H. Shim, M.I. Baskes, Phys. Rev. B 68 (2003) 144112. D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, S.B. Sinnott, J. Phys.: Condens. Matter 14 (2002) 783–802. M.D. Kluge, J.R. Ray, A. Rahman, J. Chem. Phys. 85 (1986) 4028–4031. P.M. Gullet, G. Wagner, A. Slepoy, Numerical tools for atomistic simulations, SANDIA Report 2003-8782, Sandia National Lab, 2003. K.V. Workumand, G.T. Gao, J.D. Schall, J.A. Harrison, J. Chem. Phys. 125 (2006) 144506. S.J. Plimpton, J. Comp. Phys. 117 (1995) 1–19, http://lammps.sandia.gov. K.-H. Hellwege, A. Hellwege (Eds.), Landolt–Börnstein: Crystal and Solid State Physics, vol. 11, Springer, Berlin, p. 116. G. Simmons, H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook, 2nd edition, MIT Press, Cambridge, MA, 1971. Y.A. Chang, L. Himmel, J. Appl. Phys. 37 (1966) 3567–3572. G.T. Gao, K.V. Workum, J.D. Schall, J.A. Harrison, J. Phys.: Condens. Matter 18 (2006) S1737–S1750.