Quantum Monte Carlo Calculations of Catalytic Energy Barriers in a

Article pubs.acs.org/JPCC

Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

Quantum Monte Carlo Calculations of Catalytic Energy Barriers in a Metallorganic Framework with Transition-Metal-Functionalized Nodes Anouar Benali,† Ye Luo,† Hyeondeok Shin,† Dale Pahls,‡,§ and Olle Heinonen*,‡ †

Computational Sciences Division, Argonne National Laboratory, Argonne, Illinois 60439, United States Materials Science Division, Argonne National Laboratory, Lemont, Illinois 60439, United States § Department of Chemistry, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States Downloaded via UNIV OF WINNIPEG on July 12, 2018 at 21:52:29 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: We have investigated electronic energy barriers for ethylene hydrogenation and C−H bond activation in transition-metal-functionalized Zr-based nodes in the NU-1000 metal−organic framework using quantum Monte Carlo (QMC) simulations. We compare our QMC results with results using density functional theory schemes. The density functional theory calculations can have errors up to 20 kcal/mol compared to the QMC calculations, underestimating the energy barrier for some reactions and overestimating the barrier for others. The hybrid functionals PBE0 and HSE06 generally perform best when compared with QMC, although they can still differ from the best QMC by 10 kcal/mol.

■

exchange−correlation (xc) functional.15 A well-known flaw in DFT is that band gaps in semiconductors and insulators, and highest occupied molecular orbital (HOMO)−lowest unoccupied molecular orbital (LUMO) gaps in molecular systems, are typically underestimated in simpler approximations such as the generalized gradient approximation (GGA)16 of the xc functional. While more complicated, and computationally much more expensive, approximations such as meta-GGA17 and hybrid functionals18−20 can give much better result for band gaps and HOMO−LUMO gaps, even these approximations do not achieve chemical accuracy but can have errors of the order of 10 kcal/mol. In particular, molecular systems containing transition metal atoms with partially filled 3d-shells can be especially problematic for DFT-based methods because of a poor or inaccurate description of the electronic correlations. In recent years, quantum Monte Carlo (QMC) methods, and in particular fixed-node diffusion Monte Carlo (FNDMC),21 have demonstrated a great success in describing the properties of large sets of molecules22−28 and solids,29−37 with accuracies often at or better than chemical accuracy of 1 kcal/ mol. By solving explicitly the many-body Schrödinger equation through stochastic sampling, the computational effort of FNDMC scales as ∼N3−N4, with N the number of electrons, while maintaining an almost linear scaling with the number of computer cores, making it ideal for Leadership class supercomputers. DMC achieves its favorable scaling by projecting the ground state of a trial wavefunction ΨT onto an imaginary-

INTRODUCTION A challenging and fundamental problem in heterogeneous catalysis is to design supported catalysts with high reactivity and selectivity. Metallorganic frameworks (MOFs)1−4 are porous frameworks composed of organic linkers and inorganic nodes and can provide a support system for catalytic sites. MOFs have recently been the focus of intense research as they allow both for ligand design as well as for atomically precise design and functionalization of the inorganic nodes for high catalytic reactivity.5−13 A critical step in the process of directed, mechanistic design of supported heterogeneous catalysis is the ability to calculate and model catalytic processes with chemical accuracy, ∼1 kcal/mol, especially for reaction-limiting steps. An essential component in such calculations is to accurately calculate the electronic (0 K) energy barrier from the reactant (R) to the transition state (TS), as well as from the TS to the product state (P). While quantum chemical methods, such as coupled clusters or configuration-interaction methods, can achieve such accuracy, the methods are very expensive as they scale as N7 or N!, with N the number of electrons or states. This computational expense limits the practical applicability of quantum chemical methods; this is particularly true for functionalized nodes in MOFs because of the rather large number of atoms that have to be included in such calculations. Therefore, the practical computational workhorse is density functional theory (DFT).14,15 DFT has the advantage of being rather efficient, scaling as N3, and calculations with systems composed of 103 atoms can readily be done on large computer clusters or at leadership-class facilities. However, DFT methods suffer from uncontrolled inaccuracies. These stem from the fact that electronic interactions beyond classical Hartree interactions are included only approximately in the © XXXX American Chemical Society

Received: March 10, 2018 Revised: June 8, 2018 Published: June 28, 2018 A

DOI: 10.1021/acs.jpcc.8b02368 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Ball and stick model of the R, TS, and P states from left to right for ethylene hydrogenation with Ni in silver, Zr in green, C in brown, O in red, and H in white. Prior to reaching the R state, three OH and two H2O are bound to the Ni. One OH group is released by exposure to H2 and replaced by H bound to the Ni. The complex is then activated and binds an ethylene group to the Ni with the concomitant release of two H2O to reach the R state. The lone H bound to the Ni is then transferred through the TS state and binds to the ethylene molecule in the P state.

the electron−electron distance, electron−ion distance, and a cutoff radius. While multiple choices of functions are possible, we use one-dimensional B-spline functions with a range of optimizable parameters M. The expectation value ⟨ΨT|H|ΨT⟩ is then minimized by varying the parameters M in the Jastrow factor through a variational Monte Carlo (VMC) preliminary procedure. The choice of the single-particle orbitals defining the Slater determinants becomes the only uncontrolled parameter in ΨT within the FN approximation. These orbitals can have any form, as long as they are continuous and usually come from Hartree−Fock (HF), DFT, or quantum-chemical methods such as complete active space self-consistent field (CASSCF) or full configuration interaction (CI) methods. For small molecules, multi-Slater determinant expansions from a CASSCF40 or selected CI28,41 have shown that a systematic improvement of ΨT is possible, but these approaches are not yet possible for large systems. For solids, the QMC nodal surface can be optimized by performing DMC calculations as a function of a parameter, such as the U-parameter in the socalled DFT + U methods30,31 or by modifying the level of exact exchange in hybrid functionals.42 The value of the parameter at which the DMC energy is minimized provides the optimal QMC nodal surface within that parameter space. Such optimizations typically give very good results for ground states of transition metal oxides,30,31,43 which have strong correlations and are beyond standard DFT methods, such as the GGA. A second usual approximation is to use pseudopotentials (PPs) to represent the interactions with ions and core electrons. The DMC PPs are not the same as those typically used in DFT calculationsthe DMC PPs are much “harder” (contain fewer core electrons) and are constructed to represent the many-body physics that results from interactions with core electrons. FN-DMC has been used in recent studies (see, e.g., refs44,46) to calculate reaction energy barriers in molecular systems with very good results. Zhou and Wang44 studied the 19 H-transfer reactions in the HTBH38/08 database. Using single Slater determinant trial wavefunctions based on the self-consistent DFT LDA Kohn−Sham orbitals, they found an error of about 1.6 kcal/mol or less in the calculated barrier heights when compared to “gold standard” coupled cluster CCSD(T) calculations. They also found that PPs introduced negligible errors compared to all-electron calculations. Krongchon, Busemeyer, and Wagner45 tested the accuracy of DMC forward and reverse barrier heights for 19 nonhydrogen transfer reactions on a database from Peverati and Truhlar.47

