Accuracy of Quantum Monte Carlo Methods for Point Defects in Solids

Accuracy of Quantum Monte Carlo Methods for Point Defects in Solids

William D. Parker1 1 Department of Physics, The Ohio State University, 191 W. Woodruff Ave., Columbus, Ohio 43210, USA
2 Department of Materials Science and Engineering, Cornell University, Ithaca, New York 14853, USA
   John W. Wilkins1 1 Department of Physics, The Ohio State University, 191 W. Woodruff Ave., Columbus, Ohio 43210, USA
2 Department of Materials Science and Engineering, Cornell University, Ithaca, New York 14853, USA
   Richard G. Hennig\Ast1,2 1 Department of Physics, The Ohio State University, 191 W. Woodruff Ave., Columbus, Ohio 43210, USA
2 Department of Materials Science and Engineering, Cornell University, Ithaca, New York 14853, USA
   1 Department of Physics, The Ohio State University, 191 W. Woodruff Ave., Columbus, Ohio 43210, USA
2 Department of Materials Science and Engineering, Cornell University, Ithaca, New York 14853, USA
XXXX, revised XXXX, accepted XXXX

Quantum Monte Carlo approaches such as the diffusion Monte Carlo (DMC) method are among the most accurate many-body methods for extended systems. Their scaling makes them well suited for defect calculations in solids. We review the various approximations needed for DMC calculations of solids and the results of previous DMC calculations for point defects in solids. Finally, we present estimates of how approximations affect the accuracy of calculations for self-interstitial formation energies in silicon and predict DMC values of , and  eV for the X, T and H interstitial defects, respectively, in a 16(+1)-atom supercell.


e-mail, Phone: +01-607-2546429



1 Introduction

Point defects, such as vacancies, interstitials and anti-site defects, are the only thermodynamically stable defects at finite temperatures [2]. The infinite slope of the entropy of mixing at infinitesimally small defect concentrations results in an infinite driving force for defect formation. As a result, at small defect concentrations, the entropy of mixing always overcomes the enthalpy of defect formations. In addition to being present in equilibrium, point defects often control the kinetics of materials, such as diffusion and phase transformations and are important for materials processing. The presence of point defects in materials can fundamentally alter the electronic and mechanical properties of a material. This makes point defects technologically important for applications such as doping of semiconductors [3, 4], solid solution hardening of alloys [5, 6], controlling the transition temperature for shape-memory alloys [7] and the microstructural stabilization of two-phase superalloys.

However, the properties of defects, such as their structures and formation energies, are difficult to measure in some materials due to their small sizes, low concentrations, lack of suitable radioactive isotopes, etc. Quantum mechanical first-principles or ab initio theories make predictions to fill in the gaps left by experiment [8].

The most widely used method for the calculation of defect properties in solids is density functional theory (DFT). DFT replaces explicit many-body electron interactions with quasiparticles interacting via a mean-field potential, i.e. the exchange-correlation potential, which is a functional of the electron density [9]. A universally true exchange-correlation functional is unknown, and DFT calculations employ various approximate functionals, either based on a model system or an empirical fit. The most commonly used functionals are based on DMC simulations [10] for the uniform electron gas at different densities, e.g., the local density approximation (LDA) [11, 12] and gradient expansions, e.g., the generalized gradient approximation (GGA) [13, 14, 15, 16, 17]. These local and semi-local functionals suffer from a significant self-interaction error reflected in the variable accuracy of their predictions for defect formation energies, charge transition levels and band gaps [18, 19]. Another class of functionals, called hybrid functionals, include a fraction of exact exchange to improve their accuracy [20, 21].

The seemingly simple system of Si self-interstitials exemplifies the varied accuracy of different density functionals and many-body methods. The diffusion and thermodynamics of silicon self-interstitial defects dominates the doping and subsequent annealing processes of crystalline silicon for electronics applications [4, 22, 23]. The mechanism of self-diffusion in silicon is still under debate. Open questions [24] include: (1) are the interstitial atoms the prime mediators of self-diffusion, (2) what is the specific mechanism by which the interstitials operate, and (3) what is the value of the interstitial formation energy? Quantum mechanical methods are well suited to determine defect formation energies. LDA, GGA and hybrid functionals predict formation energies for these defects ranging from about to  eV [25]. Quasiparticle methods such as the  approximation reduce the self-interaction error in DFT and are expected to improve the accuracy of the interstitial formation energies. Recent calculations [26] predict formation energies of about  eV in close agreement with HSE hybrid functionals [25] and previous DMC calculations [27, 25]. Quantum Monte Carlo methods provide an alternative to DFT and a benchmark for defect formation energies [28, 29].

