- Email: [email protected]
Quantum effects on elastic constants of diamond by path-integral Monte Carlo simulations
Computational Materials Science xxx (xxxx) xxxx
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Quantum effects on elastic constants of diamond by path-integral Monte Carlo simulations B.G.A. Britoa, G.-Q. Haib, L. Cândidoc,
⁎
a
Departamento de Física, Instituto de Ciências Exatas e Naturais e Educação (ICENE), Universidade Federal do Triângulo Mineiro – UFTM, 38064-200 Uberaba, MG, Brazil b Instituto de Física de São Carlos, Universidade de São Paulo, 13560-970 São Carlos, SP, Brazil c Instituto de Física, Universidade Federal de Goiás, 74690-900 Goiânia, GO, Brazil
ARTICLE INFO
ABSTRACT
Keywords: Diamond Quantum effects on elastic constants Path integral Monte Carlo
Using the path-integral Monte Carlo method, we investigate the quantum effects on the elastic constants C11, C12 , and C44 of the diamond crystal in a wide temperature range at ambient pressure. The Tersoff potential is used to describe the interatomic interactions and the elastic constants are determined by a direct method derived from the stress-strain curve. We find that the elastic constants C11 and C44 behave as a quadratic function of temperature at low temperatures and C12 is practically a constant being independent of temperature. The quantum effects are significant at low temperatures up to about 1000 K. The quantitative differences between the PIMC calculations and the experimental measurements of the elastic constants C11, C12 , and C44 at low temperatures are about 3%, 17% , and 9%, respectively. Using the obtained elastic constants, we also estimate the bulk modulus B, anisotropic factor A, Cauchy pressure PCauchy , and Poisson ratio of the diamond crystal.
1. Introduction Diamond is the hardest substance known to man. It has excellent transparency over an extensive wavelength range, high thermal conductivity at room temperature, and very high resistance to chemical corrosion [1,2]. For example, in the study of high pressure materials, the hardness property has been an important characteristic used for diamond anvil cell. Both the diamond and graphite are composed of the single chemical element: carbon. The difference between these two solids is the arrangement of the carbon atoms in the lattices and, at ambient conditions, graphite is thermodynamically the most stable allotrope of carbon by a small enthalpy difference of 2.9 kJ mol−1 [3]. However, there is an almost prohibitive activation barrier between them making their interconversion impossible at ambient conditions [4]. The mechanical and electronic properties of diamond have been widely investigated [5,6]. There are studies showing that the quantum nature of atomic motion affects the thermal and mechanical properties of the material at room temperature [7,8]. Quantities such as vibrational energy, kinetic energy, heat capacity, and bulk modulus deviate significantly from the classical predictions even at ambient conditions. It is known that the methods such as density functional theory (DFT) and molecular dynamics (MD) are able to accurately predict several ⁎
properties of covalent solid, in particular, the lattice parameter, phonon dispersion relation [9], phase transition [10], and dynamics of point and extended defect formation [11]. However, within these theoretical frameworks, there are difficulties in dealing with the zero-point vibrational motion of the atoms and handling quantum effects at finite temperatures. Path-integral Monte Carlo (PIMC) is a powerful method that treats the atomic nuclei as quantum particles naturally. PIMC is versatile and can use different approximations such as the empirical potential [12–14], or those based on standard electronic structure methods [15,16]. It has been considered as one of the most important methods to simulate a quantum system at finite temperature with essential contributions to our understanding of the thermodynamic properties of several systems, mainly in the condensed matter [8,13,14,17–20]. In this paper, we apply the PIMC method to investigate the quantum effects on the elastic-stiffness coefficient of the diamond. Our simulation shows that the quantum effects impact the elastic constants and related quantities of the diamond considerably even at room temperature. To evidence the quantum effects in this system, we will compare the PIMC calculation results with those obtained from the classical Monte Carlo (MC) simulations and also with available experimental results. The paper is organized as follows. In Section 2 we present a short description of the PIMC formalism deriving an expression for the
Corresponding author. E-mail address: [email protected] (L. Cândido).
https://doi.org/10.1016/j.commatsci.2019.109387 Received 16 August 2019; Received in revised form 25 October 2019; Accepted 26 October 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Please cite this article as: B.G.A. Brito, G.-Q. Hai and L. Cândido, Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109387
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
stress-tensor observable. In Section 3 we present the results for elasticstiffness coefficients C11, C12 , and C44 as a function of temperature. Using the estimated elastic constants, we also obtain other quantities such as the bulk modulus, anisotropic factor, Cauchy pressure, and Poisson ratio. Finally, we summarize our work in Section 4.
The above equation, i.e. Eq. (7), is exact. But, in the calculations we use the action within the primitive approximation, given by
2. Numerical approach and computational details
Spri=
ij
We perform both the PIMC and classical MC simulations for a diamond crystal composed of N carbon atoms in the isothermal-isobaric (NPT) and isothermal-isochoric (NVT) ensembles in a wide temperature range from 50 to 2000 K and at pressure of 1 atm. In the PIMC simulations, the carbon atom is considered as a quantum particle of mass 12.011 a.u. Each atom is modeled by a set of M beads, forming a cyclic chain or a polymer. For the interatomic interactions, we employ the widely used Tersoff potential [21]. The specific parameters for the Tersoff potential used in this work can be found in the paper by Herrero and Ramírez in Ref. [8]. The PIMC method [22] is based on the evaluation of the quantum density matrix
Sk
d R1…d RM 1exp k=1
1
dRdR (R, R , ) R|O|R ,
dR 0 dRM
1
× exp( S (R 0 RM 1, R0 )).
ij
F (N , V , T ) =
lnZ (N , V , T ).
1 det(h )
3
hjk k=1
F hik
. N,T
.
N
M
=1
=1
(7)
|r
,
r
2 , + 1|
U (r , ),
(8)
NM
ij
mM 2
N
M
3
=1
=1
l=1
hil (s
hik (s
,
s , )k
M
N
M
=1
=1
,
V (h s , )i
s , )l, 3
hjk (s , )k
.
k=1
(9) Thus, the expectation value of the stress tensor components depends on the beads position difference, and of the first derivative of the potential related to the position of the beads [26]. There are two methods for estimating the elastic constants using computational simulations: (i) The direct method [27,26] which derives elastic constant from the stress-strain curve, and (ii) the fluctuation method [28,29] which use ensemble averages of the fluctuation of the strain or stress. Here, we use the direct method to obtain the elastic stiffness constants through the relation between stress and strain by the Hooke’s law
(2)
ij
=
Cijkl
kl,
(10)
k, l
(3)
where Cijkl represents the elastic stiffness tensor and kl is the kl component of the strain tensor. The PIMC simulations are performed in a three-dimensional box using a diamond unit cell with periodic boundary conditions in the three directions. Most calculations are performed within a simulation box containing 216 atoms. A “time step” = 1.333 × 10 4 K−1 (corresponding to M = 150 at T = 50 K) is used for the discretized imaginary time path integrals. The extrapolation to the = 0 limit shows that this approximation induces an error of about 0.1% in the obtained total energy. Checks for larger systems containing 512 atoms show that the obtained results do not change within the statistical error when comparing to a system with 216 atoms. For the averages, we use 200.000 PIMC time steps.
(4)
3. Results and discussions We start the simulation using the NPT ensemble to determine the system density as a function of temperature at a constant pressure of 1 atm. Then, the obtained densities are used as input in the NVT simulation to calculate the elastic constants. The NVT simulation cell uses its vectors parallel to the edges of the unit cell of the diamond lattice. If the simulation box is scaled or deformed by some strain, we have the new vectors of the cell given by the relation
(5)
matrix h in the following way [25],
=
=1
det(h )
(1)
We can write an estimator for the components of the stress tensor ( ) as a function of the derivative of the free energy F by the elements of the
ij
=1
k=1
where det(h ) gives the volume of the simulation cell. The connection of the path-integral formalism with the thermodynamic observables in NVT ensemble is made through the Helmholtz free energy (F) obtained by the partition function [24],
1
M
Mm 2 2
3
tween 0 and 1, and h is the matrix used in the Parrinello and Rahman approximation to describe the homogeneous distortions of the periodic simulation cell [23]. In the dimensionless coordinates, the partition function given by Eq. (3) can be re-written as
det(h ) NM dR 0 dRM 1× exp( S (R 0 RM 1, R 0)),
N
1
=
×
The bead position ri, k can be written as ri, k = h si, k , where si, k is the dimensionless coordinate whose components have values ranging be-
Z (N , V , T )=
k=1
where the first two terms are kinetic energy and the last potential energy. Eqs. (7) and (8) as shown above were first obtained by Ardila et al. [26]. By substituting Spri in Eq. (8) into Eq. (7) we get the thermodynamic estimator for the stress tensor components, given by
and the partition function is written as
Z (N , V , T ) =
M
S hik
ij
2 3NM ln 2 2Mm
+
by sampling the set of paths {R 0, R1, …RM 1, RM} with R k being the set of beads {r1, k , …, rN , k} , and ri, k a bead that is the position of the ith particle in the k th time slice. The “time step” is given by = / M with = 1/ kB T (kB is the Boltzmann constant and T the absolute temperature). The action of a link is defined as Sk S (R k 1, R k; ) ln[ (Rk 1, Rk ; )] is handled in the primitive approximation. To evaluate the 3NM -dimensional integrals in Eq. (1), we employ the Metropolis bisection sampling technique [18]. The expectation value of an operator O is given by
O =Z
NM
det(h )
M
(R 0 , RM ; ) =
3
1
=
(6)
Ai
It can be re-written in terms of the expectation value of the derivative of the action S by substituting the free energy given by Eq. (4) into Eq. (6), and taking the appropriate derivatives leading to
new
= (I ±
ij )
Ai ,
(11)
where Ai are the edge vectors and I is the unit matrix. Two strain tensors 2
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
10
(a)
1100
C11 (GPa)
σ xx (GPa)
5 0
(a)
1050
1000
−5 950
−10
0
−0.01
−0.005
0
0.005
400
0.01
800 1200 1600 Temperature (K)
εxx
(b)
160
6
C12 (GPa)
(b)
σ xy (GPa)
4 2
2000
140 120
0
100
−2
80 0
−4
400
800
1200
1600
2000
Temperature (K)
−6 −0.01
−0.005
0
0.005
(c)
0.01
640
εxy
xx
e 0 0 = 0 0 0 0 0 0
C44 (GPa)
Fig. 1. Responses of (a) the tensile stress xx under the imposed strain xx and (b) the shear stress xy under the imposed strain xy at temperature 300 K and pressure 1 atm.
(12)
0
0 e /2 0 = e /2 0 0 0 0 0
yy xy
= C11 xx , = zz = C12 = C44 xy.
400
800
1200
1600
2000
Temperature (K) Fig. 2. Temperature dependence of the elastic constants (a) C11 , (b) C12 , and (c) C44 of diamond at pressure 1 atm. The red open and green filled circles are results from the PIMC and classical MC simulations, respectively, with the red and green dotted lines being guides for the eye. Different experimental measurements are given by the black stars [30], black triangles [31], black squares [32], and the black solid curves with the error range indicated by two black dotted lines [33]. The red solid curves are the quadratic fittings of the PIMC results of C11 and C44 at low temperatures in the range 0 to 1600 K and an average value of C12 shifted upward by 34.8 GPa in (a) and 21.5 GPa in (b) and shifted downward by 52.1 GPa in (c).
(13)
are used to determine the elastic constants by the relations: xx
560
520
and xy
600
xx ,
(14)
The values of ij are obtained from the ensemble average in the NVT simulation with e = 0, ± 0,005, and ± 0.01 in Eqs. (12) and (13). The elastic constants are calculated from linear strain-stress relations given in Eq. (14). Fig. 1(a) and (b) give the obtained stresses xx and xy , respectively, as a function of applied strains at temperature 300 K and pressure 1 atm. In both cases, the stress shows a linear dependence on the strain for small distortions. Fig. 1(a) shows the tensile stress response xx as a function of the imposed tensile strain xx . For xx = 0 , we have only the hydrostatic compression, which in this case is small and falls within the error bar. The slope of the line gives the elastic constant C11 for longitudinal compression. Fig. 1(b) shows a similar response curve for shear stress xy in response of the shear strain xy . As expected for this kind of strain, the linear fit passes through the origin indicating no shear stress
at 1 atm, in other words, zero shear deformation. We verified that, in this range of strain, the stress-strain dependence is almost perfect linear for all the studied temperatures. When a tensile stress is applied along the x-axis of the diamond crystal, using Eq. (14) we can estimate the elastic constants C11 and C12 which are associated with the stress responses in the directions parallel with (longitudinal compression) and perpendicular to (transverse expansion) the x-axis, respectively. While C44 is shear modulus given by xy / xy . Fig. 2 shows C11, C12 , and C44 as a function of temperature at constant pressure 1 atm. The results obtained from the PIMC and classical MC calculations are represented by red open and green filled 3
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
circles, respectively. The experimental results from Migliori et al. [30], McSckimin and Adreatch [31], and Nagakubo et al. [32] are shown by black stars, triangles, and squares, respectively. The black solid curves give the experimental results from Zouboulis et al. [33] with two black dotted curves indicating the estimate experimental error range. Comparing C11 in Fig. 2(a) with C12 in 2(b), one can notice that the quantum effects on C11 are much more significant than on C12 . The classical and quantum results of C12 are close to each other with overlaps of the error bars in the whole temperature range. The PIMC results of C12 are practically independent of temperature with an average value of 105.3 GPa with error bars varying from 2 to 6 GPa. An upward shift of this average value by 21.5 GPa is shown by the red solid line in Fig. 2(b) which is in agreement with the experimental results of Zouboulis et al. [33]. However, the temperature dependence of C11 is significant. The classical MC predicts a linear decrease of C11 with increasing temperature while the PIMC results show a quadratic function dependence at low temperatures. At high temperatures, the PIMC and classical MC results of C11 approach to each other as expected. A numerical fitting of the PIMC results in the temperature range 0 < T < 1600 K yields (0) (2) 2 (0) (2) C11 c11 c11 T with c11 = 1048.6 GPa and c11 = 33.5 × 10 6 GPa/K2. Its upward shift by 34.8 GPa is given by the red solid curve in Fig. 2(a) showing a good agreement with the available experimental results. Fig. 2(c) shows the elastic constant C44 as a function of temperature. C44 is also known as shear modulus µ . The shear modulus describes the tendency of the object suffering deformation of its shape and keeping constant volume. Being similar to C11, the classical MC simulation shows that C44 decreases linearly with increasing temperature. The quantum effects on C44 are more significant at lower temperature. For T < 1600 K, the PIMC values of C44 are well described by the quadratic (0) (2) 2 (0) C44 c44 c44 T , c44 = 630.3 GPa function with and (2) c44 = 18.1 × 10 6 GPa/K2. Its downward shift by 52.1 GPa is in very good agreement with the experimental results as shown by the red solid curve in Fig. 2(c). In general, our PIMC simulations using the Tersoff’s potential demonstrate that the theory can describe correctly the temperature dependence of the elastic constants C11, C12 , and C44 . The domination of the quantum effects at low temperatures leads to the elastic constants C11 and C44 decrease quadratically with increasing temperature. Quantitatively, considering the experimental results as references, the PIMC calculations underestimate C11 and C12 at low temperatures by about 3% and 17%, respectively, and overestimate C44 by 9% . The anisotropy factor for a cubic crystal can be defined [34] by the 2C ratio A = C 44C . Because C44 = (C11 C12 )/2 for a cubic isotropic 11 12 crystal, the value of A is expected to be equal to 1 when the crystal is isotropic, while any deviation from the unity indicates a certain degree of elastic anisotropy. There have been studies showing that elastic anisotropy of crystal is important for material engineering due to its correlation with the producing of microcracks in materials [35]. For the temperature range under investigation, we find the value of A from the PIMC calculation is practically constant. The average value is 1.35 ± 0.07 showing that diamond is an anisotropic crystal even at ambient pressure. However, the value of A also indicates a negligible influence of the temperature on the elastic anisotropy of the diamond crystal. From the above calculations, we can estimate the Cauchy pressure for the system, giving by PCauchy = C12 C44 . A negative (positive) Cauchy pressure reveals brittleness (ductility) of the crystal [36,37]. The more negative PCauchy , the more angular character of atomic bonding in the crystal. Our PIMC and MC simulations using the Tersoff potential predict correctly a very negative Cauchy pressure of the diamond for all temperatures as shown in Fig. 3. This is a consequence of the correct directional bonding description for the carbon atoms in diamond lattice given by the Tersoff potential. In fact, diamond is the hardest known of the highest strength material and at the same time, the most brittle one of the most negative Cauchy pressure [37]. The
−460
P Cauchy (GPa)
−480
−500
−520
−540
−560 0
400
800
1200
1600
2000
Temperature (K) Fig. 3. The Cauchy pressure of diamond as a function of temperature. The red open and green filled circles are the results from the PIMC and classical MC simulations, respectively.
same negative behavior cannot be observed using pure pair potentials because they have no preferred direction of bonding. We also notice that the PCauchy value obtained from the classical MC approaches to the PIMC result at high temperatures. However, the classical approximation yields more negative Cauchy pressure for temperatures below 1000 K, indicating a more directional bonding. This is a consequence of the reduction of the atomic vibrational energy in the classical approximation where only thermal energy is considered. The PIMC simulation takes into account the quantum effects due to the zero-point energy of the atomic motion and predicts correctly the weak temperature dependence of the Cauchy pressure at low temperatures. The bulk modulus (B) for cubic crystals can be obtained from the estimated stiffness constants using the expression,
B=
1 (C11 + 2C12). 3
(15)
The bulk modulus obtained from our simulation is shown in Fig. 4. We have included in the figure the computational result from Herrero and Ramírez [8] and the available experimental results obtained by Migliori et al. using resonant-ultrasound spectroscopy [30] and those by Zouboulis et al. extracted from the Brillouin scattering [33]. Also given is
B (GPa)
440
420
400
380 0
400
800
1200
1600
2000
Temperature (K) Fig. 4. The isothermal bulk modulus of diamond as a function of temperature at hydrostatic pressure 1 atm. The red open (green filled) circles are the PIMC (classical MC) results. The blue lozenges are the PIMC results from Herreo and Ramírez [8]. The black stars and the black solid curve indicate the experimental results from Refs. [30,33], respectively. The red solid curve is the quadratic fitting of our PIMC results in the temperature range 0 to 1600 K with an upward shift of 26.7 GPa. The red and green dotted lines are guides for the eyes. 4
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
ensembles to study the elastic stiffness of diamond at a constant pressure of 1 atm. The Tersoff potential is used to describe the interatomic interactions and the Hooke’s law for the stress-strain relation to determine the elastic constants. Our PIMC simulations with the Tersoff potential demonstrate that the theory can pick up the essential physics of the elastic properties, showing that C11, C44 , as well as the bulk modulus B, behave as quadratic functions of temperature up to 1600 K and C12 is temperature independent in consistence with the experimental observations. Quantitatively, the differences between theoretical estimates of the elastic constants and the experimental results vary from 3% to 16%. The obtained PIMC and classical MC results at low temperatures are very different. For temperatures below 1000 K, the quantum effects are significant on C11 and C44 , but not C12 . We find that the anisotropic factor is practically independent of temperature with the average value 1.35 ± 0.07 indicating that diamond is an anisotropic crystal even at ambient pressure. The PIMC simulation also predicts that the Cauchy pressure is negative and varies little at low temperatures up to 500 K. The Poisson ratio of the diamond is found to be about 0.09 at low temperatures and 0.102 at T = 2000 K.
0.105 0.100
ν
0.095 0.090 0.085 0.080 0
400
800
1200
1600
2000
Temperature (K) Fig. 5. The Poisson ratio as a function of the temperature. The red open (green filled) circles with error bars are the PIMC (classical MC) results. The red and green dotted lines are guides for the eyes.
CRediT authorship contribution statement
the result of the classical MC calculation. The difference between the PIMC and classical MC results indicates the quantum effects are responsible for softening the diamond bulk modulus in about 17 GPa at zero temperature. At high temperatures the quantum and classical results of the bulk modulus approach to each other. The isothermal bulk modulus from the PIMC simulations of Herrero and Ramírez were obtained from the volume mean-square fluctuation of the simulation cell using the Tersoff potential [8,14]. The differences between our simulation results and those obtained by Herrero and Ramírez are very small within the error bars. The PIMC bulk modulus results for T < 1600 K can be described by the function B (T ) = 419.3 11.5 × 10 6T 2 and its upward shift by 26.7 GPa given by the red solid curve agrees very well with the experimental data. The bulk modulus from the PIMC simulation presents a very similar temperature-dependent behavior as indicated by the experiments, but quantitatively it is about 6% smaller at low temperature. Such a quantitative difference is consistent with the above discussion on the results of the elastic constants. We end up the discussion by analyzing the Poisson’s ratio and its temperature dependence. This ratio is a measure of the Poisson effect due to the transverse deformation to any longitudinal stretching of the material [38]. For an isotropic material the Poisson ratio is defined as [39]
=
s12 C12 = , s11 C12 + C11
Braulio Gabriel Alencar Brito: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Guoqiang Hai: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Ladir Cândido: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research was supported by Brazilian agencies CNPq, FAPESP and FAPEG/PRONEX under Grant No. 201710267000503. The authors acknowledge computational resources from the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, URL: http://sdumont.lncc.br, and the CENAPAD-SP (Centro Nacional de Processamento de Alto Desempenho em São Paulo).
(16)
where s11 and s12 are the compliance of the system. The range of Poisson’s ratio for isotropic solids is from 1 to 0.5. As discussed above, the diamond is an anisotropic crystal. Thus, its Poisson’s ratio is dependent on the crystal direction. Fig. 5 shows the temperature dependence of the Poisson ratio at 1 atm calculated within the isotropic approximation. For T 600 K, the PIMC value of the obtained Poisson ratio is around 0.092 and almost independent of temperature, while the classical MC result decreases linearly to the value 0.084 at zero temperature. In the whole temperature range under study, the Poisson ratio obtained from the PIMC simulation has a small variation with minimum about 0.090 at low temperature and maximum 0.102 at T = 2000 K. According to a recent review on the subject, the Poissons ratio for a single crystalline diamond is dependent on the crystal direction. It varies accordingly between 0.00786 and 0.115, with an average value of 0.0691 [40–42].
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.commatsci.2019.109387. References [1] S. Koizumi, C. Nebel, M. Nesladek, Physics and Applications of CVD Diamond, Wiley-VCH Verlag GmbH and Co. KGaA, 2008. [2] Z. Jenei, E. O’Bannon, S. Weir, H. Cynn, M. Lipp, W. Evans, Single crystal toroidal diamond anvils for highpressure experiments beyond 5 megabar, Nat. Commun. 1 (2018), https://doi.org/10.1038/s41467-018-06071-x. [3] F.P. Bundy, The p, t phase and reaction diagram for elemental carbon, 1979, J. Geophys. Res. 85 (1980) 6930, https://doi.org/10.1029/JB085iB12p06930. [4] L.R. Radovic, Chemistry and Physics of Carbon vol. 29, Marcel Dekker, New York, 2004. [5] D.M. Gruen, O.A. Shenderova, A.Y. Vul, Synthesis, properties and applications of ultrananocrystalline diamond, NATO Science Series 192. [6] F. Mashali, E.M. Languri, J. Davidson, D. Kerns, W. Johnson, K. Nawaz,
4. Conclusions In this paper, we have carried out both the quantum path-integral Monte Carlo and classical Monte Carlo simulations in the NPT and NVT 5
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
G. Cunningham, Thermo-physical properties of diamond nanofluids: a review, Int. J. Heat Mass Transfer 129 (2019) 1123, https://doi.org/10.1016/j. ijheatmasstransfer.2018.10.033. R. Ramírez, C.P. Herrero, E.R. Hernández, Path-integral molecular dynamics simulation of diamond, Phys. Rev. B 73 (2006) 245202, , https://doi.org/10.1103/ PhysRevB.73.245202. C.P. Herrero, R. Ramírez, Structural and thermodynamic properties of diamond: a path-integral monte carlo study, Phys. Rev. B 63 (2000) 024103, , https://doi.org/ 10.1103/PhysRevB.63.024103. G. Kresse, J. Furthmüller, J. Hafner, Ab initio force constant approach to phonon dispersion relations of diamond and graphite, Europhys. Lett. 32 (2015) 729, https://doi.org/10.1209/0295-5075/32/9/005. R.Z. Khaliullin, H. Eshet, T.D. Kuhne, J. Behler, M. Parrinello, Nucleation mechanism for the direct graphite-to-diamond phase transition, Nat. Mater. 10 (2011) 693, https://doi.org/10.1038/nmat307810.1038/nmat3078. G. Thiering, A. Gali, Characterization of oxygen defects in diamond by means of density functional theory calculations, Phys. Rev. B 94 (2016) 125202, , https:// doi.org/10.1103/PhysRevB.94.125202. B.G.A. Brito, A. Antonelli, Efficient method to include nuclear quantum effects in the determination of phase boundaries, J. Chem. Phys. 137 (3) (2012) 034114. B.G.A. Brito, L. Cândido, G.-Q. Hai, F.M. Peeters, Quantum effects in a free-standing graphene lattice: path-integral against classical monte carlo simulations, Phys. Rev. B 92 (2015) 195416, , https://doi.org/10.1103/PhysRevB.92.195416. B.G.A. Brito, L. Cândido, J.N.R. Teixeira, G.-Q. Hai, Thermodynamic properties of solid molecular hydrogen by path integral monte carlo simulations, Chem. Phys. Lett. 691 (2018) 330, https://doi.org/10.1016/j.cplett.2017.11.043. D. Marx, M. Parrinello, Structural quantum effects and three-centre two-electron bonding in ch+5, Nature 375 (1995) 216, https://doi.org/10.1038/375216a0. D. Marx, M. Parrinello, The effect of quantum and thermal fluctuations on the structure of the floppy molecule c2h3+, Science 271 (1996) 179, https://doi.org/ 10.1126/science.271.5246.179. C.P. Herrero, Path-integral monte carlo study of amorphous silicon, J. Non-Cryst. Solids 271 (1) (2000) 18–28, https://doi.org/10.1016/S0022-3093(00)00086-7 URL:http://www.sciencedirect.com/science/article/pii/S0022309300000867. D.M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67 (1995) 279–355, https://doi.org/10.1103/RevModPhys. 67.279. J.M. McMahon, M.A. Morales, C. Pierleoni, D.M. Ceperley, The properties of hydrogen and helium under extreme conditions, Rev. Mod. Phys. 84 (2012) 1607–1653, https://doi.org/10.1103/RevModPhys. 84.1607. C. Pierleoni, D.M. Ceperley, B. Bernu, W.R. Magro, Equation of state of the hydrogen plasma by path integral monte carlo simulation, Phys. Rev. Lett. 73 (1994) 2145–2149, https://doi.org/10.1103/PhysRevLett. 73.2145. J. Tersoff, Empirical interatomic potential for carbon, with applications to amorphous carbon, Phys. Rev. Lett. 61 (1988) 2879, https://doi.org/10.1103/ PhysRevLett.61.2879. E.L. Pollock, D.M. Ceperley, Simulation of quantum many-body systems by pathintegral methods, Phys. Rev. B 30 (1984) 2555–2568, https://doi.org/10.1103/ PhysRevB.30.2555. M. Parrinello, R.A., Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182, https://doi.org/10.1063/1. 328693. R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill
Book Company, 1965. [25] M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2010. [26] L.A.P. Ardila, S.A. Vitiello, M. Koning, Elastic constants of hcp 4 he: path-integral monte carlo results versus experiment, Phys. Rev. B 84 (2011) 094119, , https:// doi.org/10.1103/PhysRevB.84.094119. [27] M. Sprik, R.W. Impey, M.L. Klein, Second-order elastic constants for the lennardjones solid, Phys. Rev. B 29 (1984) 4368, https://doi.org/10.1103/PhysRevB.29. 4368. [28] M. Parrinello, A. Rahman, Strain fluctuations and elastic constants, J. Chem. Phys. 76 (1982) 2662, https://doi.org/10.1063/1.443248. [29] D.R. Squire, A.C. Holt, W.G. Hoover, Isothermal elastic constants for argon theory and monte carlo calculations, Physica 42 (1969) 388, https://doi.org/10.1016/ 0031-8914(69)90031-7. [30] A. Migliori, H. Ledbetter, R.G. Leisure, C. Pantea, J.B. Betts, Diamond’s elastic stiffnesses from 322 k to 10 k, J. Appl. Phys. 104 (2008) 053512, , https://doi.org/ 10.1063/1.2975190. [31] H.J. McSkimin, P. Andreatch, Elastic moduli of diamond as a function of pressure and temperature, J. Appl. Phys. 43 (1972) 2944, https://doi.org/10.1063/1. 1661636. [32] A. Nagakubo, M. Arita, H. Ogi, H. Sumiya, N. Nakamura, M. Hirao, Elastic constant c 11 of 12 c diamond between 10 and 613 k, App. Phys. Lett. 108 (2016) 221902, , https://doi.org/10.1063/1.4952613. [33] E.S. Zouboulis, M. Grimsditch, A.K. Ramdas, S. Rodriguez, Temperature dependence of the elastic moduli of diamond: a brillouin-scattering study, Phys. Rev. B 57 (1998) 2889, https://doi.org/10.1103/PhysRevB.57.2889. [34] C. Zener, Elasticity and Anelasticity of Metals, University of Chicago, Chicago, 1948. [35] M. Rajagopalan, Electronic topological transition in nb 3 al under compression: an ab initio study, Physica B 413 (2013) 1, https://doi.org/10.1016/j.physb.2012.12. 029. [36] N. Miao, B. Sa, J. Zhou, Z. Sun, Theoretical investigation on the transition-metal borides with ta 3 b 4 -type structure: a class of hard and refractory materials, Comput. Mater. Sci. 50 (2011) 1559, https://doi.org/10.1016/j.commatsci.2010. 12.015. [37] H. Niu, X.-Q. Chen, P. Liu, W. Xing, X. Cheng, D. Li, Y. Li, Extra-electron induced covalent strengthening and generalization of intrinsic ductile-to-brittle criterion, Scientific Rep. 2 (2012) 718, https://doi.org/10.1038/srep00718. [38] R.S. Lakes, Negative-poisson’s-ratio materials: auxetic solids, Annu. Rev. Mater. Res. 47 (2017) 63, https://doi.org/10.1146/annurev-matsci-070616-124118. [39] J.F. Nye, Physical Properties of Crystals: Their Representation by Tensors and Matrices, Science Publications, Oxford, 2006. [40] P. Hess, The mechanical properties of various chemical vapor deposition diamond structures compared to the ideal single crystal, J. Appl. Phys. 111 (2012) 051101, , https://doi.org/10.1063/1.3683544. [41] J. Philip, P. Hess, T. Feygelson, J.E. Butler, S. Chattopadhyay, K.H. Chen, L.C. Chen, Elastic, mechanical, and thermal properties of nanocrystalline diamond films, J. Appl. Phys. 93 (2003) 2164, https://doi.org/10.1063/1.1537465. [42] C.A. Klein, G.F. Cardinale, Young’s modulus and poisson’s ratio of cvd diamond, Diamond Rel. Mater. 2 (1993) 918, https://doi.org/10.1016/0925-9635(93) 90250-6.
6
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Quantum effects on elastic constants of diamond by path-integral Monte Carlo simulations B.G.A. Britoa, G.-Q. Haib, L. Cândidoc,
⁎
a
Departamento de Física, Instituto de Ciências Exatas e Naturais e Educação (ICENE), Universidade Federal do Triângulo Mineiro – UFTM, 38064-200 Uberaba, MG, Brazil b Instituto de Física de São Carlos, Universidade de São Paulo, 13560-970 São Carlos, SP, Brazil c Instituto de Física, Universidade Federal de Goiás, 74690-900 Goiânia, GO, Brazil
ARTICLE INFO
ABSTRACT
Keywords: Diamond Quantum effects on elastic constants Path integral Monte Carlo
Using the path-integral Monte Carlo method, we investigate the quantum effects on the elastic constants C11, C12 , and C44 of the diamond crystal in a wide temperature range at ambient pressure. The Tersoff potential is used to describe the interatomic interactions and the elastic constants are determined by a direct method derived from the stress-strain curve. We find that the elastic constants C11 and C44 behave as a quadratic function of temperature at low temperatures and C12 is practically a constant being independent of temperature. The quantum effects are significant at low temperatures up to about 1000 K. The quantitative differences between the PIMC calculations and the experimental measurements of the elastic constants C11, C12 , and C44 at low temperatures are about 3%, 17% , and 9%, respectively. Using the obtained elastic constants, we also estimate the bulk modulus B, anisotropic factor A, Cauchy pressure PCauchy , and Poisson ratio of the diamond crystal.
1. Introduction Diamond is the hardest substance known to man. It has excellent transparency over an extensive wavelength range, high thermal conductivity at room temperature, and very high resistance to chemical corrosion [1,2]. For example, in the study of high pressure materials, the hardness property has been an important characteristic used for diamond anvil cell. Both the diamond and graphite are composed of the single chemical element: carbon. The difference between these two solids is the arrangement of the carbon atoms in the lattices and, at ambient conditions, graphite is thermodynamically the most stable allotrope of carbon by a small enthalpy difference of 2.9 kJ mol−1 [3]. However, there is an almost prohibitive activation barrier between them making their interconversion impossible at ambient conditions [4]. The mechanical and electronic properties of diamond have been widely investigated [5,6]. There are studies showing that the quantum nature of atomic motion affects the thermal and mechanical properties of the material at room temperature [7,8]. Quantities such as vibrational energy, kinetic energy, heat capacity, and bulk modulus deviate significantly from the classical predictions even at ambient conditions. It is known that the methods such as density functional theory (DFT) and molecular dynamics (MD) are able to accurately predict several ⁎
properties of covalent solid, in particular, the lattice parameter, phonon dispersion relation [9], phase transition [10], and dynamics of point and extended defect formation [11]. However, within these theoretical frameworks, there are difficulties in dealing with the zero-point vibrational motion of the atoms and handling quantum effects at finite temperatures. Path-integral Monte Carlo (PIMC) is a powerful method that treats the atomic nuclei as quantum particles naturally. PIMC is versatile and can use different approximations such as the empirical potential [12–14], or those based on standard electronic structure methods [15,16]. It has been considered as one of the most important methods to simulate a quantum system at finite temperature with essential contributions to our understanding of the thermodynamic properties of several systems, mainly in the condensed matter [8,13,14,17–20]. In this paper, we apply the PIMC method to investigate the quantum effects on the elastic-stiffness coefficient of the diamond. Our simulation shows that the quantum effects impact the elastic constants and related quantities of the diamond considerably even at room temperature. To evidence the quantum effects in this system, we will compare the PIMC calculation results with those obtained from the classical Monte Carlo (MC) simulations and also with available experimental results. The paper is organized as follows. In Section 2 we present a short description of the PIMC formalism deriving an expression for the
Corresponding author. E-mail address: [email protected] (L. Cândido).
https://doi.org/10.1016/j.commatsci.2019.109387 Received 16 August 2019; Received in revised form 25 October 2019; Accepted 26 October 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Please cite this article as: B.G.A. Brito, G.-Q. Hai and L. Cândido, Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109387
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
stress-tensor observable. In Section 3 we present the results for elasticstiffness coefficients C11, C12 , and C44 as a function of temperature. Using the estimated elastic constants, we also obtain other quantities such as the bulk modulus, anisotropic factor, Cauchy pressure, and Poisson ratio. Finally, we summarize our work in Section 4.
The above equation, i.e. Eq. (7), is exact. But, in the calculations we use the action within the primitive approximation, given by
2. Numerical approach and computational details
Spri=
ij
We perform both the PIMC and classical MC simulations for a diamond crystal composed of N carbon atoms in the isothermal-isobaric (NPT) and isothermal-isochoric (NVT) ensembles in a wide temperature range from 50 to 2000 K and at pressure of 1 atm. In the PIMC simulations, the carbon atom is considered as a quantum particle of mass 12.011 a.u. Each atom is modeled by a set of M beads, forming a cyclic chain or a polymer. For the interatomic interactions, we employ the widely used Tersoff potential [21]. The specific parameters for the Tersoff potential used in this work can be found in the paper by Herrero and Ramírez in Ref. [8]. The PIMC method [22] is based on the evaluation of the quantum density matrix
Sk
d R1…d RM 1exp k=1
1
dRdR (R, R , ) R|O|R ,
dR 0 dRM
1
× exp( S (R 0 RM 1, R0 )).
ij
F (N , V , T ) =
lnZ (N , V , T ).
1 det(h )
3
hjk k=1
F hik
. N,T
.
N
M
=1
=1
(7)
|r
,
r
2 , + 1|
U (r , ),
(8)
NM
ij
mM 2
N
M
3
=1
=1
l=1
hil (s
hik (s
,
s , )k
M
N
M
=1
=1
,
V (h s , )i
s , )l, 3
hjk (s , )k
.
k=1
(9) Thus, the expectation value of the stress tensor components depends on the beads position difference, and of the first derivative of the potential related to the position of the beads [26]. There are two methods for estimating the elastic constants using computational simulations: (i) The direct method [27,26] which derives elastic constant from the stress-strain curve, and (ii) the fluctuation method [28,29] which use ensemble averages of the fluctuation of the strain or stress. Here, we use the direct method to obtain the elastic stiffness constants through the relation between stress and strain by the Hooke’s law
(2)
ij
=
Cijkl
kl,
(10)
k, l
(3)
where Cijkl represents the elastic stiffness tensor and kl is the kl component of the strain tensor. The PIMC simulations are performed in a three-dimensional box using a diamond unit cell with periodic boundary conditions in the three directions. Most calculations are performed within a simulation box containing 216 atoms. A “time step” = 1.333 × 10 4 K−1 (corresponding to M = 150 at T = 50 K) is used for the discretized imaginary time path integrals. The extrapolation to the = 0 limit shows that this approximation induces an error of about 0.1% in the obtained total energy. Checks for larger systems containing 512 atoms show that the obtained results do not change within the statistical error when comparing to a system with 216 atoms. For the averages, we use 200.000 PIMC time steps.
(4)
3. Results and discussions We start the simulation using the NPT ensemble to determine the system density as a function of temperature at a constant pressure of 1 atm. Then, the obtained densities are used as input in the NVT simulation to calculate the elastic constants. The NVT simulation cell uses its vectors parallel to the edges of the unit cell of the diamond lattice. If the simulation box is scaled or deformed by some strain, we have the new vectors of the cell given by the relation
(5)
matrix h in the following way [25],
=
=1
det(h )
(1)
We can write an estimator for the components of the stress tensor ( ) as a function of the derivative of the free energy F by the elements of the
ij
=1
k=1
where det(h ) gives the volume of the simulation cell. The connection of the path-integral formalism with the thermodynamic observables in NVT ensemble is made through the Helmholtz free energy (F) obtained by the partition function [24],
1
M
Mm 2 2
3
tween 0 and 1, and h is the matrix used in the Parrinello and Rahman approximation to describe the homogeneous distortions of the periodic simulation cell [23]. In the dimensionless coordinates, the partition function given by Eq. (3) can be re-written as
det(h ) NM dR 0 dRM 1× exp( S (R 0 RM 1, R 0)),
N
1
=
×
The bead position ri, k can be written as ri, k = h si, k , where si, k is the dimensionless coordinate whose components have values ranging be-
Z (N , V , T )=
k=1
where the first two terms are kinetic energy and the last potential energy. Eqs. (7) and (8) as shown above were first obtained by Ardila et al. [26]. By substituting Spri in Eq. (8) into Eq. (7) we get the thermodynamic estimator for the stress tensor components, given by
and the partition function is written as
Z (N , V , T ) =
M
S hik
ij
2 3NM ln 2 2Mm
+
by sampling the set of paths {R 0, R1, …RM 1, RM} with R k being the set of beads {r1, k , …, rN , k} , and ri, k a bead that is the position of the ith particle in the k th time slice. The “time step” is given by = / M with = 1/ kB T (kB is the Boltzmann constant and T the absolute temperature). The action of a link is defined as Sk S (R k 1, R k; ) ln[ (Rk 1, Rk ; )] is handled in the primitive approximation. To evaluate the 3NM -dimensional integrals in Eq. (1), we employ the Metropolis bisection sampling technique [18]. The expectation value of an operator O is given by
O =Z
NM
det(h )
M
(R 0 , RM ; ) =
3
1
=
(6)
Ai
It can be re-written in terms of the expectation value of the derivative of the action S by substituting the free energy given by Eq. (4) into Eq. (6), and taking the appropriate derivatives leading to
new
= (I ±
ij )
Ai ,
(11)
where Ai are the edge vectors and I is the unit matrix. Two strain tensors 2
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
10
(a)
1100
C11 (GPa)
σ xx (GPa)
5 0
(a)
1050
1000
−5 950
−10
0
−0.01
−0.005
0
0.005
400
0.01
800 1200 1600 Temperature (K)
εxx
(b)
160
6
C12 (GPa)
(b)
σ xy (GPa)
4 2
2000
140 120
0
100
−2
80 0
−4
400
800
1200
1600
2000
Temperature (K)
−6 −0.01
−0.005
0
0.005
(c)
0.01
640
εxy
xx
e 0 0 = 0 0 0 0 0 0
C44 (GPa)
Fig. 1. Responses of (a) the tensile stress xx under the imposed strain xx and (b) the shear stress xy under the imposed strain xy at temperature 300 K and pressure 1 atm.
(12)
0
0 e /2 0 = e /2 0 0 0 0 0
yy xy
= C11 xx , = zz = C12 = C44 xy.
400
800
1200
1600
2000
Temperature (K) Fig. 2. Temperature dependence of the elastic constants (a) C11 , (b) C12 , and (c) C44 of diamond at pressure 1 atm. The red open and green filled circles are results from the PIMC and classical MC simulations, respectively, with the red and green dotted lines being guides for the eye. Different experimental measurements are given by the black stars [30], black triangles [31], black squares [32], and the black solid curves with the error range indicated by two black dotted lines [33]. The red solid curves are the quadratic fittings of the PIMC results of C11 and C44 at low temperatures in the range 0 to 1600 K and an average value of C12 shifted upward by 34.8 GPa in (a) and 21.5 GPa in (b) and shifted downward by 52.1 GPa in (c).
(13)
are used to determine the elastic constants by the relations: xx
560
520
and xy
600
xx ,
(14)
The values of ij are obtained from the ensemble average in the NVT simulation with e = 0, ± 0,005, and ± 0.01 in Eqs. (12) and (13). The elastic constants are calculated from linear strain-stress relations given in Eq. (14). Fig. 1(a) and (b) give the obtained stresses xx and xy , respectively, as a function of applied strains at temperature 300 K and pressure 1 atm. In both cases, the stress shows a linear dependence on the strain for small distortions. Fig. 1(a) shows the tensile stress response xx as a function of the imposed tensile strain xx . For xx = 0 , we have only the hydrostatic compression, which in this case is small and falls within the error bar. The slope of the line gives the elastic constant C11 for longitudinal compression. Fig. 1(b) shows a similar response curve for shear stress xy in response of the shear strain xy . As expected for this kind of strain, the linear fit passes through the origin indicating no shear stress
at 1 atm, in other words, zero shear deformation. We verified that, in this range of strain, the stress-strain dependence is almost perfect linear for all the studied temperatures. When a tensile stress is applied along the x-axis of the diamond crystal, using Eq. (14) we can estimate the elastic constants C11 and C12 which are associated with the stress responses in the directions parallel with (longitudinal compression) and perpendicular to (transverse expansion) the x-axis, respectively. While C44 is shear modulus given by xy / xy . Fig. 2 shows C11, C12 , and C44 as a function of temperature at constant pressure 1 atm. The results obtained from the PIMC and classical MC calculations are represented by red open and green filled 3
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
circles, respectively. The experimental results from Migliori et al. [30], McSckimin and Adreatch [31], and Nagakubo et al. [32] are shown by black stars, triangles, and squares, respectively. The black solid curves give the experimental results from Zouboulis et al. [33] with two black dotted curves indicating the estimate experimental error range. Comparing C11 in Fig. 2(a) with C12 in 2(b), one can notice that the quantum effects on C11 are much more significant than on C12 . The classical and quantum results of C12 are close to each other with overlaps of the error bars in the whole temperature range. The PIMC results of C12 are practically independent of temperature with an average value of 105.3 GPa with error bars varying from 2 to 6 GPa. An upward shift of this average value by 21.5 GPa is shown by the red solid line in Fig. 2(b) which is in agreement with the experimental results of Zouboulis et al. [33]. However, the temperature dependence of C11 is significant. The classical MC predicts a linear decrease of C11 with increasing temperature while the PIMC results show a quadratic function dependence at low temperatures. At high temperatures, the PIMC and classical MC results of C11 approach to each other as expected. A numerical fitting of the PIMC results in the temperature range 0 < T < 1600 K yields (0) (2) 2 (0) (2) C11 c11 c11 T with c11 = 1048.6 GPa and c11 = 33.5 × 10 6 GPa/K2. Its upward shift by 34.8 GPa is given by the red solid curve in Fig. 2(a) showing a good agreement with the available experimental results. Fig. 2(c) shows the elastic constant C44 as a function of temperature. C44 is also known as shear modulus µ . The shear modulus describes the tendency of the object suffering deformation of its shape and keeping constant volume. Being similar to C11, the classical MC simulation shows that C44 decreases linearly with increasing temperature. The quantum effects on C44 are more significant at lower temperature. For T < 1600 K, the PIMC values of C44 are well described by the quadratic (0) (2) 2 (0) C44 c44 c44 T , c44 = 630.3 GPa function with and (2) c44 = 18.1 × 10 6 GPa/K2. Its downward shift by 52.1 GPa is in very good agreement with the experimental results as shown by the red solid curve in Fig. 2(c). In general, our PIMC simulations using the Tersoff’s potential demonstrate that the theory can describe correctly the temperature dependence of the elastic constants C11, C12 , and C44 . The domination of the quantum effects at low temperatures leads to the elastic constants C11 and C44 decrease quadratically with increasing temperature. Quantitatively, considering the experimental results as references, the PIMC calculations underestimate C11 and C12 at low temperatures by about 3% and 17%, respectively, and overestimate C44 by 9% . The anisotropy factor for a cubic crystal can be defined [34] by the 2C ratio A = C 44C . Because C44 = (C11 C12 )/2 for a cubic isotropic 11 12 crystal, the value of A is expected to be equal to 1 when the crystal is isotropic, while any deviation from the unity indicates a certain degree of elastic anisotropy. There have been studies showing that elastic anisotropy of crystal is important for material engineering due to its correlation with the producing of microcracks in materials [35]. For the temperature range under investigation, we find the value of A from the PIMC calculation is practically constant. The average value is 1.35 ± 0.07 showing that diamond is an anisotropic crystal even at ambient pressure. However, the value of A also indicates a negligible influence of the temperature on the elastic anisotropy of the diamond crystal. From the above calculations, we can estimate the Cauchy pressure for the system, giving by PCauchy = C12 C44 . A negative (positive) Cauchy pressure reveals brittleness (ductility) of the crystal [36,37]. The more negative PCauchy , the more angular character of atomic bonding in the crystal. Our PIMC and MC simulations using the Tersoff potential predict correctly a very negative Cauchy pressure of the diamond for all temperatures as shown in Fig. 3. This is a consequence of the correct directional bonding description for the carbon atoms in diamond lattice given by the Tersoff potential. In fact, diamond is the hardest known of the highest strength material and at the same time, the most brittle one of the most negative Cauchy pressure [37]. The
−460
P Cauchy (GPa)
−480
−500
−520
−540
−560 0
400
800
1200
1600
2000
Temperature (K) Fig. 3. The Cauchy pressure of diamond as a function of temperature. The red open and green filled circles are the results from the PIMC and classical MC simulations, respectively.
same negative behavior cannot be observed using pure pair potentials because they have no preferred direction of bonding. We also notice that the PCauchy value obtained from the classical MC approaches to the PIMC result at high temperatures. However, the classical approximation yields more negative Cauchy pressure for temperatures below 1000 K, indicating a more directional bonding. This is a consequence of the reduction of the atomic vibrational energy in the classical approximation where only thermal energy is considered. The PIMC simulation takes into account the quantum effects due to the zero-point energy of the atomic motion and predicts correctly the weak temperature dependence of the Cauchy pressure at low temperatures. The bulk modulus (B) for cubic crystals can be obtained from the estimated stiffness constants using the expression,
B=
1 (C11 + 2C12). 3
(15)
The bulk modulus obtained from our simulation is shown in Fig. 4. We have included in the figure the computational result from Herrero and Ramírez [8] and the available experimental results obtained by Migliori et al. using resonant-ultrasound spectroscopy [30] and those by Zouboulis et al. extracted from the Brillouin scattering [33]. Also given is
B (GPa)
440
420
400
380 0
400
800
1200
1600
2000
Temperature (K) Fig. 4. The isothermal bulk modulus of diamond as a function of temperature at hydrostatic pressure 1 atm. The red open (green filled) circles are the PIMC (classical MC) results. The blue lozenges are the PIMC results from Herreo and Ramírez [8]. The black stars and the black solid curve indicate the experimental results from Refs. [30,33], respectively. The red solid curve is the quadratic fitting of our PIMC results in the temperature range 0 to 1600 K with an upward shift of 26.7 GPa. The red and green dotted lines are guides for the eyes. 4
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
ensembles to study the elastic stiffness of diamond at a constant pressure of 1 atm. The Tersoff potential is used to describe the interatomic interactions and the Hooke’s law for the stress-strain relation to determine the elastic constants. Our PIMC simulations with the Tersoff potential demonstrate that the theory can pick up the essential physics of the elastic properties, showing that C11, C44 , as well as the bulk modulus B, behave as quadratic functions of temperature up to 1600 K and C12 is temperature independent in consistence with the experimental observations. Quantitatively, the differences between theoretical estimates of the elastic constants and the experimental results vary from 3% to 16%. The obtained PIMC and classical MC results at low temperatures are very different. For temperatures below 1000 K, the quantum effects are significant on C11 and C44 , but not C12 . We find that the anisotropic factor is practically independent of temperature with the average value 1.35 ± 0.07 indicating that diamond is an anisotropic crystal even at ambient pressure. The PIMC simulation also predicts that the Cauchy pressure is negative and varies little at low temperatures up to 500 K. The Poisson ratio of the diamond is found to be about 0.09 at low temperatures and 0.102 at T = 2000 K.
0.105 0.100
ν
0.095 0.090 0.085 0.080 0
400
800
1200
1600
2000
Temperature (K) Fig. 5. The Poisson ratio as a function of the temperature. The red open (green filled) circles with error bars are the PIMC (classical MC) results. The red and green dotted lines are guides for the eyes.
CRediT authorship contribution statement
the result of the classical MC calculation. The difference between the PIMC and classical MC results indicates the quantum effects are responsible for softening the diamond bulk modulus in about 17 GPa at zero temperature. At high temperatures the quantum and classical results of the bulk modulus approach to each other. The isothermal bulk modulus from the PIMC simulations of Herrero and Ramírez were obtained from the volume mean-square fluctuation of the simulation cell using the Tersoff potential [8,14]. The differences between our simulation results and those obtained by Herrero and Ramírez are very small within the error bars. The PIMC bulk modulus results for T < 1600 K can be described by the function B (T ) = 419.3 11.5 × 10 6T 2 and its upward shift by 26.7 GPa given by the red solid curve agrees very well with the experimental data. The bulk modulus from the PIMC simulation presents a very similar temperature-dependent behavior as indicated by the experiments, but quantitatively it is about 6% smaller at low temperature. Such a quantitative difference is consistent with the above discussion on the results of the elastic constants. We end up the discussion by analyzing the Poisson’s ratio and its temperature dependence. This ratio is a measure of the Poisson effect due to the transverse deformation to any longitudinal stretching of the material [38]. For an isotropic material the Poisson ratio is defined as [39]
=
s12 C12 = , s11 C12 + C11
Braulio Gabriel Alencar Brito: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Guoqiang Hai: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Ladir Cândido: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research was supported by Brazilian agencies CNPq, FAPESP and FAPEG/PRONEX under Grant No. 201710267000503. The authors acknowledge computational resources from the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, URL: http://sdumont.lncc.br, and the CENAPAD-SP (Centro Nacional de Processamento de Alto Desempenho em São Paulo).
(16)
where s11 and s12 are the compliance of the system. The range of Poisson’s ratio for isotropic solids is from 1 to 0.5. As discussed above, the diamond is an anisotropic crystal. Thus, its Poisson’s ratio is dependent on the crystal direction. Fig. 5 shows the temperature dependence of the Poisson ratio at 1 atm calculated within the isotropic approximation. For T 600 K, the PIMC value of the obtained Poisson ratio is around 0.092 and almost independent of temperature, while the classical MC result decreases linearly to the value 0.084 at zero temperature. In the whole temperature range under study, the Poisson ratio obtained from the PIMC simulation has a small variation with minimum about 0.090 at low temperature and maximum 0.102 at T = 2000 K. According to a recent review on the subject, the Poissons ratio for a single crystalline diamond is dependent on the crystal direction. It varies accordingly between 0.00786 and 0.115, with an average value of 0.0691 [40–42].
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.commatsci.2019.109387. References [1] S. Koizumi, C. Nebel, M. Nesladek, Physics and Applications of CVD Diamond, Wiley-VCH Verlag GmbH and Co. KGaA, 2008. [2] Z. Jenei, E. O’Bannon, S. Weir, H. Cynn, M. Lipp, W. Evans, Single crystal toroidal diamond anvils for highpressure experiments beyond 5 megabar, Nat. Commun. 1 (2018), https://doi.org/10.1038/s41467-018-06071-x. [3] F.P. Bundy, The p, t phase and reaction diagram for elemental carbon, 1979, J. Geophys. Res. 85 (1980) 6930, https://doi.org/10.1029/JB085iB12p06930. [4] L.R. Radovic, Chemistry and Physics of Carbon vol. 29, Marcel Dekker, New York, 2004. [5] D.M. Gruen, O.A. Shenderova, A.Y. Vul, Synthesis, properties and applications of ultrananocrystalline diamond, NATO Science Series 192. [6] F. Mashali, E.M. Languri, J. Davidson, D. Kerns, W. Johnson, K. Nawaz,
4. Conclusions In this paper, we have carried out both the quantum path-integral Monte Carlo and classical Monte Carlo simulations in the NPT and NVT 5
Computational Materials Science xxx (xxxx) xxxx
B.G.A. Brito, et al.
[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
G. Cunningham, Thermo-physical properties of diamond nanofluids: a review, Int. J. Heat Mass Transfer 129 (2019) 1123, https://doi.org/10.1016/j. ijheatmasstransfer.2018.10.033. R. Ramírez, C.P. Herrero, E.R. Hernández, Path-integral molecular dynamics simulation of diamond, Phys. Rev. B 73 (2006) 245202, , https://doi.org/10.1103/ PhysRevB.73.245202. C.P. Herrero, R. Ramírez, Structural and thermodynamic properties of diamond: a path-integral monte carlo study, Phys. Rev. B 63 (2000) 024103, , https://doi.org/ 10.1103/PhysRevB.63.024103. G. Kresse, J. Furthmüller, J. Hafner, Ab initio force constant approach to phonon dispersion relations of diamond and graphite, Europhys. Lett. 32 (2015) 729, https://doi.org/10.1209/0295-5075/32/9/005. R.Z. Khaliullin, H. Eshet, T.D. Kuhne, J. Behler, M. Parrinello, Nucleation mechanism for the direct graphite-to-diamond phase transition, Nat. Mater. 10 (2011) 693, https://doi.org/10.1038/nmat307810.1038/nmat3078. G. Thiering, A. Gali, Characterization of oxygen defects in diamond by means of density functional theory calculations, Phys. Rev. B 94 (2016) 125202, , https:// doi.org/10.1103/PhysRevB.94.125202. B.G.A. Brito, A. Antonelli, Efficient method to include nuclear quantum effects in the determination of phase boundaries, J. Chem. Phys. 137 (3) (2012) 034114. B.G.A. Brito, L. Cândido, G.-Q. Hai, F.M. Peeters, Quantum effects in a free-standing graphene lattice: path-integral against classical monte carlo simulations, Phys. Rev. B 92 (2015) 195416, , https://doi.org/10.1103/PhysRevB.92.195416. B.G.A. Brito, L. Cândido, J.N.R. Teixeira, G.-Q. Hai, Thermodynamic properties of solid molecular hydrogen by path integral monte carlo simulations, Chem. Phys. Lett. 691 (2018) 330, https://doi.org/10.1016/j.cplett.2017.11.043. D. Marx, M. Parrinello, Structural quantum effects and three-centre two-electron bonding in ch+5, Nature 375 (1995) 216, https://doi.org/10.1038/375216a0. D. Marx, M. Parrinello, The effect of quantum and thermal fluctuations on the structure of the floppy molecule c2h3+, Science 271 (1996) 179, https://doi.org/ 10.1126/science.271.5246.179. C.P. Herrero, Path-integral monte carlo study of amorphous silicon, J. Non-Cryst. Solids 271 (1) (2000) 18–28, https://doi.org/10.1016/S0022-3093(00)00086-7 URL:http://www.sciencedirect.com/science/article/pii/S0022309300000867. D.M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67 (1995) 279–355, https://doi.org/10.1103/RevModPhys. 67.279. J.M. McMahon, M.A. Morales, C. Pierleoni, D.M. Ceperley, The properties of hydrogen and helium under extreme conditions, Rev. Mod. Phys. 84 (2012) 1607–1653, https://doi.org/10.1103/RevModPhys. 84.1607. C. Pierleoni, D.M. Ceperley, B. Bernu, W.R. Magro, Equation of state of the hydrogen plasma by path integral monte carlo simulation, Phys. Rev. Lett. 73 (1994) 2145–2149, https://doi.org/10.1103/PhysRevLett. 73.2145. J. Tersoff, Empirical interatomic potential for carbon, with applications to amorphous carbon, Phys. Rev. Lett. 61 (1988) 2879, https://doi.org/10.1103/ PhysRevLett.61.2879. E.L. Pollock, D.M. Ceperley, Simulation of quantum many-body systems by pathintegral methods, Phys. Rev. B 30 (1984) 2555–2568, https://doi.org/10.1103/ PhysRevB.30.2555. M. Parrinello, R.A., Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182, https://doi.org/10.1063/1. 328693. R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill
Book Company, 1965. [25] M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2010. [26] L.A.P. Ardila, S.A. Vitiello, M. Koning, Elastic constants of hcp 4 he: path-integral monte carlo results versus experiment, Phys. Rev. B 84 (2011) 094119, , https:// doi.org/10.1103/PhysRevB.84.094119. [27] M. Sprik, R.W. Impey, M.L. Klein, Second-order elastic constants for the lennardjones solid, Phys. Rev. B 29 (1984) 4368, https://doi.org/10.1103/PhysRevB.29. 4368. [28] M. Parrinello, A. Rahman, Strain fluctuations and elastic constants, J. Chem. Phys. 76 (1982) 2662, https://doi.org/10.1063/1.443248. [29] D.R. Squire, A.C. Holt, W.G. Hoover, Isothermal elastic constants for argon theory and monte carlo calculations, Physica 42 (1969) 388, https://doi.org/10.1016/ 0031-8914(69)90031-7. [30] A. Migliori, H. Ledbetter, R.G. Leisure, C. Pantea, J.B. Betts, Diamond’s elastic stiffnesses from 322 k to 10 k, J. Appl. Phys. 104 (2008) 053512, , https://doi.org/ 10.1063/1.2975190. [31] H.J. McSkimin, P. Andreatch, Elastic moduli of diamond as a function of pressure and temperature, J. Appl. Phys. 43 (1972) 2944, https://doi.org/10.1063/1. 1661636. [32] A. Nagakubo, M. Arita, H. Ogi, H. Sumiya, N. Nakamura, M. Hirao, Elastic constant c 11 of 12 c diamond between 10 and 613 k, App. Phys. Lett. 108 (2016) 221902, , https://doi.org/10.1063/1.4952613. [33] E.S. Zouboulis, M. Grimsditch, A.K. Ramdas, S. Rodriguez, Temperature dependence of the elastic moduli of diamond: a brillouin-scattering study, Phys. Rev. B 57 (1998) 2889, https://doi.org/10.1103/PhysRevB.57.2889. [34] C. Zener, Elasticity and Anelasticity of Metals, University of Chicago, Chicago, 1948. [35] M. Rajagopalan, Electronic topological transition in nb 3 al under compression: an ab initio study, Physica B 413 (2013) 1, https://doi.org/10.1016/j.physb.2012.12. 029. [36] N. Miao, B. Sa, J. Zhou, Z. Sun, Theoretical investigation on the transition-metal borides with ta 3 b 4 -type structure: a class of hard and refractory materials, Comput. Mater. Sci. 50 (2011) 1559, https://doi.org/10.1016/j.commatsci.2010. 12.015. [37] H. Niu, X.-Q. Chen, P. Liu, W. Xing, X. Cheng, D. Li, Y. Li, Extra-electron induced covalent strengthening and generalization of intrinsic ductile-to-brittle criterion, Scientific Rep. 2 (2012) 718, https://doi.org/10.1038/srep00718. [38] R.S. Lakes, Negative-poisson’s-ratio materials: auxetic solids, Annu. Rev. Mater. Res. 47 (2017) 63, https://doi.org/10.1146/annurev-matsci-070616-124118. [39] J.F. Nye, Physical Properties of Crystals: Their Representation by Tensors and Matrices, Science Publications, Oxford, 2006. [40] P. Hess, The mechanical properties of various chemical vapor deposition diamond structures compared to the ideal single crystal, J. Appl. Phys. 111 (2012) 051101, , https://doi.org/10.1063/1.3683544. [41] J. Philip, P. Hess, T. Feygelson, J.E. Butler, S. Chattopadhyay, K.H. Chen, L.C. Chen, Elastic, mechanical, and thermal properties of nanocrystalline diamond films, J. Appl. Phys. 93 (2003) 2164, https://doi.org/10.1063/1.1537465. [42] C.A. Klein, G.F. Cardinale, Young’s modulus and poisson’s ratio of cvd diamond, Diamond Rel. Mater. 2 (1993) 918, https://doi.org/10.1016/0925-9635(93) 90250-6.
6