time diffusion equation that is solved through stochastic sampling. However, to maintain an antisymmetric electronic wavefunction, the nodal surface [the (3N − 1)-dimensional hypersurface where the wavefunction is zero] in the initial trial wavefunction is preserved. This is the fixed node approximation, and it constitutes the only uncontrolled approximation in the method. However, because DMC satisfies the variational Rayleigh−Ritz theorem, it is guaranteed that the FN-DMC energy, EFN, yields a strict upper bound to the true ground state energy, and if the QMC nodal surface is the true nodal surface, the DMC energy is the exact true ground state energy, E0, so EFN ≥ E0. The challenge then becomes one of how to optimize the nodal surface. In practice, the most common trial wavefunction ΨT has the form of a single- or multireference Slater−Jastrow function, where the Slater determinant is antisymmetric and imposes the nodal surface, and the Jastrow function is a correlation term that is symmetric under exchange of electrons ND

ΨT =

∑ ciDi↑Di↓ e J

(1)

i=1

where ND is the number of the spin-dependent Slater determinants Di, ci are their corresponding coefficients (weights), and J is the Jastrow function. The Jastrow function modifies the local structure of ΨT and takes into account many-body correlations. It also significantly improves the representation of the many-body wavefunction by reducing the statistical variance of the local energy and also by improving the quality of the DMC projection operator.38,39 The Jastrow function is a many-body expansion with J=

∑ V1(ri) + i

+...

1 2

∑ V2(ri, rj) + i,j

1 6

∑ V3(ri, rj, rk) i ,j,k

(2)

where ri are the positions of the N electrons. The functions V1, Vn, and so on correspond to n-body terms describing correlations. In practice, and in this study, typically only one-, two-, and three-body Jastrow functions are used. The one-body term is approximated as a sum over atom-centered swave type functions that depend on the local ionic species; the two-body term is approximated as a spin-dependent liquid-like factor (the electron−electron term), optionally with a second factor that additionally depends on the ionic coordinates (the electron−electron−ion term). The three-body electronic term is a three-dimensional polynomial function that depends on B


Article


Figure 2. Ball and stick model of the 2-Zn (left) and 4-Zn (center) R states with Cu in blue, Zn in silver, Zr in green, C in brown, O in red, and H in white. The corresponding 2-Ni and 4-Ni states are analogous with Ni replacing Zn. Prior to activation, the complex is oxidized by removing a hydrogen atom from the bridging OH group bound to the M atom (in 2-M) or the Cu atom (in 4-M). Upon introduction of a methane molecule, the transition state (right) is activated after which a methyl radical, on top between the benzoate rings, is generated.

avoids self-interactions, and DFT xc-energies provide a good approximation to the true xc-energy when calculated on the HF density. The HF-DFT method has been shown to give considerable improvement on DFT results for the energy difference between high-spin and low-spin states of the socalled spin-crossover complexes.22 We find that not to be the case here: the HF-DFT energy barriers are no better than conventional DFT barriers. We also find that the lack of a systematic method to optimize the QMC nodal surface of the trial wavefunction ΨT (which is then inherited throughout the VMC and DMC calculations) is a bit of a problem in the QMC calculations. As mentioned above, there are multiple routes to optimize the QMC nodal surface of a small molecule (multireference Slater determinant, optimizing the value of U in DFT + U, etc...). However, in the present work, the molecular systems are so large and the DMC calculations so expensive that such a parameter-sweep (e.g., by placing a Uparameter on the Ni 3d orbital) is not feasible. Instead of a parameter sweep, we used two different nodal surfaces in each barrier calculation. The DMC energies for the two nodal surfaces were typically within two standard deviations of one another with the DMC energy based on the DFT local spin density approximation (LSDA) xc functional giving lower energy. In view of earlier work establishing the accuracy of DMC energy barrier, especially the studies by Zhou and Wang44 and by Krongchon, Busemeyer, and Wagner,45 this gives us confidence that the DMC results we obtain with error bars based on the standard deviations of the DMC energies bracket the true energy barriers. The one exception, as we will discuss later, is the C−H bond activation with the Zr6 nodes functionalized with a Cu and a Zn in a particular configuration, the 2-Zn configuration. This system appears to have degenerate or near-degenerate spin multiplets which may make it a multireference system for which the true nodal surface is not accessible by our single Slater determinant trial wavefunctions.