In this paper, we review the approximations that are made in diffusion Monte Carlo (DMC) calculations for solids and estimate how these approximations affect the accuracy of point defect calculations, using the Si self-interstitial defects as an example. Section 2 describes the quantum Monte Carlo method and its approximations. Section 3 reviews previous quantum Monte Carlo calculations for defects in solids, and Section 4 discusses the results of our calculations for interstitials in silicon and the accuracy of the various approximations.

2 Quantum Monte Carlo method

Quantum Monte Carlo (QMC) methods are among the most accurate electronic structure methods available and, in principle, have the potential to outperform current computational methods in both accuracy and cost for extended systems. QMC methods scale as with system size and can handle large systems. At the present time, calculations for as many as 1,000 electrons on 1,000 processors make effective use of available computational resources [25]. Current work is under way to develop algorithms that extend the system size accessible by QMC methods to petascale computers [30].

Continuum electronic structure calculations primarily use two QMC methods [28]: the simpler variational Monte Carlo (VMC) and the more sophisticated diffusion Monte Carlo (DMC). In VMC, a Monte Carlo method evaluates the many-dimensional integral to calculate quantum mechanical expectation values. Accuracy of the results depends crucially on the quality of the trial wave function which is controlled by the functional form of the wave function and the optimization of the wave functions parameters [31]. DMC removes most of the error in the trial wave function by stochastically projecting out the ground state using an integral form of the imaginary-time Schrödinger equation.

One of the most accurate forms of trial wave functions for quantum Monte Carlo applications to problems in electronic structure is a sum of Slater determinants of single-particle orbitals multiplied by a Jastrow factor and modified by a backflow transformation:

The Jastrow factor typically consists of a low order polynomial and a plane-wave expansion in electron coordinates and nuclear coordinates that efficiently describe the dynamic correlations between electrons and nuclei. Static (near-degeneracy) correlations are described by a sum of Slater determinants. Symmetry-adapted linear combinations of Slater determinants, so-called configuration state functions (CSF), reduce the number of determinant parameters . For extended systems, the lack of size consistency for a finite sum of CSF’s makes this form of trial wave functions impractical, and a single determinant is used instead. Finally, the backflow transformation allows the nodes of the trial wave function to be moved which can efficiently reduce the fixed-node error [32]. Since the backflow transformed coordinate of an electron depends on the coordinates of all other electrons, the Sherman-Morrison formula used to efficiently update the Slater determinant does not apply, increasing the scaling of QMC to . If a finite cutoff for the backflow transformation is used, the Sherman-Morrison-Woodbury formula [33] applies and the scaling is reduced to .

Optimization of the many-body trial wave function is crucial because accurate trial wave functions reduce statistical and systematic errors in both VMC and DMC. Much effort has been spent on developing improved methods for optimizing many-body wave functions, and this continues to be the subject of ongoing research. Energy and variance minimization methods can effectively optimize the wave function parameters in VMC calculations [31, 34]. Recently developed energy optimization methods enable the efficient optimization of CSF coefficients and orbital parameters in addition to the Jastrow parameters for small molecular systems, eliminating the dependence of the results on the input trial wave function [31].

VMC and DMC contain two categories of approximation to make the many-electron solution tractable: controlled approximations, whose errors can be made arbitrarily small through adjustable parameters, and uncontrolled approximations, whose errors are unknown exactly. The controlled approximations include the finite DMC time step, the finite number of many-electron configurations that represent the DMC wave function, the basis set approximation, e.g., spline or plane-wave representation, for the single-particle orbitals such of of the trial wave function and the finite-sized simulation cell. The uncontrolled approximations include the fixed-node approximation which constraints the nodes of the wave function in DMC to be the same as the ones for the trial wave function, the replacement of the core electrons around each atom with a pseudopotential to represent the core-valence electronic interaction and the locality approximation that uses the trial wave function to project the nonlocal angular momentum components of the pseudopotential.

2.1 Controlled approximations

Time step

Diffusion Monte Carlo is based on the transformation of the time-dependent Schrödinger equation into an imaginary-time diffusion equation with a source-sink term. The propagation of the -dimensional electron configurations (walkers) that sample the wave function requires a finite imaginary time step which introduces an error in the resulting energy [35, 36].

