Diffusion Quantum Monte Carlo study of the equation of state and point defects in aluminum
Abstract
The many-body diffusion quantum Monte Carlo (DMC) method with twist-averaged boundary conditions is used to calculate the ground-state equation of state and the energetics of point defects in fcc aluminum, using supercells up to 1331 atoms. The DMC equilibrium lattice constant differs from experiment by 0.008 Å or 0.2%, while the cohesive energy using DMC with backflow wave functions with improved nodal surfaces differs by 27 meV. DMC calculated defect formation and migration energies agree with available experimental data, except for the nearest-neighbor divacancy, which is found to be energetically unstable in agreement with previous density functional theory (DFT) calculations. DMC and DFT calculations of vacancy defects are in reasonably close agreement. Self-interstitial formation energies have larger differences between DMC and DFT, of up to 0.33eV, at the tetrahedral site. We also computed formation energies of helium interstitial defects where energies differed by up to 0.34eV, also at the tetrahedral site. The close agreement with available experiment demonstrates that DMC can be used as a predictive method to obtain benchmark energetics of defects in metals.
pacs:
61.72.J-, 64.30.Ef, 02.70.SsI Introduction
The mechanical properties of metals are dominated by the formation and migration energies of defects. Experimentally it is often difficult to measure desired defect properties directly, so it is important to have accurate theoretical approaches to calculate defect properties. The many-body diffusion quantum Monte Carlo (DMC) approach foulkes01 () is the most accurate method for systems with more than electrons but has not previously been applied to defects in metals. Until today, the most successful quantum mechanics-based defect calculations in metals use density functional theory (DFT). However, the approximate functionals used (i) lack sufficient specific and universal accuracy, (ii) can not be systematically improved, and (iii) there are now many approximate functionals to choose from, all giving different results. Thus, DMC results are an ideal candidate to wholly replace DFT when benchmark thermodynamic properties are required. Furthermore, the relative energies of reference phases (bulk metals and compounds) would be of great value in thermodynamic databases and in the subsequent prediction of phase diagrams (e.g. CALPHAD and Thermocalc).
In semiconductors, DMC calculations of the formation energies of point defects and surface energies have shown that important corrections to DFT arise when electronic correlations are fully taken into account. Deviations in formation energies of more than 1 eV were found in silicon leung99 () and diamond hood03 (). Additionally, activation energies of common chemical reactions obtained by DFT methods have been shown to differ substantially from benchmark diffusion Monte Carlo (DMC) values grossman97 () and quantum chemical results. For metallic systems, the size of the errors in DFT calculations of defects is largely unknown, as more accurate benchmark calculations do not currently exist. It is highly desirable to demonstrate the feasibility of DMC, with its increased predictive accuracy, as a replacement of DFT for challenging systems such as metals, particularly as computer power increases.
Aluminum is an ideal starting point for carrying out initial DMC calculations of defects since it is one of the simplest metals with a close-packed fcc structure that contains no 3d electrons. As a result it is considered a prototype material for testing the validity of theoretical calculations. Aluminum is well characterized experimentally so there is an abundance of data available. There have been previous quantum Monte Carlo calculations of the bulk properties of aluminum, however, these were done using the less accurate variational Monte Carlo gaudoin02a (); gaudoin02b () and the calculations had large statistical noise.
In this paper, we report well converged results for the bulk properties of fcc aluminum using DMC. We explore a larger range of volumes in order to compare the ground-state equation of state calculated with DFT. Our DMC calculations of the defect properties of aluminum include the simplest point defect, the vacancy for which numerous DFT carling00 (); mattsson06 (); delczeg09 () and experimental results erhart91 (); hehenkamp94 (); fluss78 (); seeger71 (); desorbo59 () are available. We compute the nearest-neighbor divacancy binding energy, a defect that DFT calculations carling00 (); uesugi03 () have found to be unstable. Since this instability is counter to both experimental studies hehenkamp94 (); fluss84 () and simple bond counting argumentsdamask59 () it is important to perform calculations of this defect with DMC, a method that unlike DFT, does not rely on approximations for exchange and correlation.
We also examine two other types of defects: First, self-interstitials can arise due to irradiation with energetic particles, through plastic deformation, or through their production in thermal equilibrium at high temperatures. We have obtained DMC results for the formation energies of the -dumbbell, the octahedral, and the tetrahedral self-interstitials. Second, experimental studies of irradiated aluminum show the presence of He bubbles rajainmaki88 (); birtcher94 (); hamaguchi04 (). The formation and growth of helium bubbles can alter a material’s mechanical properties through void swelling, embrittlement, and surface blistering katoh03 (); vassen91 (). Since irradiation, through He implantation or transmutation, gives rise to He atoms at substitutional or interstitial lattice sites, the energetics of these types of defects are important. Presented here are DMC calculated formation energies for He at the substitutional, octahedral and the tetrahedral interstitial sites, which are compared with previous DFT calculations yang08 (). We demonstrate that twist-averaged boundary conditions offer far superior statistics over DFT calculated single-particle finite-size corrections for DMC calculations of defects, and also remove a dependence on data from other methods. Overall the calculations reported here used several million processor hours, which is affordable on large computer clusters and supercomputers.
Ii Bulk Aluminum
The DMC approach foulkes01 () is a stochastic method for evolving a wave function using the imaginary-time Schrödinger equation. In principle, and in contrast to most electronic structure methods, the systematic errors can be measured and systematically reduced umrigarprl2007 (); bajdichprl2010 (). In practice, the most significant sources of error are (i) the fixed-node and fixed-phase approximations, a variational solution to the Fermion sign problem, (ii) adequately sampling the Brillouin zone, which in contrast with insulators is a significant problem in metals, and (iii) pseudopotential error and corresponding locality error (although these may be avoided through all-electron calculations esler10 ()).
For our DMC calculations we used the CASINO code needs10 (), with guiding wave functions formed by a product of Slater determinants for up and down spin electrons and a Jastrow correlation function. The single-particle orbitals in the determinants were obtained from DFT calculations using the generalized gradient approximation (GGA) for the exchange-correlation term since it should perform better than the local density approximation (LDA) where the electron distribution shows large spatial variations as it does at a vacancy in aluminum carling00 (). (For bulk aluminum GGA is more accurate than LDA, see Fig. 3.) For GGA we used the Perdew-Burke-Ernzerhof (PBE) form perdew96 () rather than the Perdew-Wang-91 (PW91) form perdew92 () since it gives better values mattsson06 () for the vacancy formation energy. The DFT calculations were performed using the plane-wave PWSCF code giannozzi09 () with Troullier-Martins non-local pseudopotentials. For the non-local pseudopotentials we used the locality approximation mitas91 () in the DMC calculations. The orbitals were evaluated in DMC via real-space cubic splines williamson01 (). For the DMC simulations we used 5280 walkers with a time step of 0.01 au, which our test calculations revealed gave time-step errors in Table 1 of less than 0.01 eV.
Calculations of solids using DMC require finite simulation supercells with periodic boundary conditions imposed on the Hamiltonian. This leads to finite-size effects that can be divided between single-particle and many-body contributions. In a metal the DMC kinetic energy contains large single-particle contributions due to the sharp Fermi surface, which impacts our calculations in two ways.
First, in a real metal the number of orbitals with energies below the Fermi level is usually not equal to the number of electrons required in the simulation cell. In DFT calculations one can use partial occupations of these orbitals to create a closed shell configuration guaranteeing that the charge density has the correct symmetry. In a DMC calculation with a guiding wave function containing a single determinant for the spin up and for the spin down electrons the use of partial occupations is not an option. Gaudoin, et. al. gaudoin02b () found differences in aluminum as large as 0.1 eV/atom in variational Monte Carlo total energies depending on the occupations of the orbitals at the Fermi level. As our calculations show in Fig. 3 for aluminum an error on the order of 0.1 eV/atom is too large to accurately discern the equilibrium lattice constant.
A second complication for DMC calculations of metals arises as the system size is varied, which can produce band crossings as energy levels pass through the Fermi level. This can create discontinuous changes in the nodal surface of the guiding wave function as the symmetry of the occupied orbitals change, leading to discontinuities in the DMC energy as shown in Fig. 1. For simulation cells containing 64 or 125 atoms we found that single-particle DFT corrections were ineffective in removing the large discontinuities we found in calculated DMC energies as the volume was varied. However, the use of well-converged twist-averaged boundary conditions lin01 (); teweldeberhan10 () was effective in producing smooth energy versus volume curves at these systems sizes as shown in Fig. 1. Our calculations with 64 and 125 atoms used -point centered grids with 13x13x13 and 10x10x10 twists, respectively.
With twist-averaged boundary conditions the remaining finite-size effects arise from many-body contributions drummond08 () that scale as
(1) |
As shown in Fig. 2 this equation agrees well with our DMC energies using twist-averaged boundary conditions for simulation cells containing between 27 and 1331 atoms.
We used Eq. (1) to obtain infinite-size extrapolated DMC total energies from calculations using 64 and 125 atom simulation cells with twist-averaged boundary conditions using a range of lattice constants shown in Fig. 3. The energy points were fit to a quartic and a Murnaghan equation of state murnaghan44 (). Both fits yielded an equilibrium lattice constant of 4.030(1) Å. This compares well with the experimental value, 4.022 Å, with zero-point energy and finite-temperature effects removed gaudoin02a (). In contrast DFT calculations with the local density approximation (LDA) and GGA (PBE) yield 3.960 Å and 4.046 Å, respectively.
For the DMC cohesive energy we initially obtained 3.341(1) eV using the fixed-nodes defined by the GGA orbitals, compared with the experimental value of 3.43 eV with zero-point energy and finite-temperature effects removed, 4.21 eV using LDA, and 3.52 eV using GGA (PBE). Although the fixed-node DMC results with GGA nodes are already the most accurate, to assess the possible nodal error we also performed DMC calculations with twist-averaged boundary conditions and extrapolation to an infinite-sized supercell using optimized backflow wave functions (with backflow transformations that contained electron-electron, electron-nuclei, and electron-electron-nuclei terms rios06 ()) at the optimum lattice constant and also for an isolated atom. Backflow wave functions can be substantially more accurate than single determinant non-backflow wave functions, typically yielding an additional few percent DMC correlation energy in atomic calculations and nearly 100 percent of the correlation energy in the homogeneous electron gasrios06 (). However, since backflow is too expensive to apply routinely, we have used the backflow result, similar to how corrections for all-electrons have been applied esler10 (), to shift our single determinant fixed-node energies in Fig. 3. The backflow cohesive energy is 3.403(1) eV. A complete backflow evaluation of the lattice constant might further reduce the residual differences from experiment.
We expect the backflow correction in metallic systems with atomic number higher than aluminum to be at least as large as our computed correction in aluminum. This indicates that to obtain cohesive energies to better than 0.1eV – other than by fortuitous error cancellation – backflow or other nodal optimization must be considered.
We performed additional DMC total energy calculations of bulk aluminum for atomic volumes smaller than those shown in Fig. 3 by following the same procedure of using Eq. (1) to obtain infinite-size extrapolated DMC total energies from calculations using 64 and 125 atom simulation cells with converged twist-averaged boundary conditions. A Murnaghan fit of the energy-volume DMC data was used to obtain the pressure for a range of atomic volumes. These calculated DMC pressures are shown in the upper part of Fig. 4 along with pressures calculated with DFT using GGA (PBE). At this scale the differences between the DMC and GGA pressures are not visible. The solid black curve in the lower part of Fig. 4 shows the difference in pressures between GGA and DMC. A common procedure for constructing an equation of state of a material at low temperatures is to shift the computed equilibrium lattice constant so that it coincides with the experimental equilibrium lattice constantcorrea08 (). Following a similar procedure of shifting the GGA calculated energy-volume curve so that the GGA equilibrium lattice constant of 4.046 Å coincides with the DMC equilibrium lattice constant of 4.030(1) Å yields a pressure difference between the GGA and the DMC equation of states that increases markedly at smaller atomic volumes as shown in the red dashed line in Fig.4. This demonstrates that the applying a rigid shift to DFT equation of state calculations can result in larger errors than if a shift is not applied.
Iii Defects in aluminum
The results of our calculations of point defects are presented in Table 1. The atomic positions were taken from complete structure and volume relaxed DFT calculations using GGA at zero pressure. The calculated GGA and DMC vacancy formation and migration energy agree with experiment. However, GGA no longer agrees with experiment when “surface” corrections carling00 (); sandberg02 (); mantina08 () of 0.15 eV and 0.05 eV are added to the GGA (PBE) functional producing 0.82 eV and 0.65 eV for the vacancy formation and migration energy, respectively. Previous GGA (PBE) calculations mattsson06 () of the vacancy formation energy without the “surface” correction using 4x4x4 supercells yielded values of 0.61 eV and 0.63 eV with norm-conserving and projector augmented-wave pseudopotentials, respectively. This difference with our GGA results is likely partially due to finite-size effects since we obtained 0.64 eV for the vacancy formation energy using a 4x4x4 supercell. The GGA results in Table 1 correspond to calculations using finite-size converged 7x7x7 supercells. Convergence of the GGA defect structures was established by computing the energies using , , , and supercells. The DMC results in Table 1 were done using supercells with -point centered grids of twists. The finite-size errors for our DMC calculations are likely to be small since the largest GGA energy difference among all of the defects comparing and supercells was 0.02 eV. Shown in Table 2 are GGA and DMC data demonstrating the convergence with supercell size for all the defects considered.
GGA | DMC | Experiment | |
---|---|---|---|
0.67 {0.82} | 0.67(3) erhart91 (), 0.67 hehenkamp94 (), 0.66(2) fluss78 () | ||
0.60 {0.65} | 0.62 seeger71 (), 0.61(3) erhart91 (), 0.65(6) desorbo59 () | ||
1.37 | 1.17(7)^{1}^{1}1Computed using from Ref. [erhart91, ] and from Ref. [doyama64, ] | ||
-0.03 | 0.17(5) doyama64 (), 0.20 hehenkamp94 () | ||
2.70 | 3.0 erhart91 (), 3.2(5) schilling78 () | ||
2.91 | |||
3.23 | |||
1.63 | |||
3.26 | |||
3.33 |
GGA | DMC with DFT corrections | DMC with twist averaging | ||||||
---|---|---|---|---|---|---|---|---|
supercell | ||||||||
For the nearest-neighbor divacancy with DMC we obtain a negative binding energy, -0.10 eV, which implies that two isolated vacancies are energetically preferred to a nearest-neighbor divacancy. This agrees with previous DFT calculations which also found a similar negative binding energy, -0.05 eV using uesugi03 () LDA and -0.08 eV using carling00 () GGA (PW91). Thus the disagreement between previous calculations and the original interpretation of experimental data hehenkamp94 (); fluss84 (), which gave positive binding, is likely not a result of DFT approximations. Our results are consistent with the reinterpretation carling00 () of the data.
Experimentally schilling78 () the -dumbbell was found to be the lowest energy self-interstitial in aluminum. Of the self-interstitials investigated we found that the -dumbbell has the lowest formation energy. The calculated formation energy was 2.94 eV using DMC and 2.70 eV using GGA. Our DMC value agrees with the experimental estimates of 3.0 erhart91 () and 3.2(5) eV schilling78 (). For the -dumbbell we obtain a relaxation volume of 2.25, which agrees closely with experimental estimates of 1.9(4) and 1.7(4) schilling78 (). Our GGA (PBE) result is larger than a previous GGA (PW91) result sandberg02 () of 2.43 eV. For the self-interstitials we see differences between the calculated DMC and GGA formation energies as large as 0.24 eV.
The DMC formation energy for a He substitutional defect is 1.72 eV, while the energies for He interstitials are larger than 3 eV. The ordering of these energies is consistent with sites with larger free volumes having lower energies. Previous DFT calculations yang08 () using GGA (PW91) obtained 1.53, 3.18, and 3.20 eV for He at the substitutional, octahedral, and tetrahedral site, respectively. Comparing our GGA and DMC calculations we see differences between 0.09 and 0.34 eV for the He impurity. Similar to the self-interstitials, the GGA exchange-correlation errors are larger for these defects than for the vacancy.
Relying on DFT corrections to minimize the single-particle finite-size errors instead of twist averaging, yielded poorer convergence for defect energies. For example, the difference in vacancy formation energies between 4x4x4 and 5x5x5 supercells was 0.23 eV compared with 0.01 eV with twist averaging, while the result for divacancy binding was reversed as shown in Table 2. DMC calculations in larger 6x6x6 supercells would be an order of magnitude more expensive, and may still be inferior to the twist-averaged results. The use of twist averaging is essential in metals for defects and excitations in which the fractional change in the total energy due to the presence of the defect or excitation is inversely proportional to the number of atoms in the supercell, i.e., “1/N” effects.
Iv Conclusions
In summary, DMC with twist-averaged boundary conditions can be used to obtain an accurate equation of state of aluminum. Our DMC results confirm previous DFT calculations that the nearest-neighbor divacancy is unstable in aluminum. Our calculated formation and migration energies of point defects show excellent agreement with available experiment, demonstrating that DMC can be used to obtain benchmark energetics of defects in metals and used as a baseline where no experiment is available.
Acknowledgements.
We thank J. Dubois and J. Kim for useful discussions. This material is based upon work supported as part of the CDP, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences under Award Number ERKCS99. This work performed under the auspices of the U.S. DOE by LLNL under Contract DE-AC52-07NA27344. Computing support for this work came from the LLNL Institutional Computing Grand Challenge program. Research performed at the Materials Science and Technology Division and the Center of Nanophase Material Sciences at ORNL was sponsored by the Division of Materials Sciences and the Division of Scientific User Facilities U.S. DOE.References
- (1) For a review see, W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- (2) W-K. Leung, R. J. Needs, G. Rajagopal, S. Itoh, and S. Ihara, Phys. Rev. Lett. 83, 2351 (1999).
- (3) R. Q. Hood, P. R. C. Kent, R. J. Needs, and P. R. Briddon, Phys. Rev. Lett. 91, 076403 (2003).
- (4) J. C. Grossman and L. Mitas, Phys. Rev. Lett. 79, 4353 (1997).
- (5) R. Gaudoin and W. M. C. Foulkes, Phys. Rev. B 66, 052104 (2002).
- (6) R. Gaudoin, W. M. C. Foulkes, and G. Rajagopal, J. Phys.: Condens. Matter 14, 8787 (2002).
- (7) K. Carling, G. Wahnstrom, T. R. Mattsson, A. E. Mattsson, N. Sandberg, and G. Grimvall, Phys. Rev. Lett. 85, 3862 (2000).
- (8) A. E. Mattsson, R. Armiento, P. A. Schultz, and T. R. Mattsson, Phys. Rev B 73, 195123 (2006).
- (9) L. Delczeg, E. K. Delczeg-Czirjak, B. Johansson, and L. Vitos, Phys. Rev. B 80, 205121 (2009).
- (10) P. Erhart, P. Jung, H. Schult and H. Ullmaier, Atomic Defects in Metals, Landolt-Börnstein, New Series, Group III, vol. 25 Condensed Matter (Springer-Verlag, Heidelberg, 1991).
- (11) T. Hehenkamp, J. Phys. Chem. Solids 55, 907 (1994).
- (12) M. J. Fluss, L. C. Smedskjaer, M. K. Chason, D. G. Legnini, and R. W. Siegel, Phys. Rev. B 17, 3444 (1978).
- (13) A. Seeger, D. Wolf and H. Mehrer, Phys. Status Solidi B 48, 481 (1971).
- (14) W. DeSorbo and D. Turnbull, Acta Met. 7, 83 (1959); W. DeSorbo and D. Turnbull, Phys. Rev. 115, 560 (1959).
- (15) T. Uesugi, M. Kohyama, and K. Higashi, Phys. Rev. B 68, 184103 (2003).
- (16) M. J. Fluss, S. Berko, B. Chakraborty, K. R. Hoffmann, P. Lippel, and R. W. Siegel, J. Phys. F 14, 2831 (1984).
- (17) A. C. Damask, G. J. Dienes, and V. G. Weizer, Phys. Rev. 113, 781 (1959).
- (18) H. Rajainmäki, S. Linderoth, H. E. Hansen, R. M. Nieminen, and M. D. Bentzon, Phys. Rev. B 38, 1087 (1988).
- (19) R. C. Birtcher, S. E. Donnelly, and C. Templier, Phys. Rev. B 50, 764 (1994).
- (20) D. Hamaguchi and Y. Dai, J. Nucl. Mater. 329, 958 (2004).
- (21) Y. Katoh, M. Ando, and A. Kohyama, J. Nucl. Mater. 323, 251 (2003).
- (22) R. Vassen, H. Trinkaus, and P. Jung, Phys. Rev. B 44, 4206 (1991).
- (23) L. Yang, X. T. Zu, and F. Gao, Physica B 403, 2719 (2008).
- (24) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, R. G. Hennig, Phys. Rev. Lett. 98 110201 (2007).
- (25) M. Bajdich, M. L. Tiago, R. Q. Hood, P. R. C. Kent, and F. A. Reboredo, Phys. Rev. Lett. 104 193001 (2010).
- (26) K. P. Esler, R. E. Cohen, B. Militzer, J. Kim, R. J. Needs, and M. D. Towler Phys. Rev. Lett. 104, 185702 (2010)
- (27) R. J. Needs, M. D. Towler, N. D. Drummond and P. López Ríos, J. Phys.: Condens. Matter 22, 023201 (2010).
- (28) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (29) J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992); 48, 4978 (1993).
- (30) P. Giannozzi et al., J. Phys.: Condens. Matter 21, 395502 (2009).
- (31) L. Mitàš, E. L. Shirley, and D. M. Ceperley, J. Chem. Phys. 95, 3467 (1991).
- (32) A. J. Williamson, R. Q. Hood, and J. C. Grossman, Phys. Rev. Lett. 87, 246406 (2001).
- (33) C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001).
- (34) A. M. Teweldeberhan, J. L. Dubois, and S. A. Bonev, Phys. Rev. Lett 105, 235503 (2010).
- (35) N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Phys. Rev. B 78, 125106 (2008).
- (36) F. D. Murnaghan, PNAS 30, 244 (1944).
- (37) P. López Ríos, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74, 066701 (2006).
- (38) A. A. Correa, L. X. Benedict, D. A. Young, E. Schwegler, and S. A. Bonev, Phys. Rev. B 78, 024101 (2008).
- (39) N. Sandberg, B. Magyari-Köpe, and T. R. Mattsson, Phys. Rev. Lett. 89, 065901 (2002).
- (40) M. Mantina, Y. Wang, R. Arroyave, L. Q. Chen, Z. K. Liu, and C. Wolverton, Phys. Rev. Lett. 100, 215901 (2008).
- (41) M. Doyama and J. S. Koehler, Phys. Rev. 134, A522 (1964).
- (42) W. Schilling, J. Nucl. Mater 69, 465 (1978).