Their calculations used PPs. They examined the error introduced by the fixed node approximation using single Slater determinant trial wavefunctions obtained from the HF and DFT calculations using the Perdew−Burke-Ernzerhof (PBE)16 and PBE019 xc functionals. Important conclusions from that study are first of all that the error introduced by the fixed node approximation is relatively robust with respect to the nodal surface chosen by these three methods, in particular the results based on the PBE and PBE0 nodal surfaces were statistically indistinguishable, and secondly that the errors were no more than about 1.5 kcal/mol compared to the lowest-energy DMC values. We are here investigating energy barriers in catalytic reactions of transition-metal-functionalized nodes in a particular MOF, NU-1000, composed of Zr6(μ3-O)4(μ3O)4(H2O)4(OH)4 nodes (hereafter, referred to as Zr6) and tetratopic 1,3,6,8-(p-benzoate)pyrene (TBAPy4−) linkers.48 Previous studies11 have shown that the Zr6 nodes can be functionalized by one or more transition metal atoms, such as Ni, Zn, and Cu, that may provide good catalytic selectivity and reactivity. We focus here on (i) ethylene hydrogenation (Figure 1) and (ii) C−H bond activation (Figure 2). Our calculations are based on previous studies that identified the reaction-limiting R, TS, and P states and also reduced the functionalized NU-1000 to more manageable clusters.11,12 We use these previous studies as starting points and use the coordinates of the R, TS, and P states obtained in those studies as input to our DFT and QMC calculations. Our aims are to assess the accuracy of DFT-based calculations with the DMC results as benchmarks and to investigate whether or not there are discernible trends in the DFT calculations that may make them more useful than the accuracy of single calculations would allow, e.g., by error cancellations in differences of energies. Our results show that DFT calculations of energy barriers can have errors in the tens of kcal/mol compared to the DMC results, with hybrid functionals such as PBE0 typically performing better than GGA or meta-GGA functionals. More troubling is the lack of a trend in our work: in (i), ethylene hydrogenation, the energy barriers are all underestimated by DFT calculations, while in (ii), C−H bond activation, they are overestimated. In addition to DFT calculations with a range of xc functionals, we also did some limited calculations using the so-called density-corrected DFT, or HF-DFT, method, in which HF density for the system is used as input to calculations of the xc-energy in DFT;49 this is then added to the kinetic, potential, and Hartree energy. The idea is that the HF density may provide a better density than DFT, especially for some molecular systems, as HF strictly

■

METHODS For part (i), the ethylene hydrogenation, the calculations were based on a cluster model containing a single Zr6 node functionalized with a Ni atom, and the linkers were truncated to formate groups. Here, the R, TS, and P states contained 77 atoms. The geometries and how they were obtained are detailed in ref 11. Our calculations are based on singlet states 5 (R) and 6 (P) and the singlet TS 5 → 6 in ref 11. For the DFT calculations, we use the all-electron FHI-aims code50 with the “light” default basis, radial mesh and angular momentum expansions and with zero-order relativistic corrections (zora). C


Article

The Journal of Physical Chemistry C For the HSE0620 and PBE019 xc functionals, we also expanded the “light” default basis set to include the “second tier” basis functions for all elements (with the “hydro 5 g 10.4” state commented out for Ni and Zr). The difference in calculated energy barriers between using the “light default” and expanded “light default” basis sets was no more than 0.01 eV or 0.2 kcal/ mol. We performed calculations within the LSDA15 using the Perdew−Wang parametrization,51 the PBE16 GGA functional, the M06-L52 and TPSS17 meta-GGA functionals, as well as the B3LYP,18 PBE0,19 and HSE0620 hybrid functionals, in addition to HF calculations. For the HF-DFT calculations, we first converged the system at the HF level and then did a singleshot non-selfconsistent calculation of the TPSS xc energy using the HF density. The HF-DFT energy was then obtained by adding the HF kinetic energy, HF total electrostatic energy (both Hartree electron−electron interactions and electron−ion interactions) to Exc,TPSS(HF). We also include the M06-L results for the electronic energy barriers from ref 11 based on Gaussian0953 with the Def2-SVP basis set for H, C, and O atoms, Def2TZVP for Ni and Zr, and SDD PP for Zr. For comparison, we also performed FHI-aims calculations using the M06-L xc functional with the default “tight” setting (includes second tier basis functions as well as refined radial mesh and higher-order angular momentum expansions) with and without zora. The difference between the “light” and “tight” energy barriers with zora were less than 0.015 eV (0.35 kcal/mol), and the difference between the “tight” energy barriers without zora and the Guassian09 energy barriers were less than 0.017 eV (0.39 kcal/mol). For part (ii), we based our calculations on the R and TS structures detailed in ref 12 with the pyrene linkers replaced by benzoate and with the Zr6 node functionalized by two transition metals, Cu and Ni or Cu and Zn. We chose these as ref 12 indicated that these structures have low activation energies for C−H bond activation and also because the Cu−Ni functionalized node is computationally more complex than the homobimetallic node functionalized with two Cu atoms (which had lower activation energy according to ref 12 than the Cu−Ni nodes). Because the metals can be supported by inequivalent positions at the μ3O and μ3OH positions, we considered both possibilities for each set of Cu−Ni and Cu− Zn. We calculated only the electronic energy difference between the R and TS, which we denote 2-M and 4-M with “M” Ni or Zn and TS2-M and TS4-M, respectively. The detailed geometries of the R and TS and how they were obtained are described in ref 12. On the basis of this reference, we considered spin quintets for the nodes containing Ni and spin triplets for nodes with Zn, although we did some checks with other spin states. The R and TS cluster models contained 153 and 158 atoms, respectively. We performed DFT calculations using FHI-aims with the default “tight” settings using the LSDA15 with the Perdew−Wang parametrization,51 PBE GGA,16 M11-L, and TPSS meta-GGA,17,54 and B3LYP and PBE019 hybrid functionals. We also include the M06-L results from ref 12 based on Gaussian09 with Stuttgart effective core potentials and associated basis sets for Ni, Cu, Zn, and Zr and 6-31G(d) for H, C, and O. The structures both for the ethylene hydrogenation and the C−H bond activation were optimized using M06-L.11,12 We also re-optimized, using the PBE0 xc functional and the “light” default settings and basis sets plus second tier basis sets, the R and TS structures in the ethylene hydrogenation keeping the carbon atoms fixed (as they represent the benzoate linkers that

are stiff), the ethylene group, and also the lone hydrogen that is initially bonded to the Ni. The differences between energy using the PBE0 xc functional on the PBE0-optimized structures and the energy using the PBE0 xc functional on the M06-L-optimized structures were less than 4.4 kcal/mol, but the energy difference between the R and TS structures changed by less than 0.05 kcal/mol. We also optimized the R structure with the “light” default settings and the PBE0 xc functional keeping only the carbon atoms fixed. This resulted in a change of total energy of only 1.2 kcal/mol. These reoptimizations indicate that energy differences between the R and TS states (and presumably between TS and P states) change very little when the structures are re-optimized because structures optimized using meta-GGA or hybrid functionals differ very little (see, however, later discussion on the 2-Zn node), and error cancellations further reduce energy differences. This is also consistent with evaluations of energy differences between high- and low-spin ferrocene complexes,22 where the energy difference was practically independent of the combination of xc functionals used for optimization and energy evaluation. All our QMC calculations were carried out using the QMCPACK code.55−58 The antisymmetric part of the trial wavefunction was generated using the “pwscf” code within the Quantum Espresso package,59 with both DFT LSDA and GGA-PBE xc functionals. We used plane wave (PW) orbitals to avoid any basis-set convergence issues that may arise with the use of Gaussian orbitals. With the choice of PW orbitals for molecular systems spurious effect stemming from the periodicity of the orbitals have to be removed by increasing the size of the simulation box until reaching convergence of the energy as a function of the box size. Using LSDA, the energy of the system converged with a box size of 30 × 30 × 30 Å3. Because all the studied systems contain transition metals, core electrons were approximated with PPs developed specifically for QMC, while reproducing other DFT results when used in DFT. Because QMC-PPs are “hard” PPs, a kinetic energy cutoff of 400 Ry was used in all our calculations. This cutoff corresponds to the highest converged cutoff among all the atoms in the molecules (in this case, Ni). Norm-conserving PPs for this study (H, C, N, O, Cu, Zn, Zr, and Ni) were generated with a PW basis-set using the OPIUM package.60 The accuracy of the Ni, Zn, Cu, and Zr PPs was already demonstrated in previous DMC studies.37,61−63 One-, two-, and three-body Jastrow factors were used to describe the many-body correlations and to improve the overall quality of the trial wavefunction and computational efficiency. The parameters of the Jastrow wavefunctions were optimized prior to the DMC calculations through VMC. Because DMC projects the ground state onto imaginary time, we used a shorttime approximation with time-steps of Δτ = 0.001. To test the accuracy of this approximation, we reduced the time-step to Δτ = 0.0008, resulting in negligible energy differences and therefore justifying the choice of the time-step. We did 1500 blocks of statistics using about 65k walkers, with 10 Monte Carlo steps per block after discarding the first 300 blocks. As mentioned in the introduction, the only uncontrolled approximation (beyond the use of PPs) is the quality of the QMC nodal surface inherited from DFT. Therefore, all DMC calculations were carried with both LSDA and GGA-PBE orbitals. Ideally, using more functionals would ensure selecting the best QMC nodal surface for our calculations,44,45 but due D