Controlling the time step error is simply a matter of performing calculations for a range of time steps either to determine when the total energy or defect formation energy reaches the required accuracy or to perform an extrapolation to a zero time step using a low order polynomial fit of the energy as a function of time step. Smaller time steps, however, require a larger total number of steps to sample sufficiently the probability space. Thus, the optimal time step should be small enough to add no significant error to the average while large enough to keep the total number of Monte Carlo steps manageable. In addition, the more accurate the trial wave function is the smaller the error due to the time-step will be [36].

Configuration population

In DMC, a finite number of electron configurations represent the many-body wave function. These configurations are the time-independent Schrödinger equation’s analogues to particles in the diffusion equation and have also been called psips [35] and walkers [28]. To improve the efficiency of sampling the many-body wave function, the number of configurations is allowed to fluctuated from time step to time step in DMC using a branching step. However, the total number of configurations needs to be controlled to avoid the configuration population to diverge or vanish [36]. This population control introduces a bias in the energy. In practice where tested [37], few hundreds of configurations are sufficient to reduce the population control bias in the DMC total energy below the statistical uncertainty.

The VMC and DMC calculations parallelize easily over walkers. After an initial decorrelation run, the propagation of a larger number of walkers is computationally equivalent to performing more time steps. The variance of the total energy scales like

where denotes the number of walkers, the number of time steps and the auto correlation time.

Basis set

A sum of basis functions with coefficients represents the single-particle orbitals in the Slater determinant. A DFT calculation usually determines these coefficients. Plane waves provide a convenient basis for calculations of extended systems since they form an orthogonal basis that systematically improves with increasing number of plane waves that span the simulation cell. Increasing the number of plane waves until the total energy converges within an acceptable threshold in DFT creates a basis set that has presumably the same accuracy in QMC.

Since the plane wave basis functions are extended throughout the simulation cell, the evaluation of an orbital at a given position requires a sum over all plane waves. Furthermore, the number of plane waves is proportional to the volume of the simulation cell. The computational cost of orbital evaluation can significantly be reduced by using a local basis, such as B-splines, which replaces the sum over plane waves with a sum over a small number of local basis functions. The resulting polynomial approximation reduces the computational cost of orbital evaluation at a single point from the number of plane waves (hundreds to thousands depending on the basis set) to the number of non-zero polynomials (64 for cubic splines) [38]. The wavelength of the highest frequency plane wave sets the resolution of the splines. Thus, the most important quantity to control in the basis set approximation is the size of the basis set.

Simulation cell

Simulation cells with periodic boundary conditions are ideally suited to describe an infinite solid but result in undesirable finite-size errors that need correction. There are three types of finite-size errors. First, the single-particle finite-size error arises from the choice of a single -point in the single-particle Bloch orbitals of the trial wave function. Second, the many-body finite size error arises from the non-physical self-image interactions between electrons in neighboring cells. Third, the defect creates a strain field that results in an additional finite size error for small simulation cells.

The single-particle finite size error is greatly reduced by averaging DMC calculations for single-particle orbitals at different -points that sample the first Brillouin zone of the simulation cell, so-called twist-averaging [39] Alternatively, the single-particle finite size error can also be estimated from the DFT energy difference between a calculation with a dense -point mesh and one with the same single -point chosen for the orbitals of the QMC wave function.

For the many-body finite size error, several methods aim to correct the fictitious periodic correlations between electrons in different simulation cells. The first approach, the model periodic Coulomb (MPC) interaction [40], revises the Ewald method [41] to account for the periodicity of the electrons by restoring the Coulomb interaction within the simulation cell and using the Ewald interaction to evaluate the Hartree energy. The second approach is based on the random phase approximation for long wave lengths. The resulting first-order, finite-size-correction term for both the kinetic and potential energies can be estimated from the electronic structure factor [42]. The third approach estimates the many-body finite-size error from the energy difference between DFT calculations using a finite-sized and an infinite-sized model exchange-correlation functional [43]. This approach relies on the exchange-correlation functional being a reasonable description of the system, whereas the other two approaches (MPC and structure factor) do not have this restriction. The MPC and structure factor corrections are fundamentally related and often result in similar energy corrections [44].

