Timestep targeting timedependent and dynamical density matrix renormalization group algorithms with ab initio Hamiltonians
Abstract
We study the dynamical density matrix renormalization group (DDMRG) and timedependent density matrix renormalization group (tdDMRG) algorithms in the ab initio context, to compute dynamical correlation functions of correlated systems. We analyze the strengths and weaknesses of the two methods in small model problems, and propose two simple improved formulations, DDMRG and tdDMRG, that give increased accuracy at the same bond dimension, at a nominal increase in cost. We apply DDMRG to obtain the oxygen coreexcitation energy in the water molecule in a quadruplezeta quality basis, which allows us to estimate the remaining correlation error in existing coupled cluster results. Further, we use DDMRG to compute the local density of states and gaps, and tdDMRG to compute the complex polarization function, in linear hydrogen chains with up to 50 H atoms, to study metallicity and delocalization as a function of bondlength.
]Timestep targeting timedependent and dynamical density matrix renormalization group algorithms with ab initio Hamiltonians \abbreviations
1 Introduction
The calculation of dynamical quantities is essential for the interaction between theory and experiment. Most commonly, dynamical quantities such as the singleparticle Green’s function or optical absorption are considered in the linear response regime. In the frequency domain, the linear response of a wavefunction to a field can be written as the second derivative of a Lagrangian [McWeeny (1992), Helgaker et al. (2012)] and frequencydomain response theory in quantum chemistry has closely followed the theory of analytic energy derivatives, similar to that in structural optimization. Thus algorithms exist to compute dynamical correlation functions from HartreeFock [Dalgarno and Victor (1966)], density functional theory [Casida (1995)], configuration interaction [Koch and Harrison (1991)], coupled cluster [Koch and Jo/rgensen (1990)], and JastrowSlater wavefunctions [Mussard et al. (2017), Zhao and Neuscamman (2016)] amongst others, using analytic derivative techniques. Dynamical quantities can also be calculated in the timedomain. Here, quantum chemical methods typically formulate the equation of motion for the wavefunction from the DiracFrenkel (timedependent) variational principle [Krause et al. (2007), Li et al. (2005), Hoodbhoy and Negele (1978)]. Both kinds of algorithms can be found implemented in many modern quantum chemistry codes.
Dynamical quantities have also been studied with density matrix renormalization group (DMRG) or matrix product state (MPS) wavefunctions. Here a wide range of numerical algorithms have been explored. In the frequency domain, the first dynamical correlation functions were computed in a fixed linear space of DMRG renormalized states (i.e. by optimizing a single tensor in the MPS) [Hallberg (1995)]. Subsequent algorithms, such as the dynamical DMRG (DDMRG) [Jeckelmann (2002), Ramasesha et al. (1997), Pati et al. (1999), Kühner and White (1999)] or analytic DMRG response theory [Dorando et al. (2009), Nakatani et al. (2014)], further considered the response of the DMRG renormalized basis (i.e. all tensors in the MPS). DDMRG is widely used as a benchmark method for DMRG dynamical correlation functions, but unlike the analytical DMRG response theory does not correspond to a true derivative of a Lagrangian. The analytic DMRG response theory is equivalent to the later “tangent space” formulations of DMRG dynamical correlation functions [Haegeman et al. (2013)].
Timepropagation has also been investigated in conjunction with DMRG wavefunctions. Although a wide variety of timepropagation algorithms have been discussed [Cazalilla and Marston (2002), Vidal (2003), Daley et al. (2004), Feiguin and White (2005), White and Feiguin (2004), Kinder et al. (2011), Haegeman et al. (2011), Zaletel et al. (2015)], some, such as timeevolving block decimation [Vidal (2003)], are specialized to Hamiltonians with shortrange interactions on a 1D lattice. For quantum chemistry, it is necessary to work with longrange interactions, and one of the early timedependent DMRG (tdDMRG) algorithms that supported such Hamiltonians was the timestep targeting timedependent DMRG method [Feiguin and White (2005)]. There have also been many other important developments in timedependent DMRG which we do not discuss here, including translating timepropagation algorithms such as Chebyshev expansion and Krylov space techniques to work with MPS [Holzner et al. (2011), Ganahl et al. (2014), Wolf et al. (2015)], analytic timepropagation using the timedependent variational principle [Kinder et al. (2011), Haegeman et al. (2011)], and matrix product operator representations of the timeevolution operator with improved global timestep error [Zaletel et al. (2015)].
In the current work, we explore frequencydependent and timedependent DMRG algorithms for dynamical quantities to better understand the behaviour and applicability of these algorithms in the ab initio DMRG context [White (1992), White (1993), White and Martin (1999), Daul et al. (2000), Mitrushenkov et al. (2001), Chan and HeadGordon (2002), Legeza et al. (2003), Chan (2004), Schollwöck (2005), Moritz and Reiher (2006), Marti et al. (2008), Marti and Reiher (2010), Schollwöck (2011), Sharma and Chan (2012), Wouters and Van Neck (2014), OlivaresAmaya et al. (2015)]. There has been relatively little work computing ab initio dynamical quantities with DMRG. Earlier work in our group compared dynamical DMRG and analytic DMRG response theory for computing frequency dependent polarizabilities [Dorando et al. (2009)]. Subsequent investigations exploited the analogy between the analytic DMRG response theory and the random phase approximation to obtain DMRG excitation energies and RPAlike correlation energy contributions for some small molecules [Nakatani et al. (2014)]. To our knowledge, timedependent DMRG techniques have not yet been explored with ab initio Hamiltonians, although some studies have been carried out with model Hamiltonians of conjugated systems [Ma and Schollwöck (2009)].
We will focus here on the dynamical DMRG (DDMRG) and timestep targeting timedependent DMRG (tdDMRG) methods. We concentrate on these techniques rather than the analytic DMRG response or other timedependent formulations for two reasons. First, DDMRG and tdDMRG are simple to implement in existing DMRG codes (and are thus commonly used in applications outside of quantum chemistry). Second, our work on analytic DMRG theories showed that the quality of the response functions is tied to the similarity between the excited states and the groundstate, thus excited states with quite different entanglement structure to the groundstate are poorly described except using large bond dimensions [Nakatani et al. (2014)]. Since the primary purpose of DMRG in quantum chemistry is to describe strongly correlated systems where we can often find states of different electronic character at low energies, it is of interest to work with techniques which treat states with different character in a relatively balanced way. This is the case with DDMRG and tdDMRG methods, which treat the response wavefunction or timeevolved state on an equal footing with the groundstate or initial state. In particular, we will introduce two small improvements to the techniques, that we call DDMRG and tdDMRG. Although the change to the algorithms is small and easy to implement within existing DDMRG and tdDMRG codes, the subsequent improvement in accuracy and concomitant savings in cost is significant.
The outline of the paper is as follows. In the section 2 we give a brief overview of linear response theory dynamical correlation functions as well as frequencydependent and timedependent algorithms to compute Green’s functions. We subsequently give some background on DMRG and MPS, before discussing the detailed theory of the DDMRG and tdDMRG algorithms, as well as their DDMRG and tdDMRG improvements. In section 3 we benchmark DDMRG and tdDMRG on small systems which can be exactly treated by full configuration interaction. We next use DDMRG to compute the O core excitation energy of the water molecule in realistic basis sets. Finally, we use DDMRG to compute the LDOS and gaps of hydrogen chains up to \ceH_50 within a minimal basis, and further use tdDMRG to obtain the complex polarization function to characterize the metallicity of the groundstate as a function of bondlength. We finish with some perspectives in section 4.
2 Theoretical Methods
2.1 Linear response
When the applied fields are not too strong, linear response theory underpins spectroscopy. We briefly recap the essentials here. Consider a system in an initial eigenstate of a Hamiltonian , and consider a timedependent perturbation , where is the field strength. The linear response of the observable is given by
(1) 
where and the Kubo formula[Kubo (1957)] for the generalized susceptibility is:
(2) 
The frequency dependent susceptibility is:
(3)  
where is a infinitesimal positive number, are excited states of the system, are the associated eigenvalues. The imaginary part of the susceptibility is the spectral function, which is proportional to the rate of absorption of the applied field [Coleman (2015)],
(4) 
Different spectroscopies are described by different combinations of the operators and . For example, optical spectroscopy is described by , where is the dipole operator. Likewise, photoelectron spectroscopy can be described by the retarded Green’s function,
(5) 
where respectively, are creation/annihilation operators, and is the anticommutator. Its Lehmann representation reads
(6) 
The spectral function or density of states (LDOS) becomes
(7)  
In this work, we will focus on the Green’s function and density of states as measured by photoelectron spectroscopy, but the formalism can easily be extended to other spectroscopies.
2.2 Frequency and timedomain calculations of Green’s functions
We can obtain equivalent information on the linear response in the frequency and in the timedomain. We now discuss general strategies to compute the Green’s function in these two settings. Notice that the Green’s function has two contributions, see Eq. (6). The first part corresponds to the electron addition (EA) component of the Green’s function, while the second part corresponds to the electron removal (IP) one. Computationally, we can compute the two pieces separately. Below we present explicit formulae only for the IP part, and analogous derivations apply to the EA part.
Formally, the frequency ()dependent IP Green’s function matrix element (6) can be rewritten as,
(8) 
It is convenient to compute the Green’s function from the response equation:
(9) 
where is referred to as the correction vector [Pati et al. (1999), Kühner and White (1999)], such that the Green’s function element is the expectation value
(10) 
Using real arithmetic, we solve for the real () and imaginary parts () of the correction vector separately. To compute the imaginary part from the equation,
(11) 
we can in general minimize the Hylleraaslike functional [Jeckelmann (2002)],
(12) 
From the imaginary part, the real part can be obtained as:
(13) 
In the time () domain the IP part of the Green’s function (5) is written as:
(14) 
The steady state Green’s function is obtained at sufficiently long time . From this, the frequency dependent Green’s function (8) can be obtained by Fourier transform,
(15) 
Eq. (14) can be evaluated by a realtime propagation of an initial state (). There are many methods to carry out the timepropagation [Vidal (2003), Cazalilla and Marston (2002), Daley et al. (2004), Feiguin and White (2005), White and Feiguin (2004), Zaletel et al. (2015)]; in this work we use the simple RungeKutta (RK4) algorithm, which requires calculating four vectors:
(16) 
where is the wavefunction at the initial timestep and is the timestep. From these four vectors the state at time can then be obtained as:
(17) 
The total accumulated timestep error is .
We will next see how to translate these general expressions to compute Green’s functions in the language of DMRG.
2.3 DMRG and MPS
To lay some foundations for the timedependent algorithms, we recall the main ideas of DMRG and Matrix Product States (MPS). For details, the reader is referred to the recent reviews, see Refs. 43; 50 and 51. The MPS is the underlying variational wavefunction ansatz used in DMRG algorithms, and is a nonlinear parametrization for the wave function of the form:
(18) 
where is an occupation vector in the Fock space, and is an matrix of numbers, while and are and vectors. For a given occupancy vector, the product of matrices (with vectors for the leftmost and rightmost sites) yields the scalar wavefunction amplitude. is the bond dimension (also known as the number of renormalized states) of the DMRG wavefunction. As (or in a finite Fock space , ) the MPS becomes an exact representation of any state.
In the most general sense, the DMRG algorithm provides a way to determine the tensors in the MPS one by one from to (holding all other tensors fixed at each step) from the variational principle, or equivalently the minimization of the Lagrangian,
(19) 
One such determination of all the tensors (going forwards and backwards) is called a sweep. Note that the tensors are not unique because of the product form of the MPS; gauge matrices may be inserted in between the tensors while keeping the state invariant. To properly condition the optimization, when optimizing the th tensor, we use the socalled mixed canonical gauge around site :
(20) 
where the tensors to the left and right of satisfy the orthogonality conditions respectively:
(21) 
Because of the orthogonality conditions, the and tensors collectively define orthogonal sets of manyparticle renormalized bases, recursively,
(22) 
and the MPS wavefunction may be equivalently written in the space of these renormalized states as:
(23) 
where the symbol indicates that the wave function is in the mixed canonical form at site . At each site in a DMRG sweep one performs several operations: constructing the renormalized bases and the renormalized operators in these bases at each site (blocking); determining the site wavefunction (solving); and transforming all quantities to the canonical form of the next site (decimation).
For example, in the groundstate DMRG algorithm, at each site , we build the renormalized site Hamiltonian () by projecting the Hamiltonian () into the renormalized basis of the site:
(24) 
where projects into the basis . Then, Eq. (19) becomes a quadratic function in . We then solve for the groundstate of through:
(25) 
which amounts to a standard eigenvalue problem for in Eq. (23). The final step is to transform all quantities to the mixed canonical gauge at the neighbouring site. We do so by building the density matrix in the blocked basis with matrix elements:
(26) 
where is the reshaped matrix from the tensor . The eigenvectors of the with the largest eigenvalues form a matrix with elements ; when reshaped to this becomes the tensor that replaces in the MPS. A guess for the sitewavefunction at site can be obtained by transforming Chan and HeadGordon (2002):
(27) 
where both and are the matrix versions of the site tensors, respectively.
In many DMRG algorithms, one is interested in simultaneously representing multiple states as matrix product states. It can be convenient computationally to constrain these MPS such that different states use the same renormalized bases at each site; then each state is distinguished only by its respective site wavefunction . Such algorithms are known as stateaveraged algorithms. To construct the common renormalized bases at each site, one transforms bases from site to site via the “stateaveraged” density matrix:
(28) 
where are weights and are the density matrices of the individual states entering into the average computed using Eq. (26). In this case, the density matrix has more than nonzero eigenvalues and the transformation from site to site does not precisely preserve the states unless . For finite this requires choosing a site at which to compute observables. In our case, we report observables calculated at the middle of the sweep, although other choices are possible.
Finally, we mention that in the following sections, the action of an operator on an MPS will be frequently encountered (e.g. on the right hand side of Eq. (11)). In certain cases, it is necessary to reduce the bond dimension of the state , for example in the variational compression used in the benchmark tdDMRG(G) algorithm below, or if one needs to use a smaller bond dimension in the DDMRG calculation than in the groundstate DMRG calculation. The reduction in bond dimension can in general be achieved via a variational compression by constructing the “leastsquares” functional,
(29) 
Similar to the minimization of Eq. (19) for the ground state, the MPS representation for can be obtained by minimizing this functional using analogous DMRG sweeps. The only difference is that instead of solving an eigenvalue problem (25), a linear equation needs to be solved at each site , whose solution in the mixed canonical form is simply given by the local projection .
2.4 Ddmrg
We now discuss how to determine the frequencydependent Green’s function using MPS and the dynamical DMRG (DDMRG) algorithm. As discussed earlier, the DDMRG algorithm has proven to be one of the most accurate methods to compute Green’s functions and other frequency dependent correlation functions within a MPS representation. We earlier studied its performance for chemical problems in Ref. 17. First, we recap the algorithm and then describe a modification to improve its formal properties and accuracy which we term DDMRG.
The basic path to transcribe the equations in Sec. 2.2 into a DMRG algorithm is to translate each equation to the wavefunctions and operators at each site of the DMRG sweep. The states and operators are then expressed in the renormalized basis at site . The simplest choice is to work with a stateaveraged formalism, such that all MPS share the same renormalized basis at each site. In the standard DDMRG algorithm, we first solve equation (25) at site for the groundstate wavefunction . Then, we solve for the correction vector at each site, where in Eq. (12) we additionally use the projected quantities and . Note that the Hylleraas functional of Eq. (12) involves the square of the Hamiltonian operator, and , but this approximation becomes exact in the limit . To ensure that all states continue to share the same renormalized basis throughout the sweep, we construct the density matrix for the decimation using equally weighted contributions from , , , .
The accuracy of the DDMRG procedure is controlled by the bond dimension . This governs the quality of the representation of the states such as and , as well as the quality of the resolution of the identity approximation for . In a finite system, the imaginary factor can be chosen arbitrarily, but a smaller leads to more iterations in minimizing the Hylleraas functional, and a larger bond dimension is needed to represent accurately.
Despite the established power of the DDMRG, there are a few drawbacks to the algorithm, some of which we discussed in Ref. 17. These stem from the use of the stateaveraged formalism, which means that some accuracy in the representation of each state is lost for a given bond dimension . For example, the groundstate wavefunction in DDMRG for a given is less accurate than that obtained in the standard groundstate DMRG algorithm. A related sideeffect is that even after completing a groundstate DMRG calculation, it is necessary to reoptimize the (worse) groundstate in DDMRG to accommodate the new renormalized basis. For these reasons, we have modified the original dynamical DMRG algorithm to avoid these problems; we term the modified algorithm, DDMRG. Roughly speaking, we allow each of the states appearing in the response equation to be an independent MPS (and thus to generate its own renormalized basis at each site ). More precisely, to avoid complex MPS tensors, we keep , as independent MPS, and the pair , are represented within a common renormalized basis. This means that we can reuse the solution of a groundstate DMRG sweep as and there is no loss of accuracy in the groundstate representation during the DDMRG sweeps. The modified DDMRG scheme can be summarized as follows:

A groundstate DMRG calculation is carried out to obtain and the MPS .

We compute a separate MPS, .

We carry out the DDMRG sweep where we minimize the functional in Eq. (12) at each site using the conjugate gradient algorithm. At site , this gives the correction vectors , .

and are averaged in the density matrix, which is used to transform all quantities to the next site in the sweep, and the sweeps are iterated until convergence.
2.5 tdDMRG
The timedependent DMRG (tdDMRG) algorithm that we will discuss was introduced by Feiguin and White and belongs to the family of adaptive timedependent DMRG (tdDMRG) methods. It is based on the 4 order Runge Kutta (RK4) algorithm described in Sec. 2.2. The advantage of this tdDMRG algorithm is that it is quite simple to implement for Hamiltonians with nonlocal interactions (as relevant for quantum chemistry) within a standard DMRG program. We first describe Feiguin and White’s tdDMRG algorithm and then describe an improvement to this algorithm that we will call tdDMRG.
As discussed, we can adapt the formalism in Sec. 2.2 to a DMRG algorithm by carrying out each step within the renormalized basis at each site. Again, the simplest procedure to implement is to use a stateaveraged formalism, where all MPS appearing in the equations share the same renormalized basis at site . Thus the four RungeKutta vectors in Eq. (2.2) become vectors in the space of site , , and the Hamiltonian used to construct the vectors is . Note that higher powers of are used in constructing the RungeKutta vectors. Similarly to as in DDMRG, we invoke the approximation:
(30) 
which again, introduces an error which only vanishes in the limit of infinite bond dimension. The final consideration is the decimation step to transform from one site to the next. In tdDMRG, this is done by first computing wavefunctions at the intermediate times and using linear combinations of the vectors:
(31) 
The density matrix used for the renormalization is the weighted average of all the (site) wavefunctions at different times:
(32) 
Feiguin and White Feiguin and White (2005) found by experimentation that the choice of weights
(33) 
gave the best convergence with bond dimension during the timepropagation.
The accuracy of a tdDMRG simulation is controlled by the bond dimension as well as the timestep and total propagation time . In general, it is found that as increases, the bond dimension needs to be increased to maintain accuracy in the wavefunction, due to the generic growth of entanglement during time evolution. Decreasing the timestep decreases the RungeKutta integration error, however, it also increases the number of DMRG sweeps and thus the number of compressions of the wavefunction which can also lead to an accumulated error. Feiguin and White (2005) Consequently, the timestep should be chosen to balance the intrinsic timeintegration error with the error due to DMRG compressions.
Similarly to DDMRG, the use of a stateaveraged renormalized basis at each site introduces some undesirable errors into the tdDMRG algorithm. For example, the MPS at the beginning of a timestep, represented in the renormalized basis at time , becomes approximated by the renormalized basis at time at the end of the timestep, introducing an error in the representation of the initial state. Thus, we now consider a more accurate method, where states at different times are represented by independent MPS. In the most general extension, every state appearing in the RungeKutta scheme would be represented by its own independent MPS, i.e. , , and the RungeKutta vectors . Operations that increase the bond dimension of the MPS (e.g. when applying the Hamiltonian to construct the RungeKutta vectors, or adding the RungeKutta vectors to obtain are then followed by variational MPS compression to the desired bond dimension. We call this scheme, which corresponds to the most direct implementation of time evolution with MPS in the RungeKutta context, tdDMRG(G), to denote the general extension. However, this scheme is significantly more expensive due to the many compression steps. A practical compromise is to retain only independent renormalized bases for and , and to make use of approximations such as Eq. (30) to reduce the cost. We call this method tdDMRG. In this case, we construct the four RungeKutta states as:
(34) 
where projects onto the renormalized basis of at site , and projects onto the renormalized basis of at site . The two sets of renormalized bases and are transformed to site using the density matrices of and respectively. More precisely, we use the state average of the density matrices from the real and imaginary parts of the wavefunctions, to ensure that all tensors in the MPS are real. Note that if we carried out timepropagation using a first order timestep scheme (involving only the first RungeKutta vector ) then the above procedure is the same as tdDMRG(G), as introduces no error, and can viewed as the variational MPS compression (up to the detail of averaging the real and imaginary wavefunction contributions to the density matrix). At the RK4 level, additional errors beyond tdDMRG(G) are introduced into the higher RungeKutta vectors. However, additional compressions are avoided by reusing the projected Hamiltonian to construct the additional vectors. Importantly, the cost of the tdDMRG method is only a factor of two higher than the standard tdDMRG procedure of Feiguin and White for blocking and renormalization of the operators, but as we shall see in the following section, it gives rise to significant improvements in accuracy for a fixed bond dimension, allowing for time savings in practice.
In summary, the tdDMRG algorithm consists of:

Carrying out groundstate DMRG to obtain and .

Computing the MPS for .

Propagating in realtime for a total time () as required for the desired accuracy in the spectrum. The propagation scheme consists of sweeps for each timestep. At each site , we compute the four RungeKutta vectors using the site Hamiltonians and as in Eqs. (2.5). We update the renormalized basis for using the eigenvectors of the density matrix built from . Sweeps are carried out until convergence in (typically 24 sweeps are sufficient).

If desired, is Fourier transformed using Eq. (15) to obtain the frequencydependent Green’s function.
3 Results and Discussion
3.1 Benchmarking DDMRG and tdDMRG
The DDMRG and tdDMRG algorithms above have been implemented inside the Block DMRG code. Chan and HeadGordon (2002); Chan (2004); Ghosh et al. (2008); Sharma and Chan (2012) We now examine the performance of the DDMRG and tdDMRG algorithms in the context of two simple systems where exact results can be computed. The first is a 10 atom equally spaced hydrogen chain at the equilibrium bond distance ( (Bohr)) using a minimal STO6G basis set.Hehre et al. (1969) We shall return to the hydrogen chain problem in more detail in Section 3.3. The second is an 8 site 1D Hubbard model with . Except where otherwise stated, we will use spinadapted implementations of the algorithms. We found that, similarly to groundstate simulations, spinadaptation provides roughly a factor of two gain in the effective bond dimension (see Supplementary Material).
Here we first analyze the performance of DDMRG and tdDMRG in the context of the \ceH_10 hydrogen chain. Shown in Fig. 1 is the LDOS () ( a.u.) computed with FCI compared against DDMRG and tdDMRG ( a.u., a.u.). LDOS have been calculated in this case at the central site of the chain starting from converged DMRG calculations (=500), and calculations are done in the Löwdin orthogonalized basis. To simplify visual comparisons only the IP part of the LDOS is presented here.
From Fig. 1, we see that both DDMRG and tdDMRG approach the reference FCI result as is increased towards the maximum value (=100). However, DDMRG converges much more quickly than tdDMRG toward the exact result. In particular, DDMRG is indistinguishable from FCI already at =30, while tdDMRG requires =50100 to reach the same accuracy. At =30, the tdDMRG spectrum also has small unphysical negative parts in the frequency region between 0.5 and 0.3 a.u.. The higher accuracy of the DDMRG is to be expected given that the algorithm targets a single frequency at a time.
Analyzing the computational cost of the two algorithms we have found, for the used, that the total cost of the DDMRG and tdDMRG calculations (i.e. over all frequencies and for the total propagation time) to reach a given accuracy is quite similar. However in many molecular applications, only a small range of frequencies is of interest. In that case DDMRG is particularly efficient, as tdDMRG computes the spectra over the whole frequency range. Further, the DDMRG calculations can be carried out independently for each frequency point, allowing for easy parallelization.
Both DDMRG and tdDMRG are evolutions of their parent algorithms because they do not restrict all MPS appearing in the equations to share the same stateaveraged basis. We now examine the effect of this improvement. In Fig. 2 we compare the DDMRG and DDMRG algorithms for the 10 site hydrogen chain.
While both agree at larger bond dimension (as they must) for the smaller bond dimension () we see that the DDMRG spectrum is significantly improved over the DDMRG spectrum, and in particular the DDMRG spectrum it oscillates, and this is a consequence of representing the groundstate wavefunction by an MPS in a stateaveraged basis with only a small bond dimension. In contrast, even if we use an groundstate MPS in the DDMRG algorithm, it has a consistent converged energy across the sweep which gives rise to a much more stable spectrum.
In Fig. 3 we compare the tdDMRG and tdDMRG algorithms for the 10 site hydrogen chain (\ceH_10). We see that the tdDMRG calculation is comparable in accuracy, if not better than, the tdDMRG calculation.
In both the DDMRG and tdDMRG cases, the cost of the calculations for fixed bond dimension is roughly twice the cost of the original DDMRG and tdDMRG algorithm. On the other hand, the effective bond dimension in DDMRG and tdDMRG appears to be close to twice that in DDMRG and tdDMRG respectively. Given that the scaling with bond dimension is like , we see that the DDMRG and tdDMRG algorithms offer significant savings in practice.
Additional understanding of the behaviour of tdDMRG can be obtained comparing the timedependent Green’s function matrix elements ( in this case) calculated with tdDMRG, tdDMRG, and tdDMRG(G) using with both a linear propagator, and the 4th order RungeKutta propagator. Because of the cost of the tdDMRG(G) algorithm, which requires variational MPS compression at each time step, we performed comparisons for the simpler case of the 8site Hubbard chain. Plots of the errors calculated against the exact FCI propagation are presented in Fig. 4.
The tdDMRG(G) calculations were carried out with a general purpose MPO/MPS library without
spinadaptation Li and Chan (2017), and thus all calculations in the figure did not use spin adaptation.
We see that both with the linear propagator and the 4th order propagator, the use of more flexible renormalized bases in tdDMRG
and tdDMRG(G) significantly increases the accuracy of the propagation over the simple tdDMRG scheme. In particular, tdDMRG
roughly allows for a doubling of the propagation time over tdDMRG before a noticeable buildup of error occurs, while
tdDMRG(G) allows for a further doubling.
In the case of the linear propagator, the only difference in principle between tdDMRG and tdDMRG(G) is the use
of the real and imaginary averaged density matrix to determine the renormalized bases (compression) for the wavefunction
at the next timestep (first case), rather than the exact MPS variational compression algorithm (latter method),
3.2 Coreionization potential of \ceH_2o
As a chemical application of the methods developed here we now consider the calculation of a coreionization potential. Core spectra are generally challenging to simulate as they need a flexible treatment of electron correlation as well as the inclusion of relativistic effects Coriani and Koch (2015); Dutta et al. (2015); Wenzel et al. (2014); Brabec et al. (2012). Here, we use DDMRG to calculate the ionization potential (IP) for the deepest core orbital (O ) of water examining the basis set effects and the effects of relativity. We compare against coupled cluster calculations Dutta et al. (2015); Coriani and Koch (2015), as well as experimental reference values in Table 1.
CVS  EOM  EOM  UGA  

Basis  CCSD^{a}  CCSD  CCSD(2)^{c}  SUMRCC^{b}  DDMRG  Exp.^{d} 
ccpVDZ  543.34  543.27^{b}  541.97  542.13  539.78  
ccpVTZ  540.68  540.66^{b}  539.02  539.62  
ccpVQZ  539.73  
ccpCVDZ  542.69^{c}  541.17  541.30  
ccpCVTZ  541.15  541.13^{c}  540.03  540.10  
ccpVDZDK^{◊}  542.53  
ccpVTZDK^{◊}  539.96  
ccpVQZDK^{◊}  540.16 
We estimate the IP from a DDMRG calculation by fitting three points around the excitation peak with a parabola and extracting the position of the maximum. We used an grid of 0.01 hartree and an value of 0.05 hartree. We used a bond dimension large enough to converge the DMRG energy below the milliHartree () level (=1000 for DZ basis sets and =2000 for TZ and QZ basis sets), while a bond dimension =500 has been used in DDMRG to represent the and wave functions. Calculations using smaller bond dimensions in the ccpVQZ basis indicate that our IP results are converged to better than 0.1 eV. Smaller errors are expected for the smaller basis sets.
Overall, our computed IP’s are in general agreement with previous theoretical results and, if we use a basis set larger than double zeta (DZ), they are in good agreement with the experimental value as well. As noted above, relativistic effects are important for this quantity. Four component relativistic DMRG calculations have previously been reported in Ref. 65; here we estimate scalar relativistic corrections through the sfX2C Hamiltonian. Liu (2010); Saue (2011); Peng and Reiher (2012) The inclusion of scalar relativistic effects increases the IP by 0.350.4 eV. The final result in the largest ccpVQZ basis including scalar relativistic effects is within 0.4 eV of the experimental value. The corevalence basis sets shift the ionization potential by a similar amount but with a different sign at the DZ and TZ level.
The DDMRG calculations allow for an assessment of correlation effects beyond those treated in earlier methods. Comparing to the EOMCCSD and CVSCCSD results, we find that the correlation effects beyond doubles amount to approximately 1 eV in the IP. Interestingly, the EOMCCSD(2)* method recently developed by Dutta et alDutta et al. (2015) performs very well, with errors of roughly 0.1 eV. MRCC (UGASUMRCC) calculations, as performed by Sen et al in Ref. 59 also improve on the EOMCCSD results.
3.3 Hydrogen Chains
We now use the methods developed in this work to study longer hydrogen chains. 1D equally spaced hydrogen chains were introduced in Ref. 66 as a simple model for strong correlation in an ab initio system, with the tuning parameter being the spacing between the atoms (here denoted ). They have since become a popular model system on which to benchmark strong correlation methods Stella et al. (2011); Tsuchimochi and Scuseria (2009); Lin et al. (2011); Sinitskiy et al. (2010); Mazziotti (2011); Motta et al. (2017), and have also spawned the study of analogous ring systems with heavier atoms Fertitta et al. (2014); Wouters et al. (2016). In the thermodynamic limit, the chains are thought to undergo a metalinsulator transition with the metallic phase being found at short bond distances and a Mott insulator found at long distances. 1D hydrogen chains also serve as a dimensionally reduced setting to study the hydrogen phase diagram, which is of particular interest in understanding the high pressure interiors of planets such as Jupiter and Saturn.
The metalinsulator transition in hydrogen chains can be identified in terms of different observables. Direct evidence can be obtained by computing the bandgap in the thermodynamic limit, which must vanish for a metal. Alternatively, groundstate correlation functions can be computed. For a 1D system, the delocalization of the electrons associated with the metallic phase can be established by the vanishing of the manybody complex polarization function Resta (1994, 1998); Stella et al. (2011); Hine and Foulkes (2007). Also, the algebraic decay of the offdiagonal elements of the singleparticle density matrix can also be used to establish the metallic phase Hachmann et al. (2006). This latter criterion was used in earlier DMRG studies to characterize the metallicity of hydrogen chains at different bond lengths Hachmann et al. (2006).
Here we use the DDMRG and tdDMRG algorithms to calculate the LDOS and the complex polarization function respectively as measures of metallicity, as a function of bond length for three different hydrogen chains in the minimal STO6G basis set Hehre et al. (1969) with open (OBC) and periodic boundary conditions (PBC). We also carry out groundstate DMRG and restricted and unrestricted HartreeFock calculations to further support the results. All DMRG calculations are carried out with localized Löwdin orthogonalized atomic orbitals, and LDOS are presented at one of the (two) central atoms of the chain. The PBC Hamiltonian is defined using a periodic Coulomb interaction only along the chain (1D periodicity).
In Fig. 5 we present the DDMRG LDOS at three bond distances, for 10, 30, and 50 atom hydrogen chains using open boundary conditions. For these systems is close to the equilibrium bond distance.Hachmann et al. (2006); Motta et al. (2017) The PBC spectral functions for \ceH_50 at two different geometries () are also shown. Additional OBC LDOS e.g. for intermediate bond distance can be found in the Supplementary Material.
As the chainlength is increased, the gap is reduced but does not yet close. The finite size effects for \ceH_50 at and are well converged as one can observe by comparing the \ceH_30 and \ceH_50 LDOS. However, significant finite size effects start appearing for more compressed chains, as can be seen for the chain. We note that for compressed chains, the finite size error is mainly a singleparticle effect rather than a result of Coulomb interactions. This is because the kinetic energy scales as at small while the Coulomb energy scales as .
The DDMRG bandgap decreases significantly as the bondlength is decreased from 3.6 to 1.4 . As the broadening in the LDOS blurs the gap, it is difficult to determine the gap with high precision purely from the LDOS. For this reason we also show positions of the ionization potential (IP) and electron affinity (EA) (vertical lines) computed from groundstate DMRG calculations at the same geometries. Note that (up to finite bond dimension errors) these will occur at precisely the same position as the rightmost and leftmost energy poles of the IP and EA Green’s function computed from DDMRG. Determining the gap from IPEA for the \ceH_50 chain, for instance, gives 202, 209, and 530 gaps for respectively.
The finite chain gaps with OBC and PBC are not entirely consistent, and unfortunately it is difficult to estimate the band gaps in the thermodynamic limit. With PBC in particular, there are spurious interactions between charges and the periodic images of their exchangecorrelation holes, and this produces larger finite size effects in the PBC calculations than in the OBC calculations, leading to a very poor thermodynamic limit extrapolation with PBC. Note that both the finite chain OBC and PBC gaps start to increase at very compressed distances due to the large singleparticle finite size effects discussed above.
To understand the effect of correlation on the metallicity, we show for comparison the RHF and UHF results. Both the RHF and UHF solutions display gaps, and at short distances, the RHF gap agrees well with the DMRG gap; the RHF gap at for \ceH_50, for instance, is very similar to the DMRG reference (RHF = 234 , DMRG = 202 ). At longer distances, the RHF gap is too small and is only 175 at while the DMRG gap is 530 . The behaviour of the UHF gap with bond distance is qualitatively correct, but UHF overestimates the gap at all distances (e.g. for \ceH_50 at it is 312 while at it is 734 ). Note that at longer distances, the RHF gap is not a simple finite size effect but arises from the dimerization of the RHF solution through a bondorder wave, as can be clearly seen from the offdiagonal bondorder matrix elements of the 1particle density matrix (i.e. , ) see Fig. 7. The DMRG gaps are bounded by the RHF and UHF gaps for .
Another way to characterize the metallicity of the groundstate is from the complex polarization function. This quantity, denoted ,Resta (1998); Stella et al. (2011) is defined as:
(35) 
where is the component of the electron position vector along the chain axis ( in this case) and is the longitudinal dimension of the supercell. measures electron delocalization in the groundstate and its modulus for metallic behaviour, while in an insulator. Although is a complicated manybody observable, it can be conveniently computed by carrying out a time evolution for unit time using the fictitious Hamiltonian , followed by evaluating the overlap with the groundstate. Here we compute using the tdDMRG algorithm. Note that when PBC are imposed the direct calculation of dipole integrals is not possible.Resta (1998) Given the local character of the Gaussian basis used, we define the dipole integrals as a multiplicative operator over the basis functions of the reference cell, such that: where is the dimensionless number that indexes the position of the site on the chain. In the metallic limit, where the wavefunction is a product state of Bloch functions built from a single atom unit cell, this approximation yields as an exact evaluation would, and further the approximation becomes exact in the limit of long bond distances.
In Fig. 8 we plot the DMRG complex polarization function for \ceH_10, \ceH_30, and \ceH_50 with PBC; for the \ceH_50 chain we compare with the RHF and UHF values. The absolute value of the complex polarization function is exponentially sensitive to localization length and decreases very rapidly, for \ceH_50 for instance, near , and becomes close to zero for . A similar picture is presented by the RHF and UHF complex polarization functions. Unlike the singleparticle gap, the complex polarization function can vanish in a system even when singleparticle finite size effects are large so long as the electrons are completely delocalized. The vanishing of the DMRG complex polarization function in this system at short distances, as also reflected by the similarity in the size of the gaps, thus reflects the fact that the DMRG wavefunction begins to resemble the RHF wavefunction which is a Slater determinant of planewave like orbitals. However, the scaling of the complex polarization function with system size, much like the gap, converges only slowly with system size. Thus, to definitively establish a metal insulator transition will require studies of larger systems. These studies will be discussed in a future publication.
4 Conclusions
In this work we studied two algorithms to obtain dynamical quantities from density matrix renormalization group wavefunctions in the ab initio context: the dynamical DMRG (DDMRG) algorithm, and the timestep targeting timedependent DMRG (tdDMRG) algorithm. In particular, we proposed and implemented two improved variants of these algorithms, DDMRG and tdDMRG, in the context of computing Green’s functions and the density of states. DDMRG and tdDMRG yield improved dynamical quantities with respect to their parent DDMRG and tdDMRG algorithms, at a nominal increase in cost, and they are both simple to implement within existing ab initio DMRG codes. Our analysis suggests that DDMRG and tdDMRG require a comparable amount of computation time if we desire the density of states at a similar resolution over a large energy window. However, if one is interested only in the density of states in a small energy window (e.g. when computing the principal core ionization peak) then DDMRG is advantageous.
In our applications, we showed that in the water molecule, we could use DDMRG to compute a core excitation energy in a quadruple zeta basis at a benchmark level of quality beyond that of existing correlation treatments. This suggests that DDMRG and tdDMRG will provide benchmarking capabilities for ab initio dynamical quantities similar to that provided by groundstate DMRG for groundstate properties. We also showed in larger hydrogen chains that we could use DDMRG to compute the ab initio density of states in a system large enough to consider the thermodynamic limit of the spectrum, and used tdDMRG to compute a complicated measure of delocalization, the complex polarization function. Both these capabilities will be useful in establishing the physics of the correlated metalinsulator transition in hydrogen chains, and more broadly to approach the spectral functions of other complex condensed phase problems in the future. Finally, the feasibility of these calculations suggests that DDMRG and tdDMRG may be fruitfully used to study the correlated density of states of more complex chemical systems, such as the multicentre transition metal complexes that have previously been studied with DMRG. These are directions we will pursue in the future.
This work was supported by the US National Science Foundation via NSF:CHE1657286 and NSF:CHE1650436. E.R. would like to thank Dr. Alexander Yu. Sokolov for insightful discussions and Dr. Weifeng Hu for his help with the Block DMRG code.
Supplementary materials: effects of spinadaptation and timestep size on the accuracy of tdDMRG/tdDMRG, additional spectral functions of hydrogen chains.
Footnotes
 A single complex density matrix is used in tdDMRG(G) during the compression step.
References
 McWeeny, R. Methods of molecular quantum mechanics; Academic press, 1992.
 Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Recent advances in wave functionbased methods of molecularproperty calculations. Chem. Rev. 2012, 112, 543–631.
 Dalgarno, A.; Victor, G. The timedependent coupled HartreeFock approximation. Proc. Roy. Soc. London A. 1966; pp 291–295.
 Casida, M. E. Response theory for molecules. Recent Advances in Density Functional Methods:(Part I) 1995, 1, 155.
 Koch, H.; Harrison, R. J. Analytical calculation of full configuration interaction response properties: Application to Be. J. Chem. Phys. 1991, 95, 7479–7485.
 Koch, H.; Jo/rgensen, P. Coupled cluster response functions. Journal Chem. Phys. 1990, 93, 3333–3344.
 Mussard, B.; Coccia, E.; Assaraf, R.; Otten, M.; Umrigar, C.; Toulouse, J. Timedependent linearresponse variational Monte Carlo. arXiv preprint arXiv:1705.09813 2017,
 Zhao, L.; Neuscamman, E. Equation of motion theory for excited states in variational Monte Carlo and the Jastrow antisymmetric geminal power in Hilbert space. J. Chem. Theory Comput. 2016, 12, 3719–3726.
 Krause, P.; Klamroth, T.; Saalfrank, P. Molecular response properties from explicitly timedependent configuration interaction methods. J. Chem. Phys. 2007, 127, 034107.
 Li, X.; Smith, S. M.; Markevitch, A. N.; Romanov, D. A.; Levis, R. J.; Schlegel, H. B. A timedependent Hartree–Fock approach for studying the electronic optical response of molecules in intense fields. Phys. Chem. Chem. Phys. 2005, 7, 233–239.
 Hoodbhoy, P.; Negele, J. Timedependent coupledcluster approximation to nuclear dynamics. I. Application to a solvable model. Phys. Rev. C 1978, 18, 2380.
 Hallberg, K. A. Densitymatrix algorithm for the calculation of dynamical properties of lowdimensional systems. Phys. Rev. B 1995, 52, R9827.
 Jeckelmann, E. Dynamical densitymatrix renormalizationgroup method. Phys. Rev. B 2002, 66, 045114.
 Ramasesha, S.; Pati, S. K.; Krishnamurthy, H.; Shuai, Z.; Brédas, J. Lowlying electronic excitations and nonlinear optic properties of polymers via symmetrized density matrix renormalization group method. Synth. Met. 1997, 85, 1019 – 1022.
 Pati, S. K.; Ramasesha, S.; Shuai, Z.; Brédas, J. L. Dynamical nonlinear optical coefficients from the symmetrized densitymatrix renormalizationgroup method. Phys. Rev. B 1999, 59, 14827–14830.
 Kühner, T. D.; White, S. R. Dynamical correlation functions using the density matrix renormalization group. Phys. Rev. B 1999, 60, 335–343.
 Dorando, J. J.; Hachmann, J.; Chan, G. K.L. Analytic response theory for the density matrix renormalization group. J. Chem. Phys. 2009, 130, 184111.
 Nakatani, N.; Wouters, S.; Neck, D. V.; Chan, G. K.L. Linear response theory for the density matrix renormalization group: Efficient algorithms for strongly correlated excited states. J. Chem. Phys. 2014, 140, 024108.
 Haegeman, J.; Osborne, T. J.; Verstraete, F. Postmatrix product state methods: To tangent space and beyond. Phys. Rev. B 2013, 88, 075133.
 Cazalilla, M. A.; Marston, J. B. TimeDependent DensityMatrix Renormalization Group: A Systematic Method for the Study of Quantum ManyBody OutofEquilibrium Systems. Phys. Rev. Lett. 2002, 88, 256403.
 Vidal, G. Efficient Classical Simulation of Slightly Entangled Quantum Computations. Phys. Rev. Lett. 2003, 91, 147902.
 Daley, A. J.; Kollath, C.; SchollwÃ¶ck, U.; Vidal, G. Timedependent densitymatrix renormalizationgroup using adaptive effective Hilbert spaces. J. Stat. Mec. 2004, 2004, P04005.
 Feiguin, A. E.; White, S. R. Timestep targeting methods for realtime dynamics using the density matrix renormalization group. Phys. Rev. B 2005, 72, 020404.
 White, S. R.; Feiguin, A. E. RealTime Evolution Using the Density Matrix Renormalization Group. Phys. Rev. Lett. 2004, 93, 076401.
 Kinder, J. M.; Ralph, C. C.; Chan, G. K. Analytic time evolution, random phase approximation, and Green functions for matrix product states. arXiv preprint arXiv:1103.2155 2011,
 Haegeman, J.; Cirac, J. I.; Osborne, T. J.; Pižorn, I.; Verschelde, H.; Verstraete, F. Timedependent variational principle for quantum lattices. Phys. Rev. Lett. 2011, 107, 070601.
 Zaletel, M. P.; Mong, R. S. K.; Karrasch, C.; Moore, J. E.; Pollmann, F. Timeevolving a matrix product state with longranged interactions. Phys. Rev. B 2015, 91, 165112.
 Holzner, A.; Weichselbaum, A.; McCulloch, I. P.; Schollwöck, U.; von Delft, J. Chebyshev matrix product state approach for spectral functions. Phys. Rev. B 2011, 83, 195115.
 Ganahl, M.; Thunström, P.; Verstraete, F.; Held, K.; Evertz, H. G. Chebyshev expansion for impurity models using matrix product states. Phys. Rev. B 2014, 90, 045144.
 Wolf, F. A.; Justiniano, J. A.; McCulloch, I. P.; Schollwöck, U. Spectral functions and time evolution from the Chebyshev recursion. Phys. Rev. B 2015, 91, 115144.
 White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 1992, 69, 2863–2866.
 White, S. R. Densitymatrix algorithms for quantum renormalization groups. Phys. Rev. B 1993, 48, 10345–10356.
 White, S. R.; Martin, R. L. Ab initio quantum chemistry using the density matrix renormalization group. J. Chem. Phys. 1999, 110, 4127–4130.
 Daul, S.; Ciofini, I.; Daul, C.; White, S. R. FullCI quantum chemistry using the density matrix renormalization group. Int. J. Quantum Chem. 2000, 79, 331–342.
 Mitrushenkov, A. O.; Fano, G.; Ortolani, F.; Linguerri, R.; Palmieri, P. Quantum chemistry using the density matrix renormalization group. J. Chem. Phys. 2001, 115, 6815–6821.
 Chan, G. K.L.; HeadGordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. J. Chem. Phys. 2002, 116, 4462–4476.
 Legeza, O.; Röder, J.; Hess, B. A. Controlling the accuracy of the densitymatrix renormalizationgroup method: The dynamical block state selection approach. Phys. Rev. B 2003, 67, 125114.
 Chan, G. K.L. An algorithm for large scale density matrix renormalization group calculations. J. Chem. Phys. 2004, 120, 3172–3178.
 Schollwöck, U. The densitymatrix renormalization group. Rev. Mod. Phys. 2005, 77, 259–315.
 Moritz, G.; Reiher, M. Construction of environment states in quantumchemical densitymatrix renormalization group calculations. J. Chem. Phys. 2006, 124, 034103.
 Marti, K. H.; OndÃk, I. M.; Moritz, G.; Reiher, M. Density matrix renormalization group calculations on relative energies of transition metal complexes and clusters. J. Chem. Phys. 2008, 128, 014104.
 Marti, K. H.; Reiher, M. The Density Matrix Renormalization Group Algorithm in Quantum Chemistry. Z. Phys. Chem. 2010, 224, 583–599.
 Schollwöck, U. The densitymatrix renormalization group in the age of matrix product states. Ann. Phys. 2011, 326, 96–192.
 Sharma, S.; Chan, G. K.L. Spinadapted density matrix renormalization group algorithms for quantum chemistry. J. Chem. Phys. 2012, 136, 124121.
 Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantum chemistry. Eur. Phys. J. D 2014, 68, 272.
 OlivaresAmaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.L. The abinitio density matrix renormalization group in practice. J. Chem. Phys. 2015, 142, 034102.
 Ma, H.; Schollwöck, U. Dynamical Simulations of Polaron Transport in Conjugated Polymers with the Inclusion of Electron Electron Interactions. J. Phys. Chem. A 2009, 113, 1360–1367.
 Kubo, R. StatisticalMechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. J. Phys. Soc. Jap. 1957, 12, 570–586.
 Coleman, P. Introduction to ManyBody Physics:; Cambridge University Press: Cambridge, 2015.
 Szalay, S.; Pfeffer, M.; Murg, V.; Barcza, G.; Verstraete, F.; Schneider, R.; Legeza, Ö. Tensor product methods and entanglement optimization for ab initio quantum chemistry. Int. J. Quantum Chem. 2015, 115, 1342–1391.
 Chan, G. K.L.; Keselman, A.; Nakatani, N.; Li, Z.; White, S. R. Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms. The Journal of Chemical Physics 2016, 145, 014102.
 Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.L. Orbital optimization in the density matrix renormalization group, with applications to polyenes and Î²carotene. J. Chem. Phys. 2008, 128, 144117.
 Hehre, W. J.; Stewart, R. F.; Pople, J. A. SelfConsistent MolecularOrbital Methods. I. Use of Gaussian Expansions of SlaterType Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657–2664.
 Li, Z.; Chan, G. K.L. Spinprojected matrix product states (SPMPS): a versatile tool for strongly correlated systems. J. Chem. Theory Comput. 2017, DOI:10.1021/acs.jctc.7b00270.
 Coriani, S.; Koch, H. Communication: Xray absorption spectra and coreionization potentials within a corevalence separated coupled cluster framework. J. Chem. Phys. 2015, 143, 181103.
 Dutta, A. K.; Vaval, N.; Pal, S. EOMIPCCSD(2)*: An Efficient Method for the Calculation of Ionization Potentials. J. Chem. Theory Comput. 2015, 11, 2461–2472.
 Wenzel, J.; Wormit, M.; Dreuw, A. Calculating corelevel excitations and xray absorption spectra of mediumsized closedshell molecules with the algebraicdiagrammatic construction scheme for the polarization propagator. Journal of Computational Chemistry 2014, 35, 1900–1915.
 Brabec, J.Y.; BhaskaranNair, K.; Govind, N.; Pittner, J.Y.; Kowalski, K. Communication: Application of statespecific multireference coupled cluster methods to corelevel excitations. J. Chem. Phys. 2012, 137, 171101.
 Sen, S.; Shee, A.; Mukherjee, D. A study of the ionisation and excitation energies of core electrons using a unitary group adapted state universal approach. Mol. Phys. 2013, 111, 2625–2639.
 Liu, W. Ideas of relativistic quantum chemistry. Mol. Phys. 2010, 108, 1679–1706.
 Saue, T. Relativistic Hamiltonians for Chemistry: A Primer. ChemPhysChem 2011, 12, 3077–3094.
 Peng, D.; Reiher, M. Exact decoupling of the relativistic Fock operator. Theo. Chem. Acc. 2012, 131, 1081.
 Coriani, S.; Koch, H. Erratum: “Communication: Xray absorption spectra and coreionization potentials within a corevalence separated coupled cluster framework” [J. Chem. Phys. 143, 181103 (2015)]. J. Chem. Phys. 2016, 145, 149901.
 Ohtsuka, Y.; Nakatsuji, H. Innershell ionizations and satellites studied by the openshell reference symmetryadapted cluster/symmetryadapted cluster configurationinteraction method. J. Chem. Phys. 2006, 124, 054110.
 Knecht, S.; Legeza, Ö.; Reiher, M. Communication: Fourcomponent density matrix renormalization group. J. Chem. Phys. 2014, 140, 041101.
 Hachmann, J.; Cardoen, W.; Chan, G. K.L. Multireference correlation in long molecules with the quadratic scaling density matrix renormalization group. J. Chem. Phys. 2006, 125, 144101.
 Stella, L.; Attaccalite, C.; Sorella, S.; Rubio, A. Strong electronic correlation in the hydrogen chain: A variational Monte Carlo study. Phys. Rev. B 2011, 84, 245117.
 Tsuchimochi, T.; Scuseria, G. E. Strong correlations via constrainedpairing meanfield theory. J. Chem. Phys. 2009, 131, 121102.
 Lin, N.; Marianetti, C. A.; Millis, A. J.; Reichman, D. R. Dynamical MeanField Theory for Quantum Chemistry. Phys. Rev. Lett. 2011, 106, 096402.
 Sinitskiy, A. V.; Greenman, L.; Mazziotti, D. A. Strong correlation in hydrogen chains and lattices using the variational twoelectron reduced density matrix method. J. Chem. Phys. 2010, 133, 014104.
 Mazziotti, D. A. LargeScale Semidefinite Programming for ManyElectron Quantum Mechanics. Phys. Rev. Lett. 2011, 106, 083001.
 Motta, M.; Ceperley, D. M.; KinLic Chan, G.; Gomez, J. A.; Gull, E.; Guo, S.; JimenezHoyos, C.; Nguyen Lan, T.; Li, J.; Ma, F.; Millis, A. J.; Prokof’ev, N. V.; Ray, U.; Scuseria, G. E.; Sorella, S.; Stoudenmire, E. M.; Sun, Q.; Tupitsyn, I. S.; White, S. R.; Zgid, D.; Zhang, S. Towards the solution of the manyelectron problem in real materials: equation of state of the hydrogen chain with stateoftheart manybody methods. arXiv preprint arXiv:1705.01608 2017,
 Fertitta, E.; Paulus, B.; Barcza, G.; Legeza, O. Investigation of metalinsulatorlike transition through the ab initio density matrix renormalization group approach. Phys. Rev. B 2014, 90, 245129.
 Wouters, S.; JiménezHoyos, C. A.; Sun, Q.; Chan, G. K.L. A practical guide to density matrix embedding theory in quantum chemistry. J. Chem. Theory Comput. 2016, 12, 2706–2719.
 Resta, R. Macroscopic polarization in crystalline dielectrics: the geometric phase approach. Rev. Mod. Phys. 1994, 66, 899–915.
 Resta, R. QuantumMechanical Position Operator in Extended Systems. Phys. Rev. Lett. 1998, 80, 1800–1803.
 Hine, N. D. M.; Foulkes, W. M. C. Localization lengths over metal to band insulator transitions. J. Phys.: Condensed. Matter 2007, 19, 506212.