Article


results for the hybrid functionals B3LYP, PBE0, HSE06, and the meta-GGA M06-L are similar, while the PBE and TPSS values are lower. As a sanity check, the QE(PBE) and FHIaims PBE results are in good agreement with each other. The DMC results clearly show that the choice of QMC nodal surface makes some difference, although the difference between the DMC(LSDA) and DMC(PBE) energy differences are within two standard deviations of one another. In any case, the differences between the DMC(LSDA) and DMC(PBE) energy differences underscore the need to be able to more systematically optimize the QMC nodal surface for these kinds of calculations. Because the DMC(LSDA) energies for the R, TS, and P states were all lower than the DMC(PBE) ones (see Table 1), in some cases by more than two standard deviations,

to the very high cost of the calculations and the obtained results, we limited our study to these two functionals. It is important to mention that when studying molecular systems with transition metals and Gaussian orbitals, the choice of the basis set is often of paramount importance to ensure quality of the results, even though DMC significantly removes the dependence on the size of the basis set.64 In this study, we circumvented risks associated with the choice of basis set and used PWs (infinite basis set). Because of the very large PW cutoff, the memory needed to store all the orbitals on every compute core reaches almost 1 TB, which is much beyond what most supercomputers can handle. For this reason, all the orbitals were transformed into a hybrid representation,65 recently added in QMCPACK, to save memory without compromising the accuracy. This method expresses orbitals in spherical harmonics near the atomic centers while leaving them in B-splines in the interstitial regions. As a result, the memory requirement decreased from 1 TB to 110 GB, which fits most machines. In our calculations, we chose 1.2 Bohr radii for the spherical regions around every element and set the largest angular momentum to 4 (7) for C (H) and 5 for O, Zr, Cu, Zn, and Ni atoms. Computational data, such as input files and relevant output are curated on the Materials Data Facility, https:// materialsdatafacility.org/, DOI: 10.18126/M2J06G.

Table 1. DMC Energies in kcal/mol, Relative to the DMC(LSDA) Energy for the R, TS, and P States for Ethylene Hydrogenation Based on LSDA and PBE Nodal Surfaces ±1 Standard Deviation

■

state

LSDA

PBE

R TS P

0 ± 1.6 14.7 ± 1.5 −10.8 ± 1.6

5.6 ± 1.2 15.8 ± 1.1 −3.2 ± 1.1

we can invoke the variational principle of DMC and confidently assert that the DMC(LSDA) energies are more accurate than the DMC(PBE) ones. We can see from the figure that all DFT-based calculations vastly underestimate the energy differences, both for R → TS as well as for TS → R. The results are very much different for the 2-Zn and 4-Zn energy differences to C−H bond activation (Figure 4). There is now a very large variation between the different DFT-based results with the hybrid functionals PBE0 and B3LYP and the meta-GGA M11-L giving the smallest differences and with M06-L, TPSS, and PBE giving larger differences. The DFT LSDA and PBE energy differences obtained from FHI-aims

RESULTS AND DISCUSSION Figure 3 shows the calculated energy difference in kcal/mol between R and TS (red bars) and TS and P (blue bars) for