The defect strain finite size error, can be estimated at the DFT level using extrapolations of large simulation cells. Also since QMC force calculations are expensive and still under development [45], QMC calculations for extended systems typically start with DFT-relaxed structures. Energy changes due to small errors in the ionic position as well as thermal disorder are expected to be quite small because of the quadratic nature of the minima and will largely cancel when taking energy differences for the defect energies.

2.2 Uncontrolled approximations

Fixed-node approximation

The Monte Carlo algorithm requires a probability distribution, which is non-negative everywhere, but fermions, such as electrons, are antisymmetric under exchange, and therefore any wave function of two or more fermions has regions of positive and negative value. For quantum Monte Carlo to take the wave function as the probability distribution, Anderson [35] fixed the zeros or nodes of the wave function and took the absolute value of the wave function as the probability distribution. If the trial wave function has the nodes of the ground state, then DMC projects out the ground state. However, if the nodes differ from the ground state, then DMC finds the closest ground state of the system within the inexact nodal surface imposed by the fixed-node condition. This inexact solution has an energy higher than that of the ground state.

Three methods estimate the size of the fixed node approximation: (1) In the Slater-Jastrow form of the wave function, the single-particle orbitals in the Slater determinant set the zeroes of the trial wave function. Since these orbitals come from DFT calculations, varying the exchange-correlation functional in DFT changes the trial wave function nodes and provides an estimate of the size of the fixed-node error. (2) López Ríos et al. [32] applied backflow to the nodes by modifying the interparticle distances, enhancing electron-electron repulsion and electron-nucleus attraction. The expense of the method has thus far limited its application in the literature to studies of second and third-row atoms, the water dimer and the 1D and 2D electron gases. (3) Because the eigenfunction of the Hamiltonian has zero variance in DMC, a linear extrapolation from the variances of calculations with and without backflow to zero variance estimates the energy of the exact ground state of the Hamiltonian.


Valence electrons play the most significant roles in determining a composite system’s properties. The core electrons remain close to the nucleus and are largely inert. The separation of valence and core electron energy scales allows the use of a pseudopotential to describe the core-valence interaction without explicitly simulating the core electrons. However, there is often no clear boundary between core and valence electrons, and the core-valence interaction is more complicated than a simple potential can describe. Nonetheless, the computational demands of explicitly simulating the core electrons and the practical success of calculations with pseudopotentials in reproducing experimental values promote their continued use in QMC. Nearly all solid-state and many molecular QMC calculations to date rely on pseudopotentials to reduce the number of electrons and the time requirement of simulating the core-electron energy scales.

Comparing DMC energies using pseudopotentials constructed with different energy methods (DFT and Hartree-Fock[HF]) provides an estimate of the error incurred by the pseudopotential approximation. Additionally, the difference between density functional pseudopotential and all-electron energies estimates the size of the error introduced by the pseudopotential and is used as a correction term.

Pseudopotential locality

DMC projects out the ground state of a trial wave function but does not produce a wave function, only a distribution of point-like configurations. However, the pseudopotential contains separate potentials (or channels) for different angular-momenta of electrons. One channel, identified as local, does not require the wave function to evaluate, but the nonlocal channels require an angular integration to evaluate, and such an integration requires a wave function. Mitáš et al. [46] introduced use of the trial wave function to evaluate the nonlocal components requiring integration. This locality approximation has an error that varies in sign. While there are no good estimates of the magnitude of this error, Casula [47] developed a lattice-based technique that makes the total energy using a nonlocal potential an upper bound on the ground-state energy. Pozzo and Alfè [48] found that, in magnesium and magnesium hydride, the errors of the locality approximation and the lattice-regularized method are comparably small, but the lattice method requires a much smaller time step ( Ha vs.  Ha in Mg and  Ha vs.  Ha in MgH) to achieve the same energy. Thus, they chose the nonlocal approximation.

While all-electron calculations would, in principle, make the pseudopotential and locality errors controllable, in practice, the increase in number of electrons, required variational parameters and variance of the local energy makes such calculations currently impractical for anything but small systems and light elements [49].

