Density-matrix based determination of low-energy model Hamiltonians from ab initio wavefunctions
We propose a way of obtaining effective low energy Hubbard-like model Hamiltonians from ab initio quantum Monte Carlo (QMC) calculations for molecular and extended systems. The Hamiltonian parameters are fit to best match the ab initio two-body density matrices and energies of the ground and excited states, and thus we refer to the method as ab initio density matrix based downfolding (AIDMD). For benzene (a finite system), we find good agreement with experimentally available energy gaps without using any experimental inputs. For graphene, a two dimensional solid (extended system) with periodic boundary conditions, we find the effective on-site Hubbard to be , comparable to a recent estimate based on the constrained random phase approximation. For molecules, such parameterizations enable calculation of excited states that are usually not accessible within ground state approaches. For solids, the effective Hamiltonian enables large-scale calculations using techniques designed for lattice models.
The reliable simulation of systems for which the large-scale physics is not well-approximated by a non-interacting model, is a major challenge in physics, chemistry, and materials science. These systems appear to require a multi-scale approach in which the effective interactions between electrons at a small distance scale are determined, which then leads to a coarse-grained description of emergent correlated physics. This reduction of the Hilbert space is often known as ”downfolding”. In strongly-correlated systems, the correct effective Hamiltonian is strongly dependent on material-specific properties, motivating the need for a generic accurate method to determine it.
One can loosely categorize downfolding techniques into two strategies. The first strategy is based on performing ab initio calculations and then matching them state by state to the effective model. Alternately, some approaches employ a model for the screening of Coulomb interactions, for which the ab initio single particle wavefunctions provide the relevant inputs. For the purposes of this manuscript, the umbrella term we use for these strategies is ”fitting”. Techniques that fall into this class include the constrained density functional theory Pavirini (); Dasgupta (), the constrained random phase approximation (cRPA) Aryasetiawan (), fitting spin models to energies Valenti_kagome (); Spaldin (), and efforts by one of us Wagner_JCP () using reduced density matrices of quantum Monte Carlo (QMC) calculations. The second class is based on Löwdin downfolding Freed (); Zhou_Ceperley (); Tenno () and canonical transformation theory Glazek_Wilson (); Wegner (); White_CT (); Yanai_CT (); Watson_Chan (), which involves a sequence of unitary transformations on the bare Hamiltonian, chosen in a way that minimize the size of the matrix elements connecting the relevant low energy (valence) space to the high energy one.
Downfolding by fitting has the advantage that it is conceptually straightforward to perform, although it demands an a priori parameterization of the effective Hamiltonian. The methods have been applied to complex bulk systems Imada1 (); Imada2 (); Arya1 (); Arya2 (); Scriven (); Wehling_graphene (), but it is only recently that their accuracy is being rigorously checked RPA_Troyer (). On the other hand, canonical transformations do not need such parameterizations and can discover the relevant terms in an automated way. However, their application to complex materials remains to be carried out and tested.
Here we propose a scheme which aims to capture the ”best of both worlds”. On the one hand, we retain the simplicity of fitting and on the other we use information from accurate many-body wavefunctions to determine which terms are important. The deviations between the ab initio and model properties allows us to assess the quality of the resultant model and to discover relevant physics from the calculation. Simultaneously, the method cannot depend too much on the quality of the ab initio solution because of the inherent limitations of accuracy, especially for very big system sizes.
Once an effective model Hamiltonian in the reduced Hilbert space is obtained, as is depicted in Fig. 1, it can be used to perform a calculation on a larger system using techniques designed for lattice models dmrg_white (); tps_nishino (); Vidal_MERA (); TPS_review (); Changlani_CPS (); Neuscamman_CPS (); mezzacapo (); Marti (); DMET_2012 (); Chen_DMET (); Sandvik_loops (); Blankenbecler (); Alavi_FCIQMC (); SQMC (). This multi-step modeling procedure is needed since the ab initio calculations for a given system size are, in general, computationally more expensive than the equivalent lattice calculations. Large sizes are crucial to study finite size effects, and in turn theoretically establish the presence of a phase. In addition, excited states and dynamical correlation functions have traditionally been difficult in ab initio approaches, but have seen progress for lattice model methods dynamical_dmrg (); Daley_tDMRG (); White_tDMRG ().
In this paper, we demonstrate a downfolding method that uses data from ab initio QMC techniques to derive an effective coarse-grained Hamiltonian. This method, which we call non-eigenstate ab initio density matrix downfolding (N-AIDMD), uses many-body simulations of non-eigenstates to fit an effective low-energy Hamiltonian. We demonstrate that the method can use wavefunctions of medium quality to derive highly accurate effective Hamiltonians. After demonstrating a simple example, we downfold benzene from a 30-electron problem to a 6-electron one and show that the resulting Hamiltonian reproduces the experimental spectrum well. We also show that the method works for extended systems, by applying it to graphene.
In the present section we discuss the methods we used to generate our ab initio data. While most of our discussion is specific to QMC, the quantities used can also be calculated in almost any other wavefunction-based quantum chemistry method.
ii.1 Variational and Fixed Node Diffusion Quantum Monte Carlo
Ab initio QMC comprises of a suite of methods that efficiently sample the phase space of electrons each moving in 3-dimensional real space. When the wavefunction as a function of coordinates is known, the phase space can be sampled with variational Monte Carlo (VMC) using Metropolis algorithms. For ground state calculations, the Diffusion Monte Carlo (DMC) method, based on imaginary time evolution of the Schroedinger equation, is formally exact but in practice limited by the numerical sign problem. This problem ceases to exist if one knew the exact location of the nodes (zeroes) of the many-body wavefunction. Thus, the optimal strategy for very accurate calculations is to generate a good trial wavefunction, and optimize its parameters to minimize its variational energy. Then, we use this wavefunction to ”fix the nodes” (which may only approximately correspond to the exact nodes) and perform a DMC calculation under this constraint. This last variant is called the fixed-node DMC (FN-DMC) method and is known to be very accurate for a large variety of systems. While some more details are discussed here, we refer the interested reader to Ref. Foulkes_review, for an exhaustive review of concepts and applications. All the ab initio QMC calculations were carried out with the QWalk package QWALK ().
A typical QMC calculation was carried out as follows. First, we perform DFT calculations with the B3LYP B3LYP () or PBE functionals PBE () using GAMESS GAMESS () for molecules or CRYSTAL CRYSTAL () for solids. The lowest energy DFT orbitals provide the Slater determinant part of the trial wavefunction. For molecules, a multi-determinantal wavefunction is generated by performing a configuration interaction with singles and doubles excitations (CISD) calculation from the reference Slater determinant. Once this is done, a Jastrow factor is introduced, resulting in the ansatz for the trial wavefunction ,
where refers to the coordinates of the electrons (the spin indices, being fixed, have been suppressed), are determinants and their corresponding coefficients. When we desire eigenstates, the parameters in the Jastrow and the coefficients are optimized to get the best possible variational energy within the ansatz chosen using a technique introduced by Umrigar and coworkers Umrigar_optimization () with an efficient algorithm by Clark et al. Clark_multidet ().
The variational energy is calculated via Metropolis sampling of ,
where is a compact notation for the coordinates of the electrons and is the ”local energy”. With this trial wavefunction we perform FN-DMC to calculate the energy using the mixed (or projected) estimator,
where is obtained by a stochastic projection of under the constraint that and have the same sign everywhere.
We now discuss measurements in the QMC methods. For a generic operator , the pure (VMC) and mixed estimators are computed as,
The mixed estimator of an operator is equal to the pure estimator in two cases; (1) when is the exact wavefunction or (2) when the operator commutes with the Hamiltonian. In more general situations, higher accuracy can be obtained with the extrapolated estimatorCeperley_Kalos_Binder () for approximate eigenstates,
For accurate wavefunctions, all these estimators must approach the same value.
We will construct effective Hamiltonians using the two-body reduced density matrix (2-RDM) elements, given by the estimator (the normalization has been omitted),
where refers to the set of electron coordinates obtained by changing the location of two electrons and are a chosen set of one-particle wavefunctions (orbitals) indexed by label . The mixed estimator equivalent of Eq. (6) is obtained like that for the energy. More details of this computation have been previously discussed elsewhere by one of us Wagner_JCP (). The chosen set of orbitals is often localized on the atoms; this property makes it convenient to derive Hubbard-like models. We explain their construction next.
ii.2 Localized orbitals
Localized orbitals often provide an intuitive way of understanding an electronic system in terms of electron hops and on-site or inter-site repulsions. Thus, many works have been devoted to this subject; ranging from the Linearized Muffin-Tin Orbital (LMTO) method Andersen () to the maximally localized Wannier function construction Marzari_Wf (). Orbital localization has also been widely discussed in the quantum chemistry literature.
The idea is to first select a set of orbitals in a certain energy window. For solids these correspond to bands close to the Fermi level, and for molecules these are partially occupied orbitals which constitute the active space. Then, a unitary transformation is performed to optimize a pre-decided metric for localization. In this work, we minimize the spread ,
where are the desired localized orbitals related to the chosen set of orbitals by a unitary transformation, .
For some systems, as we will see in the case of benzene and graphene, schematically shown in Fig. 2, it is necessary to include unoccupied orbitals to get well-localized orbitals of the right symmetry Hansen (). Thus the construction of localized orbitals is not a black-box procedure and may need adaptations based on the specifics of the system.
We note that the optimized parameters of the effective Hamiltonian may, in general, depend on the choice of localized orbitals. However, we have not explored this dependence - our main objective in this paper is to assess the validity of model Hamiltonians with respect to ab initio calculations for a particular choice of single particle basis.
ii.3 Lattice model calculations
The lattice model calculations for Hubbard models of benzene and graphene at half-filling were carried out with a combination of our own codes and the freely available QUEST determinantal quantum Monte Carlo package QUEST (). For the honeycomb lattice half-filled Hubbard model, the determinantal QMC method is sign problem free and the results are exact up to statistical errors. A time step of was chosen and (the imaginary time) was set to for every calculation. Measurements were performed for sweeps, with an additional sweeps being used for equilibration.
Iii Criteria for a low energy effective Hamiltonian
Our aim is to obtain a low energy effective Hamiltonian defined in the active space of electrons. In this basis, the criteria it must satisfy are:
The reduced density matrices (RDM) of the ground and excited states obtained from the ab initio calculation must match with that of the model calculation.
The energy spectra of the ab initio and model systems must match in the energy window of interest.
The model must be detailed enough to capture the essential physics and yet simple enough to avoid over-parameterization and over-fitting.
The concept of matching RDMs [criterion (a)] has previously appeared in related contexts Acioli (); Zhou_Ceperley (); Changlani_percolation () and in work by one of us Wagner_JCP (). Most physical properties, such as the charge and spin structure factors, are functions of the 2-RDM. Since it is computationally expensive to calculate high-order RDMs, we use the matching condition only on the 2-RDM, where are orbital indices (including space and spin). This criterion automatically ensures that the combined number of electrons occupying the orbitals is equal to those in the model Hamiltonian. If any input state does not have the expected electron number in the active space, it can not be described by the effective Hamiltonian.
The importance of excited state energies used in the fitting [criterion (b)] is highlighted by the fact that the wavefunctions, and their corresponding two-body density matrices, are invariant to many kinds of terms that enter the Hamiltonian. For example, the transformation,
is, by construction, consistent with all the 2-RDM data for any , , for systems which have definite spin symmetry and particle number. Imposing certain physical constraints on the form of the interactions can reduce the need for this criterion. To give a concrete example, consider wavefunction data generated from the ground state of an unfrustrated Heisenberg spin Hamiltonian on a bipartite lattice Note1 (), where refer to nearest neighbor pairs and is the spin operator on site . Then adding gives the same reduced density matrices for the ground state, as long as is small enough to not cause energy crossings i.e. not make an original excited state the new ground state. This additional term has the effect of introducing long-range Heisenberg couplings. Moreover, the effective Hamiltonian is not unique; the Lieb-Mattis model LM () (where () and refer to the sublattice and corresponding spin), is also known to reproduce the low-energy limit of the Heisenberg model. Thus, imposing the requirement that the Hamiltonian has the nearest-neighbor form constrains to zero and picks one particular model. Similar arguments should apply to extended Hubbard models in homogeneous systems where a physical constraint is that the density-density interaction must decrease monotonically with distance between orbitals.
Iv Ab Initio Density matrix based Downfolding (AIDMD) procedures
We now discuss two procedures that use density matrices and energies to calculate parameters entering the effective Hamiltonian; both have been schematically depicted in Fig. 3.
iv.1 Eigenstate AIDMD method
In the first method, schematically depicted in Fig. 3(a), eigenstates from an ab initio calculation are used to match density matrices and energies of the corresponding model. The QMC extrapolated estimator is taken to be an accurate representation of the true one. Then the parameters of the model Hamiltonian are obtained by minimizing a cost function that is a linear combination of the energy and density matrix errors,
where the subscript is an eigenstate index, are orbital indices and the superscripts and refer to ab initio and model calculations respectively. There is no definite prescription for choosing the weight ; a good heuristic is to choose a value that gives roughly the same size of errors for the two terms that enter the cost function. The cost minimization is performed with the Nelder Mead simplex algorithm.
In practice we found that since the number of available eigenstates and the accuracy of true estimators is limited, a second method discussed next is more suited for downfolding.
iv.2 Non-eigenstate AIDMD method
Consider a set of ab initio energy averages , i.e. expectation values of the Hamiltonian, and corresponding 1- and 2-RDMs , for arbitrary low-energy states characterized by index . Assume a model 2-body Hamiltonian with effective parameters (1-body part) and (2-body part) along with a constant term ; the total number of parameters being . Then for each state , we have the equation,
where we have made the assumption that the chosen set of operators corresponding to single particle wavefunctions or orbitals, explain all energy differences seen in the ab initio data. The constant is from energetic contributions of all other orbitals which are not part of the chosen set.
We then perform calculations for low-energy states which are not necessarily eigenstates. These states are not arbitrary in the sense that they have similar descriptions of the core and virtual spaces. Each state satisfies the criteria (1) its energy average does not lie outside the energy window of interest and (2) the trace of its 1-RDM matches the electron number expected in the effective Hamiltonian.
The objective of choosing a sufficiently big set of states is to explore parts of the low-energy Hilbert space that show variations in the RDM elements. Since the same parameters describe all states, they must satisfy the linear set of equations,
which is compactly written as,
where is the dimensional vector of energies, is the matrix composed of density matrix elements and is a dimensional vector of parameters. This problem is over-determined for , which is the regime we expect to work in.
In the case of any imperfection in the model, which is the most common case, the equality will not hold exactly and one must then instead minimize the norm of the error, :
This cost function can be minimized in a single step by using the method of least squares, employing the singular value decomposition of matrix . This matrix also encodes exact (or near-exact) linear dependences. Thus, the quality of the fit is directly judged by assessing (1) the singular values of the matrix and (2) the value of the cost function itself i.e. the deviations of the input and fitted energies. We will refer to this as the non eigenstate (N)-AIDMD method throughout the rest of the paper. This idea is schematically depicted in Fig. 3(b).
The matrix gives a very natural basis for understanding renormalization effects. For example, consider a set of wavefunctions that show that the correlator does not change significantly. This would lead to the corresponding column of matrix being identical (up to a scale factor) to the first column of 1’s. Physically, this would correspond to the coupling constant being irrelevant for the low-energy physics; it can take any value including zero and can be absorbed into the constant shift term. This could also alternatively mean that the input data is correlated and does not provide enough information about , so care must be taken in constructing the set of wavefunctions.
In summary, the N-AIDMD method performs the following operation. The 1- and 2-RDMs and energy expectation values of many non-eigenstate correlated states are calculated. Then, given an effective Hamiltonian parameterization, linear equations (12) are constructed and solved. Standard model fitting principles apply, and we can evaluate the goodness of fit to determine whether a given effective Hamiltonian can sufficiently describe the data.
iv.3 Generating states for AIDMD methods
We now address the central issue of generating states to be used as inputs for the AIDMD methods.
For the E-AIDMD, the near-eigenstates were obtained by performing CISD calculations with multiple roots and optimizing a multi-determinant Jastrow wavefunction with each CISD guess as a starting point. This is known to be approximate, especially for higher excited states. It is the inherent uncertainty about the accuracy of this procedure, along with the fact that only a small number of eigenstates are accessible, that limits the utility of E-AIDMD.
From the point of view of the N-AIDMD method too, automating the construction of the database of wavefunctions may not be completely straightforward and here we offer some heuristics for doing this within QMC methods. We re-emphasize that any state described by the effective Hamiltonian must be one that does not involve large contributions from the core and virtual orbitals i.e. single particle degrees of freedom outside the chosen active space. This check can be imposed at the ab initio level by monitoring the RDMs, for example, the trace of the 1-RDM taken over the active orbitals must equal the number of electrons in the effective Hamiltonian description.
One way to generate new states is to perturb near-eigenstates. For example, after optimizing the multi-determinantal-Jastrow trial wavefunction, we artificially change the determinantal coefficients by small amounts. This procedure changes the nodal surface and gives energies close to, but different from, the optimized ground state. A second source of data is spin excitations of the DFT reference Slater determinant, generated within the space of orbitals that play an important role in the active space; for benzene and graphene these involve the Kohn-Sham orbitals with symmetry. Finally, in the case of extended systems, we chose a linear combination of determinants which, in spite of being not size-extensive, reveal additional properties of the effective Hamiltonian.
iv.4 Quantum Monte Carlo specific adaptation
The formalism introduced above applies to any method that calculates energies and density matrices. In this paper, all the expectation values entering the matrix are calculated for the chosen low-energy wavefunctions by Monte Carlo sampling.
Once this database of states has been generated, we perform two independent calculations to estimate the parameters, one using variational and the other using mixed estimators. In the latter case, we modify the linear equations (12) by using and the projected estimates of the density matrix elements i.e. and in the construction of the matrix. The implicit normalization of these mixed estimates by is assumed. This projector formulation is also amenable to coupled-cluster calculations which work with projected energies and density matrices.
The bias arising out of fixing the nodes of the projected wavefunction does not affect our formulation. This is because the method regards as some arbitrary positive sampling function associated with a low energy state, and it is this same distribution that is used for the evaluation of the density matrix elements. Each such distribution provides a linear equation encoding the relationship between the FN-DMC energy and the projected density matrix elements and the unknown parameters. Up to errors from statistical uncertainties and from the assumption of the form of the Hamiltonian, this relationship is exact.
However, the value of the optimal parameters does depend on the choice of method. For example, for the benzene molecule presented in section VI, the VMC and FN-DMC parameter values agree with each other up to ; the largest discrepancy is due to different constant terms. This discrepancy is expected, because VMC does not provide an accurate description of the core and virtual spaces.
Ideally, only the FN-DMC calculations should be used to estimate the parameters. However, the mixed estimator in FN-DMC is biased because of the inaccuracy of the trial wavefunction. Thus, we propose a better estimator for the parameters,
where is the true parameter, and and are the corresponding parameters obtained from FN-DMC and VMC calculations. The details of this result are explained in Appendix A.
V Simple application: Hubbard to Heisenberg model
To demonstrate our formalism for a simple example, we consider the two site Hubbard model and fit information from the lowest two states to a Heisenberg model.
We analytically solve for all four eigenstates of the Hamiltonian,
for two opposite spin electrons on two sites, where is the hopping, is the Hubbard on-site interaction. The Hilbert space on a single site (orbital) is spanned by four states (unoccupied), (single up occupied), (single down occupied), (doubly occupied). For completeness, we discuss some features of the solution method below.
First notice that the triplet state with energy and the state with energy , are exact eigenstates of the problem independent of the values of and .
To get the other two states, write the Hamiltonian in the basis of and ,
Then diagonalizing it, we get the energies to be,
The lowest energy corresponds to the singlet, and the corresponding eigenvector is
with the next excited state being the triplet .
We choose the Heisenberg form to fit to
To determine the parameters and , form the matrix with the lowest two energy states,
Using derived values of and , we get,
which to lowest order in is .
Observe that the correlator for is not exactly 1/4 but only approximately so. This is expected since the fluctuations from the high-energy states are not exactly zero, if it were, it would be equivalent to exactly block-diagonalizing the Hamiltonian. This exact block diagonalization is not possible in general, unless it is also accompanied with a change in the low energy degrees of freedom entering the model.
If we now rotate the two low energy eigenstates and define the orthogonal linear combinations,
and calculate their energies and construct the corresponding matrix, we get to be independent of and . This property is desirable as the method does not hinge on the requirement of eigenstates as inputs.
In terms of canonical transformations (here equivalent to second order perturbation theory), the matrix element of the effective Hamiltonian between the singly occupied states is,
Since this matrix element must equal in the Heisenberg model we arrive at the same result, namely Note2 ().
Vi Application to a molecule: Benzene
We show the workings of the described AIDMD methods for the benzene molecule. Our choice of system is motivated by its simplicity (one band model and presence of many symmetries) and the availability of experimental energies to compare to.
DFT calculations were first performed in the TZVP basis with Burkatzki-Filippi-Dolg pseudopotentials BFD (), using the molecular geometry that corresponds to a B3LYP optimized calculation. These serve as a starting point for the QMC calculations, discussed in the Methods section. In the charge neutral sector, there are a total of 30 electrons (for example, 15 and 15 for the spin-singlet state) and our objective is to downfold this system to an effective one with 6 electrons (3, 3 ).
The model Hamiltonian is defined in the space of six orbitals; a representative localized orbital has been shown in Fig. 2. These orbitals were obtained by localizing the highest three occupied and the lowest three unoccupied B3LYP orbitals (from the DFT calculation) with orbital symmetry, a well established procedure in the literature Hansen (). The overall phase of these orbitals is adjusted to enable use of parameter symmetries directly when fitting.
Appendix B discusses the QMC data used for fitting and the initial pre-processing to determine the eligibility of states that can be described by a six- orbital Hamiltonian.
vi.1 On-site Hubbard model
We consider the Hubbard model for six orbitals of benzene, given by Eq. (15), where is used for the nearest-neighbor hopping and is the effective on-site Coulomb repulsion. We will discuss multiple ways of using reduced density matrix elements to estimate these parameters.
vi.1.1 Hubbard from the E-AIDMD method
An estimate of is obtained by directly matching the half-filled ground state () 2-RDM element corresponding to the ”double occupancy” correlator () of the ab initio and lattice-model calculations. This element equals for the non interacting case () and its value reduces for .
Fig. 4 shows the dependence of , computed in the ground state of the Hubbard model of a six-site ring at half filling, on . The plot also indicates the value of this correlator computed from various wavefunctions and estimators from ab initio QMC calculations of the benzene molecule. The trends are consistent with our expectations; the Slater-Jastrow (SJ) wavefunction at the VMC level underestimates the strength of the effective interactions, which is partly remedied by the extrapolated estimator from FN-DMC. However, the bias (systematic error) is expected to be large because of the considerable difference in the two estimates. This bias is reduced by the multi-determinantal-Jastrow (CISDJ) wavefunction we employed; the difference between the variational and extrapolated estimator is about 5%.
The value of is found to be extremely sensitive to the precise value of the double occupancy correlator, a change of a few percent (i.e. from to ) changes our estimate from to (i.e. a factor of almost ). In general, this observation suggests that it is crucial to look at various other elements of the 2-RDM and to look at alternate ways of estimating Hubbard parameters.
In Fig. 5 we compare results of other correlators from extrapolated QMC estimates with those from the on-site Hubbard model on a six site ring for varying values of . We focus on the one-body density matrix (the values are the same for both spin indices ), density-density correlators and spin-spin correlators , all as a function of distance () with respect to a reference site (labelled ). The value of gives the smallest errors for most observables, except for the nearest-neighbor density-density correlator which favors a value of . In the limit that the model is perfect, all estimates must yield the same value; the differences reflect an inadequacy of the on-site Hubbard model in describing all the data. This is evidence for the need for long-range interactions.
vi.1.2 Hubbard from the N-AIDMD method
As mentioned previously, the idea of matching density matrix elements is useful only for comparing exact eigenstates. However, it is difficult to construct eigenstates with very high accuracy in the ab initio calculations and at times also for the equivalent model for large system sizes. This is why we appeal to the N-AIDMD method, introduced and explained in Sec. IV.2, which is relatively insensitive to the nature of the states input to the method. For charge-neutral benzene, we construct the matrix by taking various states in a eV energy window above the ground state using VMC and DMC methods.
Fig. 6 shows the comparison of the fitted energy and the input VMC or DMC energy. The former is obtained by taking the linear combination of the ab initio VMC or DMC density matrices weighted by the optimized parameters of the effective Hamiltonian. A perfect agreement between the fitted and ab initio input data would correspond to all energies falling exactly on the line. By this measure, the Hubbard model for benzene is reasonable, though not accurate, as is seen in Fig. 6(a) and (d). The presence of significant deviations of the order of eV from the line indicates the need for a more refined model, which we discuss in section VI.2.
The VMC data yields optimal parameters of eV and eV () and the FN-DMC data gives eV and eV (). The extrapolated estimate of the optimal parameters, eV agrees with the value of eV reported by Bursill et al. Bursill (). The extrapolated value of is also broadly consistent with a recently reported estimate Schuler_graphene () to within 10-20%.
|Parameter||PPP [Ref. Bursill, ]||PPP-VMC||PPP-DMC||PPP-Extrap||[Ref. Schuler_graphene, ]||-VMC||-DMC|
vi.2 Extended Hubbard model
Having established the need for long-range interactions in benzene, we consider the extended-Hubbard or Parisier-Pople-Parr (PPP) model,
where is the on-site Hubbard interaction and and are inter-orbital hopping and density-density interactions. We compare our results with Bursill et al Bursill (), who considered only a nearest neighbor hopping and took to be of the Ohno form Ohno ()
where is a fit parameter and is the spatial separation between nuclei. This parameterization has been widely used in the modelling of various organic polymers. Here do not make any assumptions about the form of the interactions and instead use the N-AIDMD method to determine these parameter values.
We repeat analyses similar to those for the Hubbard model, in addition to carefully looking at the variations in the 1- and 2-RDM matrix elements. This data has been discussed as part of Appendix B and has been shown in Fig. 11. For a parameter to be reliably estimated there should be a large variation in the corresponding density matrix element for different wavefunctions in the low energy space. By this metric, we find that the next nearest neighbor hopping is irrelevant in the charge-neutral sector. We thus attempt to fit to a model only with the nearest neighbor along with , and ; is not needed as it simply sets the chemical potential.
The inadequacies of the on-site Hubbard model, shown in Figs. 6(a) and (d), are rectified by the extended one, shown in Figs. 6(b) and (e); the maximum energy errors of about eV are reduced to eV. The root mean square errors are much smaller as well, reducing from eV to about eV. Adding the next-next nearest neighbor hopping , shown in Figs. 6(c) and (f), only marginally improves the accuracy of the fit; is found to be only about a tenth of the value of suggesting that its effects can be largely accounted for by .
To assess the accuracy of these parameters, we compare the results of the model with experimentally available energies. First, as Fig. 7(a) shows for the extended-Hubbard or PPP model, the extrapolated parameters give an improved agreement with experiment compared to the VMC or FN-DMC parameters. Most energy gaps of this model are in excellent agreement with available experimental energies, the errors are eV or less. The largest outlier at about eV is within of the experimental result.
Next, we compare the experimental energies with energy gaps from various model Hamiltonians. While the Hubbard and the Ohno parameterizations reproduce most experimental gaps, especially at low energies, they have at least one significant outlier. These outliers are correctly accounted for by the fitted PPP Hamiltonian, and improved upon by the introduction of . Owing to the small value of , more data in the N-AIDMD method may be needed to precisely estimate its value. The extreme sensitivity of the high energy eigenstates to may explain the deviation of the model gap from the experimental one at about eV.
We now discuss our parameter values and the errors associated with them; these have been summarized in Table I. While our PPP parameters are generally consistent with the Ohno form, there are some differences of the order of eV, that improve the quality of the fitted energies. We emphasize that we have not provided any experimental inputs; rather we have used only the QMC data (energies and density matrices) from multiple states to obtain the Hamiltonian parameters.
In order to check the robustness of the fit, we estimated errors in our parameters from a Jackknife analysis. In this scheme half the input states to the N-AIDMD method were randomly discarded and the fit performed with the retained half. Many such randomly generated ensembles of input data were taken and the resultant parameters were averaged over all of these. The difference of the parameters of this reduced data set and those obtained from the full data set provides an estimate of the systematic error, which we report in table I. The other source of error is from statistical noise, which in the present case was found to be much smaller than the systematic error.
Vii Application to a solid: Graphene
As an application of the AIDMD methods to solid materials, we consider graphene, a 2D solid of carbon atoms arranged on a honeycomb lattice. Graphene has great potential technological applications, which has spurred much work devoted to understanding it thoroughly graphene_review (). In addition, there have been several proposals for engineering exotic phases in graphene Herbut (); Duric (); Ghaemi (). That said, it is only recently that systematic studies to estimate the role of electron-electron interactions Wehling_graphene (); Schuler_graphene (); Abbamonte () have been carried out. While some of its long-distance properties appear to be adequately described by a tight-binding model, the short range features, crucial for phenomena such as magnetism, require more refined modeling.
Early studies modelled graphene as a honeycomb lattice Hubbard model with a estimated to be . This would put graphene on the verge of a metal-insulator transition Sorella_Tosatti (); Sorella_Nature (). However, recent results by Wehling et al. realized the importance of long range interactions Wehling_graphene (), which renormalize the on-site interaction to an effectively lower value. Schuler et al. Schuler_graphene () report the effective to be , which means graphene lies well in the semimetal phase of the honeycomb lattice Hubbard model.
Setting aside the question of determining all the long range interactions in graphene, we ask what best describes our ground state QMC data. To do so, we first generated the -like Wannier functions within QWalk QWALK (), a representative of which has been shown in Fig. 2.
Just as in the case of benzene, we used the fact that the effective strength of the Coulomb interaction is most sensitive to the 2-RDM element . For the unit cell with periodic boundary conditions, and using optimized Slater Jastrow wavefunctions, the extrapolated value is found to be corresponding to a . This estimate of is obtained from comparisons to lattice determinantal QMC calculations, which were carried out for sizes ranging from to to check for finite-size effects. As Fig. 8 shows, the unit cell is distinctly different from the larger unit cells and the finite-size errors in the double occupation correlator are negligible beyond sizes . The finite size effects for other short range correlation functions (not shown) are also negligible beyond .
We also calculated many non-eigenstates for graphene in an energy window eV above the ground state. Fig. 9 shows the VMC and FN-DMC energy fits to tight-binding and Hubbard models. The hopping for the tight-binding model is found to be in the range of to eV, as is indicated in Figs. 9(a) and (c). However, a precise estimate of this parameter is not particularly meaningful because the model is inadequate at capturing many states, particularly spinful excitations, in this energy window. We note that a value of eV Schuler_graphene () has been previously calculated with DFT methods.
Figs. 9(b) and (d) show that the Hubbard model reduces the errors of the tight-binding model, and the value of is found to be in the range of (VMC) to (FN-DMC). The latter estimate is expected to be more accurate and hence closer to the true value of in a small energy window associated with the ground state. We also note that this value is within of the value derived from the constrained RPA parameters Schuler_graphene (). However the Hubbard model too has outliers of about eV, which are large for an accurate model of a solid material.
This inadequacy is confirmed by assessing various correlators in the half filled ground state. Fig. 10 shows the up density-up density and up density- down density correlators as a function of distance between carbon atoms. On the scale of Fig. 10(a) (and well within the accuracy of our calculations) the like spin correlations were captured well for all values of in the range from to . However, as Fig. 10(b) shows, the Hubbard model for large tends to exaggerate the the effective interaction between the two electron spin flavors at small separations. In particular, the nearest neighbor unlike spin density-density correlator is found to be in better agreement with than any finite value. This, just like the case of benzene, suggests the need for longer range interactions in the model. There are also small deviations between the ab initio QMC and the Hubbard model results at longer distances. These correlations do not depend significantly on and these data in isolation do not rule out any model.
We have demonstrated ab initio density matrix downfolding (AIDMD) methods where ab initio quantum Monte Carlo (QMC) data is used to fit simple effective Hamiltonians. We have elaborated on the fitting procedures and the intricacies of the QMC method needed to perform calculations. The limitations of the model were judged by assessing the quality of the fitted energies and 2-body density matrices. This feature is useful for constructing refined models needed for the accurate simulation of real materials.
For the benzene molecule, while the on-site Hubbard model with was able to capture most features of the QMC ground state data, the deviations of the density matrices revealed the need for longer range interactions. Including these interactions improved the agreement of the model with both the QMC results and the experimental data. This effective Hamiltonian parameterization could be used to calculate low-frequency response functions and to check semi-empirical methods.
Since QMC calculations use size-consistent wavefunctions for extended systems and scale favorably, we believe the type of calculations presented here will be a promising alternative to DFT-based downfolding approaches for solid materials. Our demonstration for the single band model of graphene yielded an effective , in the same range as a recently reported estimate based on the constrained-RPA method Schuler_graphene (). We leave a more detailed characterization of interactions in graphene to future work.
In more complicated materials, where the form of the Hamiltonian is unclear, we suggest that the dominant terms can be obtained from canonical transformation theory followed with an accurate fit to the QMC data. In this spirit, it will also be useful to compare the predictions of the proposed AIDMD schemes with other complementary proposals for downfolding Zhang (); Zgid (); Watson_Chan ().
Finally, we remark that previously unsolved model Hamiltonians are now being accurately treated with tensor network methods Stoudenmire_White_2D_DMRG (); Corboz_Rice_Troyer (). With parallel advances in the ab initio QMC simulation of high-temperature superconductors Wagner_cuprates (); Krogel_cuprates (), a clear future direction is to deduce more refined models for these compounds, using the ideas discussed in the paper.
We thank David Ceperley, Cyrus Umrigar, Garnet Chan, Shiwei Zhang, Steven White, Gabriel Kotliar, Bryan Clark, Norman Tubman, Victor Chua and Miles Stoudenmire for discussions. We would also like to thank Cyrus Umrigar for running very useful checks against his quantum Monte Carlo code (CHAMP) in the early stages of this work. We acknowledge support from SciDAC grant DE-FG02-12ER46875. Computational resources were provided by the DOE INCITE SuperMatSim program and the Taub campus cluster at the University of Illinois at Urbana-Champaign.
X Appendix A: Estimation of true effective Hamiltonian parameters from VMC and FN-DMC calculations
In section IV.4, we discussed the discrepancy between parameters obtained from VMC and FN-DMC. This is attributed to the inability of the trial wavefunctions used in VMC to provide an accurate description of the core and virtual spaces.
As mentioned in the text, ideally, only the FN-DMC calculations should be used to estimate the parameters. However, the mixed estimator in FN-DMC is biased because of the inaccuracy of the trial wavefunction. To minimize this bias, we propose estimators which combine FN-DMC and VMC parameters.
Consider a trial VMC wavefunction which deviates from the DMC wavefunction by a small amount orthogonal to i.e.,
To obtain the parameter coupled to an operator in the effective Hamiltonian i.e.
we take partial derivatives with respect to the change of a density matrix element i.e.
Since the pure estimators for projected (DMC) wavefunctions are not easily evaluated in QMC, we use other estimators using which we indirectly obtain .
The parameter obtained from the mixed estimators within FN-DMC, is formally defined as,
On substituting the relation between and wavefunctions, we get,
which to linear order in is,
A similar expression for parameters obtained from VMC calculations,
is derived and leads to the result,
where we have used the hermiticity of the Hamiltonian and the operator .
Xi Appendix B: QMC data for the benzene molecule
Table II shows the QMC data for several eigenstates of benzene. Most of these states along with many other non-eigenstates constitute our data set for fitting a model.
|Spin||DFT||SJ-VMC||SJ-DMC||CISDJ-VMC||CISDJ-DMC||Used in Fit?|
The calculations confirm the general expectation that significant energy gains are obtained by improving wavefunctions going from the single-Slater-Jastrow form to the multi-determinantal-Jastrow form (the determinants being selected from a CISD calculation). Moreover, the DMC calculations improve total energies significantly; typically, the DMC values are eV lower than the corresponding VMC value.
The total electron count from the one-body density matrix is assessed to verify the validity of fitting to a six- orbital Hamiltonian. Table II shows that for charge-neutral benzene, the singlet state () has up and down electron counts of 2.96 each, which are close to the expected values of 3,3. For the state, roughly eV above the ground state, the deviations were slightly larger; the summed occupation numbers were 4.92 and 1.0 in comparison to the expected values of 5 and 1. Since there is a slight deviation of these numbers from integers, we rescale the one and two body density matrices by factors (all slightly greater than 1) that correctly accounts for sum rules for each individual state used in the AIDMD methods.
However, for the state, the electron occupation numbers deviate significantly from the corresponding value in the model; almost by one integer. This indicates that the state is inadequately described by the proposed Hamiltonian and hence can not be used in the fitting procedure. This deviation is not completely unexpected since this state is eV above the ground state and a potential high-energy state. Said differently, the active space at this energy scale is considerably different from that assumed for the ground state and its low-energy excitations. Thus, this QMC data suggests that it is only reasonable for the effective Hamiltonian concept to hold only in an energy window of the order of eV above the ground state.
We now assess some aspects of our non-eigenstate data. As mentioned in the main text, heuristics were used to construct these states. For example, for every near-eigenstate in a symmetry sector that was represented by multi-determinantal Jastrow form, we changed the determinantal coefficients to generate new wavefunctions. We checked the 1-RDM to make sure that it had the correct total electron number (6 electrons) on the localized orbitals that constitute our active space. Moreover, if the change in determinantal coefficients led to an energy-average outside our pre-decided energy window, the new state was discarded from the N-AIDMD procedure.
For the N-AIDMD method, we desire large variations in the density matrix elements for different wavefunctions in the low-energy space, in order to accurately estimate the Hamiltonian parameters. Fig. 11 shows these variations for all relevant density-matrix elements (within FN-DMC) that were needed for estimating the parameters of the extended Hubbard model. was found to be irrelevant, as the summed density matrix elements coupling to it were found to not vary significantly (not shown in the plot).
- (1) E. Pavarini et al., Phys. Rev. Lett. 87, 047003 (2001).
- (2) O. K. Andersen and T. Saha-Dasgupta, Phys. Rev. B 62, R16219 (2000).
- (3) F. Aryasetiawan et al., Phys. Rev. B 70, 195104 (2004).
- (4) H. O. Jeschke, F. Salvat-Pujol, and R. Valentí, Phys. Rev. B 88, 075106 (2013).
- (5) N. S. Fedorova, C. Ederer, N. A. Spaldin, A. Scaramucci, arXiv:1412.3702.
- (6) L. K. Wagner, The Journal of Chemical Physics 138, 094106 (2013).
- (7) K. F. Freed, Accounts of Chemical Research 16, 137 (1983).
- (8) S. Q. Zhou and D. M. Ceperley, Phys. Rev. A 81, 013402 (2010).
- (9) S. Ten-no, The Journal of Chemical Physics 138, (2013).
- (10) S. D. Głazek and K. G. Wilson, Phys. Rev. D 48, 5863 (1993).
- (11) F. Wegner, Annalen der Physik 506, 77 (1994).
- (12) S. R. White, The Journal of Chemical Physics 117, (2002).
- (13) T. Yanai and G. K.-L. Chan, The Journal of Chemical Physics 124, (2006).
- (14) T. J. Watson Jr., G. K-L Chan, arXiv:1502.04698.
- (15) M. Hirayama, T. Miyake, and M. Imada, Phys. Rev. B 87, 195144 (2013).
- (16) K. Nakamura, Y. Yoshimoto, Y. Nohara, and M. Imada, Journal of the Physical Society of Japan 79, 123708 (2010).
- (17) F. Nilsson, R. Sakuma, and F. Aryasetiawan, Phys. Rev. B 88, 125123 (2013).
- (18) F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schönberger, Phys. Rev. B 74, 125106 (2006).
- (19) E. P. Scriven and B. J. Powell, Phys. Rev. Lett. 109, 097206 (2012).
- (20) T. O. Wehling et al., Phys. Rev. Lett. 106, 236805 (2011).
- (21) Hiroshi Shinaoka, Rei Sakuma, Philipp Werner, Matthias Troyer, arXiv:1410.1276.
- (22) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- (23) A. Gendiar, N. Maeshima, and T. Nishino, Prog. Theor. Phys. 110, 691 (2003).
- (24) G. Vidal, Phys. Rev. Lett. 101, 110501 (2008).
- (25) F. Verstraete, V. Murg, and J. Cirac, Advances in Physics 57, 143 (2008).
- (26) H. J. Changlani, J. M. Kinder, C. J. Umrigar, and G. K.-L. Chan, Phys. Rev. B 80, 245116 (2009).
- (27) E. Neuscamman, H. Changlani, J. Kinder, and G. K.-L. Chan, Phys. Rev. B 84, 205132 (2011).
- (28) F. Mezzacapo, N. Schuch, M. Boninsegni, and J. I. Cirac, New J. Phys. 11, 083026 (2009).
- (29) K. H. Marti et al., New Journal of Physics 12, 103008 (2010).
- (30) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
- (31) Q. Chen et al., Phys. Rev. B 89, 165134 (2014).
- (32) O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
- (33) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- (34) G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 131, 054106 (2009).
- (35) F. R. Petruzielo et al., Phys. Rev. Lett. 109, 230201 (2012).
- (36) E. Jeckelmann, Phys. Rev. B 66, 045114 (2002).
- (37) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
- (38) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- (39) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- (40) L. K. Wagner, M. Bajdich, and L. Mitas, Journal of Computational Physics 228, 3390 (2009).
- (41) A.D. Becke, J.Chem.Phys. 98 (1993) 5648-5652; C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785-789; S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200-1211; P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J.Phys.Chem. 98 (1994) 11623-11627.
- (42) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (43) M.W.Schmidt, K.K.Baldridge, J.A.Boatz, S.T.Elbert, M.S.Gordon, J.H.Jensen, S.Koseki, N.Matsunaga, K.A.Nguyen, S.Su, T.L.Windus, M.Dupuis, J.A.Montgomery J. Comput. Chem., 14, 1347-1363(1993).
- (44) R.Dovesi, R.Orlando, A.Erba, C.M. Zicovich-Wilson, B.Civalleri, S.Casassa, L.Maschio, M.Ferrabone, M. De La Pierre, P. D’ Arco, Y. Noel, M. Causa, M. Rerat, B. Kirtman. Int. J. Quantum Chem. 114, 1287 (2014).
- (45) J. Toulouse and C. J. Umrigar, The Journal of Chemical Physics 126, (2007).
- (46) B. K. Clark et al., The Journal of Chemical Physics 135, (2011).
- (47) Quantum Many-Body Problems in Monte Carlo Methods in Statistical Physics, ed. K. Binder, Springer-Verlag, 1979.
- (48) O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
- (49) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
- (50) K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. Lett. 94, 026405 (2005).
- (51) QUEST: Quantum Electron Simulation Toolbox, http://quest.ucdavis.edu/index.html.
- (52) P. H. Acioli and D. M. Ceperley, The Journal of Chemical Physics 100, (1994).
- (53) H. J. Changlani, S. Ghosh, S. Pujari, and C. L. Henley, Phys. Rev. Lett. 111, 157201 (2013).
- (54) A bipartite lattice is one with two sublattices and with only connections but no or ones.
- (55) E. Lieb and D. Mattis, Journal of Mathematical Physics 3, (1962).
- (56) To map a fermion model to a spin model, a Jordan-Wigner phase is introduced which remains fixed, as there are no charge fluctuations. In this notation the relative signs between basis states entering in the singlet and triplet states is reversed and .
- (57) M. Burkatzki, C. Filippi, and M. Dolg, The Journal of Chemical Physics 126, (2007).
- (58) R. J. Bursill, C. Castleton, and W. Barford, Chemical Physics Letters 294, 305 (1998).
- (59) M. Schüler et al., Phys. Rev. Lett. 111, 036601 (2013).
- (60) K. Ohno, Theoretica chimica acta 2, 219 (1964).
- (61) A. H. Castro Neto et al., Rev. Mod. Phys. 81, 109 (2009).
- (62) I. F. Herbut, Phys. Rev. Lett. 97, 146401 (2006).
- (63) T. Đurić, N. Chancellor, and I. F. Herbut, Phys. Rev. B 89, 165123 (2014).
- (64) P. Ghaemi, J. Cayssol, D. N. Sheng, and A. Vishwanath, Phys. Rev. Lett. 108, 266801 (2012).
- (65) J. P. Reed et al., Science 330, 805 (2010).
- (66) S. Sorella and E. Tosatti, EPL (Europhysics Letters) 19, 699 (1992).
- (67) S. Sorella, Y. Otsuka and Seiji Yunoki, Scientific Reports 2, Article number:992 (2012).
- (68) W. Purwanto, S. Zhang, and H. Krakauer, Journal of Chemical Theory and Computation 9, 4825 (2013).
- (69) Alexander A. Rusakov, Jordan J. Phillips, Dominika Zgid, arXiv:1409.2921.
- (70) E. Stoudenmire and S. R. White, Annual Review of Condensed Matter Physics 3, 111 (2012).
- (71) P. Corboz, T. M. Rice, and M. Troyer, Phys. Rev. Lett. 113, 046402 (2014).
- (72) L. K. Wagner and P. Abbamonte, Phys. Rev. B 90, 125129 (2014).
- (73) K. Foyevtsova et al., Phys. Rev. X 4, 031003 (2014).