Figure 3. Energy differences in kcal/mol for ethylene hydrogenation with E(TS)-E(R) in red and E(TS)-E(P) in blue; the error bars on the DMC values are ±1 standard deviation.

ethylene hydrogenation. The QE results are labeled QE(LSDA) and QE(PBE) for the LSDA and PBE xc functionals, respectively, and DMC(LSDA) and DMC(PBE) show the DMC values for LSDA and PBE initial trial wavefunctions. The result from ref 11 is labeled MN M06-L; the rest are results from FHI-aims as described in the Methods section. The figure clearly shows that there is not much difference between the various DFT xc functionals for R → TSthe differences between the various DFT functionals are smaller than the differences between DFT and DMC(LSDA). For TS → P, the

Figure 4. Energy differences in kcal/mol for Zn functionalized Zr6 nodes with E(TS)−E(R) in red for the 2-Zn node and blue for the 4Zn node; the error bars on the DMC values are ±1 standard deviation. E


Article


correct DMC energies of these states. Furthermore, the optimized spin structures clearly depend on which xc functional was used to optimize the structures. Exhaustively examining optimized spin structures for different xc functionals and constructing a multireference trial wavefunctions are very much beyond the scope of the present work. Finally, Figure 5 shows the results for the Ni-functionalized nodes. Again, there is a large variability in the energy

and QE are in good agreement with one another. The DMC(LSDA) and DMC(PBE) energy differences are within two standard deviations of one another. Again, the DMC(LSDA) energies for the R and TS states are lower than the DMC(PBE) ones, and so, we consider the DMC(LSDA) energy differences more accurate than the DMC(PBE) ones (see Table 2). The main striking difference between the DFTTable 2. DMC Energies in kcal/mol for the R and TS States for C−H Bond Activation Using 2-Zn, 4-Zn, 2-Ni, and 4-Ni Nodes Based on LSDA and PBE Nodal Surfaces ±1 Standard Deviationa state 2-Zn 4-Zn TS2-Zn TS4-Zn 2-Ni 4-Ni TS2-Ni TS4-Ni

LSDA 7.7 0 0 4.1 0 18.7 0 5.1

± ± ± ± ± ± ± ±

1.3 1.8 1.8 1.8 2.0 2.5 3.2 2.5

PBE 10.0 5.3 8.96 2.2 7.0 18.6 3.8 11.2

± ± ± ± ± ± ± ±

reference energy 2.1 1.1 1.7 2.1 2.0 1.5 1.7 3.6

4-Zn DMC(LSDA) TS2-Zn DMC(LSDA) DMC(LSDA) 2-Ni DMC(LSDA) TS2-Ni

a

The reference energies are the DMC(LSDA) energies for the 4-Zn, TS2-Zn, 2-Ni, and TS2-Ni states, respectively. The CH4 energies were −5097.45 ± 0.06 and −5097.54 ± 0.07 kcal/mol for DMC(LSDA) and DMC(PBE), respectively.

based differences and the DMC(LSDA) ones is that the DFT differences are all positive and now overestimate the energy differences relative to DMC. The PBE0 energy difference is smallest at 2.5 and 8.5 kcal/mol for the 2-Zn and 4-Zn nodes, respectively, with the 4-Zn energy difference in good agreement with the DMC(LSDA) value. In contrast with the 4-Zn energy difference, the DMC(LSDA) 2-Zn energy difference is negative or approximately zero to within two standard deviations [both for the DMC(LSDA) and the DMC(PBE) energy differences]. A negative energy difference for C−H bond activation implies that the reaction would be exothermic, which seems unphysical. We therefore performed some more calculations to investigate this further. We first re-optimized the 2-Zn node for a singlet structure using FHI-aims with the PBE xc functional and the default “tight” settings keeping the carbon atoms in the ligand fixed and with the maximum force on any (movable) atom less than 0.006 eV/Å. The singlet FHI-aims DFT PBE energy on this singlet-optimized structure was about 11 kcal/mol lower than the triplet DFT(PBE) energy on the triplet-optimized structure, and the QE DFT LSDA singlet energy was about 7 kcal/mol lower than the QE DFT LSDA triplet energy on the triplet-optimized structure. However, the DMC(LSDA) singlet energy on the singlet structure was −43 690.80(7) eV, compared to −43 690.79(3) for the triplet state on the triplet structurethese states were degenerate within DMC(LSDA). We also calculated the QE DFT LSDA and the DMC(LSDA) singlet energies on the triplet-optimized TS2-Zn structurethe DFT LSDA triplet energy was lower than the singlet one by about 1 kcal/mol, but the DMC(LSDA) energy of the singlet state on this structure was lower than the triplet state energy by about 16 kcal/mol, a considerable amount. These results suggest that the 2-Zn and TS2-Zn systems have degenerate or near-degenerate spin multiplets, with a consequence that a multireference trial wavefunction may be necessary to accurately capture the

Figure 5. Energy differences in kcal/mol for Ni functionalized Zr6 nodes with E(TS)−E(R) in red for the 2-Ni node and blue for the 4Ni node; the error bars on the DMC values are ±1 standard deviation.

differences from the DFT calculations with PBE0 equal to the DMC(LSDA) energy differences both for 2-Ni and 4-Ni to within the error bars of DMC(LSDA). Also, the FHI-aims LSDA and PBE energy differences for 2-Ni and 4-Ni are in good agreement with the corresponding QE(LSDA) and QE(PBE) values.

■

SUMMARY AND CONCLUSIONS We have here compared catalytic electronic energy differences obtained from DFT and DMC calculations for transitionmetal-functionalized heterogeneous Zr6 nodes linked in the NU-1000 MOF. There is a large variability, up to ∼10 kcal/ mol, in the DFT values for the energy differences. From our calculations, we cannot discern any particular trend in the DFT valuessometimes they overestimate the energy differences compared to the best DMC value and sometimes they underestimate the differences. The overall best DFT results are obtained using the hybrid PBE0 and HSE06 xc functionals. This is particularly true for the bimetallic-functionalized (Cu− Ni and Cu−Zn) Zr6 nodes, where the PBE0 values are in excellent agreement with the DMC(LSDA) values. In contrast, for the ethylene hydrogenation on a single-Ni-functionalized Zr6 node, even the PBE0 energy differences are grossly underestimated relative to the DMC(LSDA) values. For the HF-DFT calculations we did, the results are no better, in fact worse, than those obtained using standard DFT schema. This seems to indicate that the errors in the DFT calculations are not primarily driven by self-interaction errors, and the systems F