Energy DFT DMC Exp. Ref.
type LDA GGA Hybrid
C diamond vacancy Formation 6.98 7.51 - - 5.96(34) - [50, 51]
Migration 2.83 - - - 4.40(36) 2.3(3)
MgO Schottky defect Formation 5.97, 6.99, 6.684 - - - 7.50(53) 5 - 7 [52, 53]
self- X 3.31 3.64 4.69 4.40 5.0(2), 4.94(5) - [27, 25]
Si interstitial T Formation 3.43 3.76 4.95 4.51 5.5(2), 5.13(5) -
defect H 3.31 3.84 4.80 4.46 4.7(2), 5.05(5) -
Figure 1: DMC,  and DFT energies (in eV) for neutral defects in three materials. DMC and experimental values have an estimated uncertainty indicated by numbers in parenthesis. For the diamond vacancy, DFT-LDA and DMC include a  eV Jahn-Teller relaxation energy. LDA relaxation produced the structures and transition path so the DMC value for migration energy is an upper bound on the true value. The Schottky energy in MgO is the energy to form a cation-anion vacancy pair. DFT-LDA produces a range from - eV depending on the representation of the orbitals and treatment of the core electrons. DMC using a plane-wave basis and pseudopotentials results in a value on the upper end of the experimental range. For Si interstitial defects, DFT values of the formation energy range from  eV below up to the DMC values, depending on the exchange-correlation functional(LDA, GGA[PBE] or hybrid[HSE]), and the  values lie within the two-standard-deviation confidence level of DMC.

3 Review of previous DMC defect calculations

To date, there have been DMC calculations for defects in three materials: the vacancy in diamond, the Schottky defect in MgO and the self-interstitials in Si.

3.1 Diamond vacancy

Diamond’s high electron and hole mobility and its tolerance to high temperatures and radiation make it a technologically important semiconductor material. Diffusion in diamond is dominated by vacancy diffusion [54], and the vacancy is also associated with radiation damage [55]. Table 1 shows the range of vacancy formation and migration energies calculated by LDA [51] and DMC [50]. DMC used structures from LDA relaxation and single-particle orbitals employing a Gaussian basis. LDA pseudopotential removed core electrons. The DMC calculations predict a lower formation energy than LDA. The DMC value for the migration energy is an upper bound on the actual number since the structures have not been relaxed in DMC. Furthermore, DMC estimates the experimentally observed dipole transition and provides an upper bound on the migration energy. [50] The GR1 optical transition is not a transition between one-electron states but between spin states E and T. DMC calculates a transition energy of  eV from E to T, close to the experimentally observed value of  eV. LDA cannot distinguish these states. For the cohesive energy, DMC predicts a value of  eV in excellent agreement with the experimental result of  eV while LDA overbinds and yields  eV.

3.2 MgO Schottky defect

MgO is an important test material for understanding oxides. Its rock-salt crystal structure is simple, making it useful for computational study. Schottky defects are one of the main types of defects present after exposure to radiation, according to classical molecular dynamics simulations [56]. Table 1 shows that DMC predicts a Schottky defect formation energy in MgO at the upper end of the range of experimental values [52].

3.3 Si interstitial defects

Table 1 shows that DFT and DMC differ by up to  eV in their predictions of the formation energies of these defects [27, 25]. We compare the DMC values with our results including tests on the QMC approximations in Section 4.

Defect formation energy (eV) Jastrow factor Pseudo- Plane-wave Finite-size Reference
-body terms plane- potential cutoff energy correction
X T H e-e e-n e-e-n waves method (eV)
5.0(2) 5.5(2) 4.7(2) 16 16 0 0 LDA 245 DFT k-pt.  [27]
4.94(5) 5.13(5) 5.05(5) 5 5 0 0 HF 435 DFT k-pt.  [25]
4.9(1) 5.2(1) 4.9(1) 8 8 3 19 DF 1088 DFT k-pt.+struc. fac. (this work)
Slater-Jastrow backflow
4.5(1) 5.1(1) 4.7(1) 8 8 3 19 DF 1088 DFT k-pt.+struc. fac. (this work)
4.4(1) 5.1(1) 4.7(1) 8 8 3 19 DF 1088 DFT k-pt.+struc. fac. (this work)
Figure 2: DMC Si defect formation energies. Varying parameters and improved methods produce values for each defect that lie within two standard deviations of each other although the energetic ordering of the defects varies. All calculations use DFT-LDA to produce the orbitals in the Slater determinant.

4 Results

