Densitymatrix based determination of lowenergy model Hamiltonians from ab initio wavefunctions
We propose a way of obtaining effective low energy Hubbardlike 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 twobody 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 onsite 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 largescale calculations using techniques designed for lattice models.
I Introduction
The reliable simulation of systems for which the largescale physics is not wellapproximated by a noninteracting model, is a major challenge in physics, chemistry, and materials science. These systems appear to require a multiscale approach in which the effective interactions between electrons at a small distance scale are determined, which then leads to a coarsegrained description of emergent correlated physics. This reduction of the Hilbert space is often known as ”downfolding”. In stronglycorrelated systems, the correct effective Hamiltonian is strongly dependent on materialspecific 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 manybody 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 multistep 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 coarsegrained Hamiltonian. This method, which we call noneigenstate ab initio density matrix downfolding (NAIDMD), uses manybody simulations of noneigenstates to fit an effective lowenergy 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 30electron problem to a 6electron 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.
Ii Methods
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 wavefunctionbased 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 3dimensional 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 manybody 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 fixednode DMC (FNDMC) 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 multideterminantal 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 ,
(1) 
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 ,
(2) 
where is a compact notation for the coordinates of the electrons and is the ”local energy”. With this trial wavefunction we perform FNDMC to calculate the energy using the mixed (or projected) estimator,
(3) 
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,
(4) 
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,
(5) 
For accurate wavefunctions, all these estimators must approach the same value.
We will construct effective Hamiltonians using the twobody reduced density matrix (2RDM) elements, given by the estimator (the normalization has been omitted),
(6) 
where refers to the set of electron coordinates obtained by changing the location of two electrons and are a chosen set of oneparticle 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 Hubbardlike 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 onsite or intersite repulsions. Thus, many works have been devoted to this subject; ranging from the Linearized MuffinTin 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 predecided metric for localization. In this work, we minimize the spread ,
(7) 
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 welllocalized orbitals of the right symmetry Hansen (). Thus the construction of localized orbitals is not a blackbox 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 halffilling 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 halffilled 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 overparameterization and overfitting.
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 2RDM. Since it is computationally expensive to calculate highorder RDMs, we use the matching condition only on the 2RDM, 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 twobody density matrices, are invariant to many kinds of terms that enter the Hamiltonian. For example, the transformation,
(8) 
is, by construction, consistent with all the 2RDM 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 longrange Heisenberg couplings. Moreover, the effective Hamiltonian is not unique; the LiebMattis model LM () (where () and refer to the sublattice and corresponding spin), is also known to reproduce the lowenergy limit of the Heisenberg model. Thus, imposing the requirement that the Hamiltonian has the nearestneighbor 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 densitydensity 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,
(9) 
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 Noneigenstate AIDMD method
Consider a set of ab initio energy averages , i.e. expectation values of the Hamiltonian, and corresponding 1 and 2RDMs , for arbitrary lowenergy states characterized by index . Assume a model 2body Hamiltonian with effective parameters (1body part) and (2body part) along with a constant term ; the total number of parameters being . Then for each state , we have the equation,
(10) 
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 lowenergy 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 1RDM 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 lowenergy Hilbert space that show variations in the RDM elements. Since the same parameters describe all states, they must satisfy the linear set of equations,
(11) 
which is compactly written as,
(12) 
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 overdetermined 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, :
(13) 
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 nearexact) 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 lowenergy 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 NAIDMD method performs the following operation. The 1 and 2RDMs and energy expectation values of many noneigenstate 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 EAIDMD, the neareigenstates were obtained by performing CISD calculations with multiple roots and optimizing a multideterminant 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 EAIDMD.
From the point of view of the NAIDMD 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 reemphasize 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 1RDM 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 neareigenstates. For example, after optimizing the multideterminantalJastrow 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 KohnSham orbitals with symmetry. Finally, in the case of extended systems, we chose a linear combination of determinants which, in spite of being not sizeextensive, 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 lowenergy 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 coupledcluster 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 FNDMC 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 FNDMC 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 FNDMC calculations should be used to estimate the parameters. However, the mixed estimator in FNDMC is biased because of the inaccuracy of the trial wavefunction. Thus, we propose a better estimator for the parameters,
(14) 
where is the true parameter, and and are the corresponding parameters obtained from FNDMC 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,
(15) 
for two opposite spin electrons on two sites, where is the hopping, is the Hubbard onsite 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 ,
(16) 
Then diagonalizing it, we get the energies to be,
(17) 
The lowest energy corresponds to the singlet, and the corresponding eigenvector is
(18) 
with the next excited state being the triplet .
We choose the Heisenberg form to fit to
(19) 
To determine the parameters and , form the matrix with the lowest two energy states,
(20) 
Using derived values of and , we get,
(21) 
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 highenergy states are not exactly zero, if it were, it would be equivalent to exactly blockdiagonalizing 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,
(22a)  
(22b) 
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,
(23)  
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 BurkatzkiFilippiDolg 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 spinsinglet 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 preprocessing to determine the eligibility of states that can be described by a six orbital Hamiltonian.
vi.1 Onsite Hubbard model
We consider the Hubbard model for six orbitals of benzene, given by Eq. (15), where is used for the nearestneighbor hopping and is the effective onsite Coulomb repulsion. We will discuss multiple ways of using reduced density matrix elements to estimate these parameters.
vi.1.1 Hubbard from the EAIDMD method
An estimate of is obtained by directly matching the halffilled ground state () 2RDM element corresponding to the ”double occupancy” correlator () of the ab initio and latticemodel 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 sixsite 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 SlaterJastrow (SJ) wavefunction at the VMC level underestimates the strength of the effective interactions, which is partly remedied by the extrapolated estimator from FNDMC. 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 multideterminantalJastrow (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 2RDM 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 onsite Hubbard model on a six site ring for varying values of . We focus on the onebody density matrix (the values are the same for both spin indices ), densitydensity correlators and spinspin 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 nearestneighbor densitydensity 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 onsite Hubbard model in describing all the data. This is evidence for the need for longrange interactions.
vi.1.2 Hubbard from the NAIDMD 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 NAIDMD method, introduced and explained in Sec. IV.2, which is relatively insensitive to the nature of the states input to the method. For chargeneutral 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 FNDMC 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 1020%.
Parameter  PPP [Ref. Bursill, ]  PPPVMC  PPPDMC  PPPExtrap  [Ref. Schuler_graphene, ]  VMC  DMC 

2.54  2.65(2)  2.54  
10.06  9.89(6)  3.04  
7.18  6.78(8)      
5.11  5.40(4)      
4.57  4.57  4.57  4.57      
3.96  4.16(2)  3.96(2)  3.73(3)  1.20  1.7(1)  1.39(5) 
vi.2 Extended Hubbard model
Having established the need for longrange interactions in benzene, we consider the extendedHubbard or ParisierPopleParr (PPP) model,
(24) 
where is the onsite Hubbard interaction and and are interorbital hopping and densitydensity 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 ()
(25) 
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 NAIDMD 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 2RDM 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 chargeneutral 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 onsite 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 nextnext 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 extendedHubbard or PPP model, the extrapolated parameters give an improved agreement with experiment compared to the VMC or FNDMC 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 NAIDMD 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 NAIDMD 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 electronelectron interactions Wehling_graphene (); Schuler_graphene (); Abbamonte () have been carried out. While some of its longdistance properties appear to be adequately described by a tightbinding 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 metalinsulator transition Sorella_Tosatti (); Sorella_Nature (). However, recent results by Wehling et al. realized the importance of long range interactions Wehling_graphene (), which renormalize the onsite 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 2RDM 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 finitesize effects. As Fig. 8 shows, the unit cell is distinctly different from the larger unit cells and the finitesize 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 noneigenstates for graphene in an energy window eV above the ground state. Fig. 9 shows the VMC and FNDMC energy fits to tightbinding and Hubbard models. The hopping for the tightbinding 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 tightbinding model, and the value of is found to be in the range of (VMC) to (FNDMC). 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 densityup 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 densitydensity 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.
Viii Conclusion
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 2body density matrices. This feature is useful for constructing refined models needed for the accurate simulation of real materials.
For the benzene molecule, while the onsite 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 lowfrequency response functions and to check semiempirical methods.
Since QMC calculations use sizeconsistent wavefunctions for extended systems and scale favorably, we believe the type of calculations presented here will be a promising alternative to DFTbased 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 constrainedRPA 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 hightemperature 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.
Ix Acknowledgement
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 DEFG0212ER46875. Computational resources were provided by the DOE INCITE SuperMatSim program and the Taub campus cluster at the University of Illinois at UrbanaChampaign.
X Appendix A: Estimation of true effective Hamiltonian parameters from VMC and FNDMC calculations
In section IV.4, we discussed the discrepancy between parameters obtained from VMC and FNDMC. 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 FNDMC calculations should be used to estimate the parameters. However, the mixed estimator in FNDMC is biased because of the inaccuracy of the trial wavefunction. To minimize this bias, we propose estimators which combine FNDMC and VMC parameters.
Consider a trial VMC wavefunction which deviates from the DMC wavefunction by a small amount orthogonal to i.e.,
(26) 
To obtain the parameter coupled to an operator in the effective Hamiltonian i.e.
(27) 
we take partial derivatives with respect to the change of a density matrix element i.e.
(28) 
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 FNDMC, is formally defined as,
(29) 
On substituting the relation between and wavefunctions, we get,
(30) 
which to linear order in is,
(31) 
A similar expression for parameters obtained from VMC calculations,
(32) 
is derived and leads to the result,
(33) 
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 noneigenstates constitute our data set for fitting a model.
Spin  DFT  SJVMC  SJDMC  CISDJVMC  CISDJDMC  Used in Fit?  

0  37.6303  37.6229(6)  37.7213(9)  37.6352(6)  37.7259(9)  2.96,2.96  Yes 
1  37.4634  37.4546(6)  37.5555(7)  37.4814(6)  37.5707(7)  3.94,1.98  Yes 
37.4561(6)  37.5479(6)  3.94,1.98  Yes  
37.4531(6)  37.5470(6)  3.94,1.98  Yes  
2  37.3203  37.2987(6)  37.3974(7)  37.3141(6)  37.4020(7)  4.92,1.00  Yes 
3  37.0378  37.0116(4)  37.1074(7)  37.0118(4)  37.1083(7)  4.88,0.02  No 
The calculations confirm the general expectation that significant energy gains are obtained by improving wavefunctions going from the singleSlaterJastrow form to the multideterminantalJastrow 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 onebody density matrix is assessed to verify the validity of fitting to a six orbital Hamiltonian. Table II shows that for chargeneutral 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 highenergy state. Said differently, the active space at this energy scale is considerably different from that assumed for the ground state and its lowenergy 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 noneigenstate data. As mentioned in the main text, heuristics were used to construct these states. For example, for every neareigenstate in a symmetry sector that was represented by multideterminantal Jastrow form, we changed the determinantal coefficients to generate new wavefunctions. We checked the 1RDM 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 energyaverage outside our predecided energy window, the new state was discarded from the NAIDMD procedure.
For the NAIDMD method, we desire large variations in the density matrix elements for different wavefunctions in the lowenergy space, in order to accurately estimate the Hamiltonian parameters. Fig. 11 shows these variations for all relevant densitymatrix elements (within FNDMC) 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).
References
 (1) E. Pavarini et al., Phys. Rev. Lett. 87, 047003 (2001).
 (2) O. K. Andersen and T. SahaDasgupta, Phys. Rev. B 62, R16219 (2000).
 (3) F. Aryasetiawan et al., Phys. Rev. B 70, 195104 (2004).
 (4) H. O. Jeschke, F. SalvatPujol, 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. Tenno, 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. KL 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) 56485652; C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785789; S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 12001211; P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J.Phys.Chem. 98 (1994) 1162311627.
 (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, 13471363(1993).
 (44) R.Dovesi, R.Orlando, A.Erba, C.M. ZicovichWilson, 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 ManyBody Problems in Monte Carlo Methods in Statistical Physics, ed. K. Binder, SpringerVerlag, 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 JordanWigner 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).