The Journal of Physical Chemistry C are “normal” in the context of density-driven errors.49 Indeed, this is consistent with the HOMO−LUMO gaps in these systems typically being a few eV. Instead, we suspect that the errors relative to DMC(LSDA) are in a large part due to van der Waals interactions, which are either not included in the DFT calculations here or are included approximately based on fitting to reference structures, but are included essentially exactly in the QMC calculations, and can be important in these large metallorganic clusters. We did run a test on the ethylene hydrogenation using PBE0 with the Tkatchenko and Scheffler van der Waals correction added in postprocessing (“vdw_correction_hirshfeld” in FHI-aims),66 but that resulted in a change of energy differences between the R and TS, and TS and P states of no more than 0.7 kcal/mol. It has been shown in previous studies that van der Waals corrections to DFT functional can be very sensitive to geometries and charge density distribution,67,68 even when materials are known to be pure van der Waals systems. The long-range van der Waals corrections to DFT functionals are usually fitted to reproduce a specific class of systems. Therefore, moving away from the class system for which the correction was generated reduces the confidence in or applicability of the corrections. In a study by Shin et al.67 on the energetics of graphene and graphyne bilayers, van der Waals corrected DFT functionals performed differently in graphene and graphyne when compared to DMC, despite being both carbon allotropes. This suggests that the role of geometry is important when correcting for dispersion forces, which could certainly be an issue in the systems studied here compared to the ones used when fitting the van der Waals corrections. The main concern about the DMC values is the quality of the QMC nodal surface. Because of the sheer size of these systems, it was not feasible to systematically optimize the QMC nodal surface, and we instead had to settle for comparing energies obtained from DFT-LSDA and DFTPBE nodal surfaces. The differences between DMC(LSDA) and DMC(PBE) energies were small, often within one standard deviation of one another but sometimes much larger than that. The DMC(LSDA) energies were always smaller than or equal to the DMC(PBE) energies (to within statistical error), which is why we have confidence that the energies, and therefore also the energy barriers, obtained from DMC(LSDA) are better than those from DMC(PBE). Finally, the uncovered degenerate or near-degenerate spin multiplet for the 2-Zn node within DMC suggests that it may sometimes be necessary to optimize different spin structures with a few different DFT xc functionals, such as PBE, metaGGA, and hybrid functionals to clearly map out the energy landscape of the spin multiplets; when different spin states are close in energy, the optimized structures may be sensitive to which xc functional is used to optimize the structure.

■

Article

■

ACKNOWLEDGMENTS

■

REFERENCES

We thank M. Ortuño J. Ye, A. B. League, C. J. Cramer, and L. Gagliardi for helpful discussions and comments on the manuscript and V. Blum for helpful discussions about FHIaims. A.B., H.S., Y.L., and O.H. were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program and by the ASCR Leadership Computing Challenge. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC0206CH11357. We gratefully acknowledge the computing resources provided on Blues and Bebop, high-performance computing clusters operated by the Laboratory Computing Resource Center at Argonne National Laboratory. D.P. gratefully acknowledges the financial support from the Inorganometallic Catalyst Design Center, an EFRC funded by the DOE, Office of Basic Energy Sciences (DESC0012702).

(1) Yaghi, O. M.; O’Keeffe, M.; Ockwig, N. W.; Chae, H. K.; Eddaoudi, M.; Kim, J. Reticular synthesis and the design of new materials. Nature 2003, 423, 705−714. (2) Horike, S.; Shimomura, S.; Kitagawa, S. Soft porous crystals. Nat. Chem. 2009, 1, 695−704. (3) Eddaoudi, M.; Sava, D. F.; Eubank, J. F.; Adil, K.; Guillerm, V. Zeolite-like metal−organic frameworks (ZMOFs): design, synthesis, and properties. Chem. Soc. Rev. 2015, 44, 228−249. (4) Férey, G. Hybrid porous solids: past, present, future. Chem. Soc. Rev. 2008, 37, 191−214. (5) Islamoglu, T.; Goswami, S.; Li, Z.; Howarth, A. J.; Farha, O. K.; Hupp, J. T. Postsynthetic tuning of metal−organic frameworks for targeted applications. Acc. Chem. Res. 2017, 50, 805−813. (6) Peters, A. W.; Li, Z.; Farha, O. K. Enhancing the catalytic activity in the solid state: Metal-organic frameworks to the rescue. ACS Cent. Sci. 2017, 3, 367−368. (7) Rogge, S. M. J.; Bavykina, A.; Hajek, J.; Garcia, H.; OlivosSuarez, A. I.; Sepúlveda-Escribano, A.; Vimont, A.; Clet, G.; Bazin, P.; Kapteijn, F.; et al. Metal−organic and covalent organic frameworks as single-site catalysts. Chem. Soc. Rev. 2017, 46, 3134−3184. (8) Lee, J.; Farha, O. K.; Roberts, J.; Scheidt, K. A.; Nguyen, S. T.; Hupp, J. T. Metal−organic framework materials as catalysts. Chem. Soc. Rev. 2009, 38, 1450−1459. (9) Ma, L.; Abney, C.; Lin, W. Enantioselective catalysis with homochiral metal−organic frameworks. Chem. Soc. Rev. 2009, 38, 1248−1256. (10) Vermoortele, F.; Vandichel, M.; Van de Voorde, B.; Ameloot, R.; Waroquier, M.; Van Speybroeck, V.; De Vos, D. E. Electronic effects of linker substitution on Lewis acid catalysis with metal− organic frameworks. Angew. Chem., Int. Ed. 2012, 51, 4887−4890. (11) Li, Z.; Schweitzer, N. M.; League, A. B.; Bernales, V.; Peters, A. W.; Getsoian, A.; Wang, T. C.; Miller, J. T.; Vjunov, A.; Fulton, J. L.; et al. Sintering-resistant single-site nickel catalyst supported by metal− organic framework. J. Am. Chem. Soc. 2016, 138, 1977−1982. (12) Pahls, D. R.; Ortuño, M. A.; Winegar, P. H.; Cramer, C. J.; Gagliardi, L. Computational screening of bimetal-functionalized Zr6O8 MOF nodes for methane C−H bond activation. Inorg. Chem. 2017, 56, 8739−8743. (13) Kim, I. S.; Borycz, J.; Platero-Prats, A. E.; Tussupbayev, S.; Wang, T. C.; Farha, O. K.; Hupp, J. T.; Gagliardi, L.; Chapman, K. W.; Cramer, C. J.; et al. Targeted single-site MOF node modification:

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Anouar Benali: 0000-0002-2133-0338 Olle Heinonen: 0000-0002-3618-6092 Notes