We specifically test the time-step, pseudopotential and fixed-node approximations for the formation energies of three silicon self-interstitial defects, the split- interstitial (X), the tetrahedral interstitial (T) and the hexagonal interstitial (H). The QMC calculations are performed using the casino [57] code. Density functional calculations in this work used the Quantum ESPRESSO [58] and WIEN2k [59] codes. The defect structures are identical to those of Batista et al.  [25]. The orbitals of the trial wave function come from DFT calculations using the LDA exchange-correlation functional. The plane-wave basis set with a cutoff energy of  eV ( Ha) converges the DFT total energies to  meV. A 777 Monkhorst-Pack -point mesh centered at the L-point (0.5,0.5,0.5) converges the DFT total energy to  meV. A population of walkers ensured that the error introduced by the population control is negligible small. Due to the computational cost of backflow, we perform the simulations for a supercell of 16(+1) atoms and estimate the finite-size corrections using the structure factor method [42]. The final corrected DMC energies for the X, T and H defects are shown in the bottom line of Table 2.

Figure 3: DMC total energies with varying (imaginary) time steps for bulk silicon and the X defect. The error due to the finite-sized time step is smaller than the statistical uncertainty in the total energies for values from to  Ha. Note that these energies include no finite-size or pseudopotential corrections and thus differ in value from those in Figure 4.

4.1 Time step

Figure 3 shows the total energies of bulk silicon and the X defect as a function of time step in DMC. A time step of  Ha reduces the time step error to within the statistical uncertainty of the DMC total energy.

4.2 Pseudopotential

In our calculations, a Dirac-Fock (DF) pseudopotential represents the core electrons for each silicon atom [60, 61, 62]. To estimate the error introduced by the pseudopotential, we compare the defect formation energies in DFT using this pseudopotential with all-electron DFT calculations using the linearized augmented plane-wave method [59]. This comparison gives corrections of , and  eV for the H, T and X defects respectively.

Figure 4: DMC total energies calculated with and without backflow transformation and linearly extrapolated to zero variance. The backflow transformation reduces the total energies by 0.19(3), 0.62(5), 0.25(5) and 0.35(5) eV for bulk (solid line) and the X (dash-dot line), T (dashed line) and H (dotted line) defects respectively. The extrapolation to zero variance only slightly reduces the total energies by about 0.1 eV. The reduction of the total energy from the backflow transformation indicates the size of the errors due to the fixed-node and pseudopotential locality approximation.

4.3 Fixed-node approximation

The fixed-node approximation is the main source of error of diffusion Monte-Carlo calculations. To estimate the size of the error introduced by the fixed-node approximation, we perform the calculations using the backflow transformation which allows the nodes of the trial wave function to be moved and reduces the fixed-node error [32]. We estimate the error due to the fixed-node approximation by performing calculations with and without the backflow transformation and by extrapolating the resulting defect formation energies to zero variance.

Applying the backflow transformation to electron coordinates using polynomials of electron-electron, electron-nucleus and electron-electron-nucleus separations, we include polynomial terms to eighth-order for each spin type in electron-electron separation, to sixth-order in electron-nucleus separation and to third-order for each spin type in electron-electron-nucleus separation. Figure 4 shows the linear extrapolation of the DMC energies for the Slater-Jastrow and Slater-Jastrow-backflow trial wave function to zero variance. The total energy decrease for the bulk and interstitial cells due to the backflow transformation ranges from 0.20(5) to 0.62(5) eV. The backflow transformation results in a significantly improved nodal surface of the trial wave function which is reflected in the reduced the variance of the local energy.

Table 2 list the Si interstitials formation energies in DMC for the Slater-Jastrow, Slater-Jastrow-backflow wave function and the extrapolation. Applying the backflow transformation reduces the formation energies for the X, T and H interstitial by 0.42(5), 0.05(5) and 0.15(5) eV, respectively. The linear extrapolation provides a simple estimate of the remaining fixed node error. The extrapolation lowers the interstitial formation energy by a negligible amount of 0.06(10), 0.00(10), 0.07(10) eV for the X, T and H interstitial, respectively. The resulting Si interstitial formation energies are 4.4(1), 5.1(1) and 4.7(1) eV for the X, T and H interstitial, respectively, in close agreement with recent  [26] and previous HSE and DMC calculations [27, 25].

5 Conclusion

QMC methods present an accurate tool for the calculation of point defect formation energies, provided care is taken to control the accuracy of all the underlying approximations. Including corrections for the approximations yields DMC values for the Si interstitial defects on par with and hybrid-functional DFT calculations. While backflow transformation and zero-variance extrapolation to remove the fixed-node error modify the energies slightly, further work remains to carefully control for finite-size effects known to plague defect supercell calculations.