The authors declare no competing financial interest. G


Article

The Journal of Physical Chemistry C Trivalent metal loading via atomic layer deposition. Chem. Mater. 2015, 27, 4772−4778. (14) Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. 1964, 136, B864. (15) Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133. (16) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (17) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta− generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (18) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785. (19) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (20) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207−8215. (21) Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, G. Quantum Monte Carlo simulations of solids. Rev. Mod. Phys. 2001, 73, 33−83. (22) Song, S.; Kim, M.-C.; Sim, E.; Benali, A.; Heinonen, O.; Burke, K. Benchmarks and reliable DFT results for spin gaps of small ligand Fe (II) complexes. J. Chem. Theory Comput. 2018, 14, 2304−2311. (23) Wagner, L.; Mitas, L. A quantum Monte Carlo study of electron correlation in transition metal oxygen molecules. Chem. Phys. Lett. 2003, 370, 412−417. (24) Koseki, J.; Maezono, R.; Tachikawa, M.; Towler, M. D.; Needs, R. J. Quantum Monte Carlo study of porphyrin transition metal complexes. J. Chem. Phys. 2008, 129, 085103. (25) Horváthová, L.; Dubecký, M.; Mitas, L.; Š tich, I. Spin multiplicity and symmetry breaking in vanadium-benzene complexes. Phys. Rev. Lett. 2012, 109, 053001. (26) Benali, A.; Shulenburger, L.; Romero, N. A.; Kim, J.; von Lilienfeld, O. A. Application of diffusion Monte Carlo to materials dominated by van der Waals interactions. J. Chem. Theory Comput. 2014, 10, 3417−3422. (27) Mostaani, E.; Szyniszewski, M.; Price, C. H.; Maezono, R.; Danovich, M.; Hunt, R. J.; Drummond, N. D.; Fal’ko, V. I. Diffusion quantum Monte Carlo study of excitonic complexes in twodimensional transition-metal dichalcogenides. 2017, arXiv:1706.04688 [cond-mat]. (28) Caffarel, M.; Applencourt, T.; Giner, E.; Scemama, A. Communication: Toward an improved control of the fixed-node error in quantum Monte Carlo: The case of the water molecule. J. Chem. Phys. 2016, 144, 151103. (29) Shin, H.; Kang, S.; Koo, J.; Lee, H.; Kim, J.; Kwon, Y. Cohesion energetics of carbon allotropes: Quantum Monte Carlo study. J. Chem. Phys. 2014, 140, 114702. (30) Luo, Y.; Benali, A.; Shulenburger, L.; Krogel, J. T.; Heinonen, O.; Kent, P. R. C. Phase stability of TiO2 polymorphs from diffusion Quantum Monte Carlo. New J. Phys. 2016, 18, 113049. (31) Benali, A.; Shulenburger, L.; Krogel, J. T.; Zhong, X.; Kent, P. R. C.; Heinonen, O. Quantum Monte Carlo analysis of a charge ordered insulating antiferromagnet: the Ti4O7 Magnéli phase. Phys. Chem. Chem. Phys. 2016, 18, 18323−18335. (32) Wagner, L. K.; Abbamonte, P. Effect of electron correlation on the electronic structure and spin-lattice coupling of high-Tc cuprates: Quantum Monte Carlo calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 125129. (33) Yu, J.; Wagner, L. K.; Ertekin, E. Fixed-node diffusion Monte Carlo description of nitrogen defects in zinc oxide. Phys. Rev. B 2017, 95, 075209. (34) Yu, J.; Wagner, L. K.; Ertekin, E. Towards a systematic assessment of errors in diffusion Monte Carlo calculations of