The work was supported by the U.S. Department of Energy under Contract No. DE-FG02-99ER45795 and DE-FG05-08OR23339 and the National Science Foundation under Contract No. EAR-0703226. This research used computational resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, the National Center for Supercomputing Applications under grant DMR050036, the Ohio Supercomputing Center and the Computation Center for Nanotechnology Innovation at Rensselaer Polytechnic Institute. We thank Cyrus Umrigar and Ann Mattsson for helpful discussion.


  • [1]
  • [2] R. J. D. Tilley, Defects in Solids (Wiley, 2008).
  • [3] P. M. Fahey, P. B. Griffin, and J. D. Plummer, Rev. Mod. Phys. 61(2), 289–384 (1989).
  • [4] D. J. Eaglesham, P. A. Stolk, H. J. Gossmann, and J. M. Poate, Applied Physics Letters 65(18), 2305–2307 (1994).
  • [5] L. M. Pike, Y. A. Chang, and C. T. Liu, Acta Materialia 45(9), 3709 – 3719 (1997).
  • [6] J. H. Zhu, L. M. Pike, C. T. Liu, and P. K. Liaw, Acta Materialia 47(7), 2003 – 2018 (1999).
  • [7] K. Otsuka and X. Ren, Intermetallics 7(5), 511 – 528 (1999).
  • [8] L. Pelaz, L. A. Marqués, M. Aboy, P. López, and I. Santos, Eur. Phys. J. B 72(3), 323–359 (2009).
  • [9] R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61(3), 689–746 (1989).
  • [10] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45(7), 566–569 (1980).
  • [11] J. P. Perdew and A. Zunger, Phys. Rev. B 23(10), 5048–5079 (1981).
  • [12] J. P. Perdew and Y. Wang, Phys. Rev. B 45(23), 13244–13249 (1992).
  • [13] J. P. Perdew, Unified Theory of Exchange and Correlation Beyond the Local Density Approximation, in: Electronic Structure of Solids ’91, edited by P. Ziesche and H. Eschrig, (Akademie Verlag, Berlin, 1991), pp. 11–20.
  • [14] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77(18), 3865–3868 (1996).
  • [15] R. Armiento and A. E. Mattsson, Phys. Rev. B 72(8), 085108 (2005).
  • [16] Z. Wu and R. E. Cohen, Phys. Rev. B 73(23), 235116 (2006).
  • [17] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100(13), 136406 (2008).
  • [18] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, J. Chem. Phys. 123(17), 174101 (2005).
  • [19] R. M. Nieminen, Model. Sim. Mat. Sci. Eng. 17(8), 084001 (2009).
  • [20] A. D. Becke, J. Chem. Phys. 98(7), 5648–5652 (1993).
  • [21] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118(18), 8207–8215 (2003).
  • [22] D. A. Richie, J. Kim, S. A. Barr, K. R. A. Hazzard, R. Hennig, and J. W. Wilkins, Phys. Rev. Lett. 92(4), 045501 (2004).
  • [23] Y. A. Du, R. G. Hennig, T. J. Lenosky, and J. W. Wilkins, Eur. Phys. J. B 57(3), 229–234 (2007).
  • [24] R. Vaidyanathan, M. Y. L. Jung, and E. G. Seebauer, Phys. Rev. B 75, 195209 (2007).
  • [25] E. R. Batista, J. Heyd, R. G. Hennig, B. P. Uberuaga, R. L. Martin, G. E. Scuseria, C. J. Umrigar, and J. W. Wilkins, Phys. Rev. B 74(12), 121102 (2006).
  • [26] P. Rinke, A. Janotti, M. Scheffler, and C. G. V. de Walle, Phys. Rev. Lett. 102(2), 026402 (2009).
  • [27] W. K. Leung, R. J. Needs, G. Rajagopal, S. Itoh, and S. Ihara, Phys. Rev. Lett. 83(12), 2351–2354 (1999).
  • [28] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73(1), 33–83 (2001), Sections V. and VI. contrast QMC and DFT results. Section X.E. discusses scaling with computer time. Section III. introduces VMC and DMC.
  • [29] R. Needs, Quantum Monte Carlo Techniques and Defects in Semiconductors, in: Theory of Defects in Semiconductors, edited by D. Drabold and S. Estreicher, , Topics Appl. Physics Vol. 104 (Springer Verlag, Berlin, Heidelberg, 2006), p. 141.
  • [30] K. P. Esler, J. Kim, D. M. Ceperley, W. Purwanto, E. J. Walter, H. Krakauer, S. Zhang, P. R. C. Kent, R. G. Hennig, C. Umrigar, M. Bajdich, J. Kolorenc, L. Mitas, and A. Srinivasan, J. Phys. Conf. Ser. 125, 012057 (15pp) (2008).
  • [31] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett 98(11), 110201 (2007).
  • [32] P. López Ríos, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74(6), 066701 (2006).
  • [33] G. H. Golub and C. F. V. Loan, Matrix Computations, 3rd edition, (The Johns Hopkins University Press, Baltimore, 1996), chap. 2, p. 51.
  • [34] C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60(17), 1719–1722 (1988).
  • [35] J. B. Anderson, J. Chem. Phys. 63(4), 1499–1503 (1975).
  • [36] C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. Phys. 99(4), 2865–2890 (1993).
  • [37] D. Alfè, M. J. Gillan, M. D. Towler, and R. J. Needs, Phys. Rev. B 70(21), 214102 (2004).
  • [38] D. Alfè and M. J. Gillan, Phys. Rev. B 70(16), 161101 (2004).
  • [39] C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64(1), 016702 (2001).
  • [40] A. J. Williamson, G. Rajagopal, R. J. Needs, L. M. Fraser, W. M. C. Foulkes, Y. Wang, and M. Y. Chou, Phys. Rev. B 55(8), R4851–R4854 (1997).
  • [41] P. P. Ewald, Annalen der Physik 369(3), 253–287 (1921).
  • [42] S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Phys. Rev. Lett. 97(7), 076404 (2006).
  • [43] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100(12), 126404 (2008).
  • [44] N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Phys. Rev. B 78(12), 125106 (2008).
  • [45] A. Badinski, P. D. Haynes, J. R. Trail, and R. J. Needs, Journal of Physics: Condensed Matter 22(7), 074202 (2010).
  • [46] L. Mitáš, E. L. Shirley, and D. M. Ceperley, J. Chem. Phys. 95(5), 3467–3475 (1991).
  • [47] M. Casula, Phys. Rev B 74(16), 161102 (2006).
  • [48] M. Pozzo and D. Alfè, Phys. Rev. B 77(10), 104103 (2008).
  • [49] K. P. Esler, R. E. Cohen, B. Militzer, J. Kim, R. J. Needs, and M. D. Towler, Phys. Rev. Lett. 104(18), 185702 (2010).
  • [50] R. Q. Hood, P. R. C. Kent, R. J. Needs, and P. R. Briddon, Phys. Rev. Lett. 91(7), 076403 (2003).
  • [51] J. Shim, E. K. Lee, Y. J. Lee, and R. M. Nieminen, Phys. Rev. B 71(3), 035206 (2005).
  • [52] D. Alfè and M. J. Gillan, Phys. Rev. B 71(22), 220101 (2005).
  • [53] C. A. Gilbert, S. D. Kenny, R. Smith, and E. Sanville, Phys. Rev. B 76(18), 184103 (2007).
  • [54] J. Bernholc, A. Antonelli, T. M. Del Sole, Y. Bar-Yam, and S. T. Pantelides, Phys. Rev. Lett. 61(23), 2689–2692 (1988).
  • [55] A. T. Collins and I. Kiflawi, Journal of Physics: Condensed Matter 21(36), 364209 (2009).
  • [56] B. P. Uberuaga, R. Smith, A. R. Cleave, G. Henkelman, R. W. Grimes, A. F. Voter, and K. E. Sickafus, Phys. Rev. B 71(10), 104102 (2005).
  • [57] R. J. Needs, M. D. Towler, N. D. Drummond, and P. López Ríos, J. Phys. Cond. Matt. 22(2), 023201 (15pp) (2010).
  • [58] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Cond. Matt. 21(39), 395502 (19pp) (2009).
  • [59] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, WIEN2K, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, Techn. Universität Wien, Austria, 2001).
  • [60] J. R. Trail and R. J. Needs, J. Chem. Phys. 122(1), 014112 (2005).
  • [61] J. R. Trail and R. J. Needs, J. Chem. Phys. 122(17), 174109 (2005).
  • [62] J. R. Trail and R. J. Needs, CASINO pseudopotential library,
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description