semiconductors: Case study of zinc selenide and zinc oxide. J. Chem. Phys. 2015, 143, 224707. (35) Drummond, N. D.; Needs, R. J.; Sorouri, A.; Foulkes, W. M. C. Finite-size errors in continuum quantum Monte Carlo calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 125106. (36) Santana, J. A.; Krogel, J. T.; Kent, P. R. C.; Reboredo, F. A. Cohesive energy and structural parameters of binary oxides of groups IIA and IIIB from diffusion quantum Monte Carlo. J. Chem. Phys. 2016, 144, 174707. (37) Santana, J. A.; Krogel, J. T.; Kim, J.; Kent, P. R. C.; Reboredo, F. A. Structural stability and defect energetics of ZnO from diffusion quantum Monte Carlo. J. Chem. Phys. 2015, 142, 164705. (38) Grimm, R. C.; Storer, R. G. Monte-Carlo solution of Schrödinger’s equation. J. Comput. Phys. 1971, 7, 134−156. (39) Kalos, M. H.; Levesque, D.; Verlet, L. Helium at zero temperature with hard-sphere and other forces. Phys. Rev. A: At., Mol., Opt. Phys. 1974, 9, 2178−2195. (40) Morales, M. A.; McMinis, J.; Clark, B. K.; Kim, J.; Scuseria, G. E. Multideterminant wave functions in quantum Monte Carlo. J. Chem. Theory Comput. 2012, 8, 2181−2188. (41) Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-bath configuration interaction: An efficient selected configuration interaction algorithm inspired by heat-bath sampling. J. Chem. Theory Comput. 2016, 12, 3674−3680. (42) Zheng, H.; Wagner, L. K. Computation of the correlated metalinsulator transition in vanadium dioxide from first principles. Phys. Rev. Lett. 2015, 114, 176401. (43) Shin, H.; Luo, Y.; Ganesh, P.; Balachandran, J.; Krogel, J. T.; Kent, P. R. C.; Benali, A.; Heinonen, O. Electronic properties of doped and defective NiO: A quantum Monte Carlo study. Phys. Rev. Mater. 2017, 1, 073603. (44) Zhou, X.; Wang, F. Barrier heights of hydrogen-transfer reactions with diffusion quantum monte carlo method. J. Comput. Chem. 2017, 38, 798−806. (45) Krongchon, K.; Busemeyer, B.; Wagner, L. K. Accurate barrier heights using diffusion Monte Carlo. J. Chem. Phys. 2017, 146, 124129. (46) Swann, E. T.; Coote, M. L.; Barnard, A. S.; Per, M. C. Efficient protocol for quantum Monte Carlo calculations of hydrogen abstraction barriers: Application to methanol. Int. J. Quantum Chem. 2017, 117, No. e25361. (47) Peverati, R.; Truhlar, D. G. Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics. Philos. Trans. R. Soc., A 2014, 372, 20120476. (48) Mondloch, J. E.; Bury, W.; Fairen-Jimenez, D.; Kwon, S.; DeMarco, E. J.; Weston, M. H.; Sarjeant, A. A.; Nguyen, S. T.; Stair, P. C.; Snurr, R. Q.; et al. Vapor-phase metalation by atomic layer deposition in a metal−organic framework. J. Am. Chem. Soc. 2013, 135, 10294−10297. (49) Kim, M.-C.; Park, H.; Son, S.; Sim, E.; Burke, K. Improved DFT potential energy surfaces via improved densities. J. Phys. Chem. Lett. 2015, 6, 3802−3807. (50) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. Comput. Phys. Commun. 2009, 180, 2175. (51) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 13244. (52) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (53) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision E.01; Gaussian Inc.; Wallingford CT, 2009. H


Article

The Journal of Physical Chemistry C (54) Peverati, R.; Truhlar, D. G. M11-L: A local density functional that provides improved accuracy for electronic structure calculations in chemistry and physics. J. Phys. Chem. Lett. 2011, 3, 117−124. (55) Kim, J.; Esler, K. P.; Mcminis, J.; Morales, M. A.; Clark, B. K.; Shulenburger, L.; Ceperley, D. M. Hybrid algorithms in quantum Monte Carlo. J. Phys.: Conf. Ser. 2012, 402, 012008. (56) Mathuriya, A.; Luo, Y.; Benali, A.; Shulenburger, L.; Kim, J. Optimization and parallelization of B-spline based orbital evaluations in QMC on multi/many-core shared memory processors. 2017 IEEE International Parallel and Distributed Processing Symposium (IPDPS), 2017; pp 213−223. (57) Mathuriya, A.; Luo, Y.; Clay, R. C., III; Benali, A.; Shulenburger, L.; Kim, J. Embracing a new era of highly efficient and productive quantum Monte Carlo simulations. Proceedings Of The International Conference for High Performance Computing, Networking, Storage And Analysis, New York, NY, USA, 2017; pp 38:1−38:12. (58) Kim, J.; Baczewski, A. D.; Beaudet, T. D.; Benali, A.; Bennett, M. C.; Berrill, M. A.; Blunt, N. S.; Borda, E. J. L.; Casula, M.; Ceperley, D. M.; et al. QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids. J. Phys.: Condens. Matter 2018, 30, 195901. (59) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, C. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (60) OPIUM Package v. 3.7 2014. http://opium.sourceforge.net. (61) Mitra, C.; Krogel, J. T.; Santana, J. A.; Reboredo, F. A. Manybody ab initio diffusion quantum Monte Carlo applied to the strongly correlated oxide NiO. J. Chem. Phys. 2015, 143, 164710. (62) Foyevtsova, K.; Krogel, J. T.; Kim, J.; Kent, P. R. C.; Dagotto, E.; Reboredo, F. A. Ab initio quantum Monte Carlo calculations of spin superexchange in cuprates: The benchmarking case of Ca2CuO3. Phys. Rev. X 2014, 4, 031003. (63) Shin, H.; Benali, A.; Luo, Y.; Crabb, E.; Lopez-Bezanilla, A.; Ratcliff, L. E.; Jokisaari, A. M.; Heinonen, O. Zirconia and hafnia polymorphsground state structural properties from diffusion Monte Carlo. 2017, arXiv:1708.09424 [cond-mat.mtrl-sci]. (64) Dubecký, M. Noncovalent interactions by fixed-node diffusion Monte Carlo: Convergence of nodes and energy differences vs Gaussian basis-set size. J. Chem. Theory Comput. 2017, 13, 3626− 3635. (65) Luo, Y.; Esler, K. P.; Kent, P. R.; Shulenburger, L. An efficient hybrid orbital representation for quantum Monte Carlo calculations. 2018, arXiv:1805.07406. (66) Tkatchenko, A.; Scheffler, M. Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. Phys. Rev. Lett. 2009, 102, 073005. (67) Shin, H.; Kim, J.; Lee, H.; Heinonen, O.; Benali, A.; Kwon, Y. Nature of interlayer binding and stacking of hybridized carbon layers: A quantum Monte Carlo Study. J. Chem. Theory Comput. 2017, 13, 5639−5646. (68) Shulenburger, L.; Baczewski, A. D.; Zhu, Z.; Guan, J.; Tománek, D. The nature of the interlayer interaction in bulk and few-layer phosphorus. Nano Lett. 2015, 15, 8170−8175.

I


Quantum Monte Carlo Calculations of Catalytic Energy Barriers in a

Recommend Documents