Time-step targeting time-dependent and dynamical density matrix renormalization group algorithms with ab initio Hamiltonians


We study the dynamical density matrix renormalization group (DDMRG) and time-dependent density matrix renormalization group (td-DMRG) 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 td-DMRG, that give increased accuracy at the same bond dimension, at a nominal increase in cost. We apply DDMRG to obtain the oxygen core-excitation energy in the water molecule in a quadruple-zeta 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 td-DMRG 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 bond-length.

]Time-step targeting time-dependent 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 single-particle 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 frequency-domain 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 Hartree-Fock [Dalgarno and Victor (1966)], density functional theory [Casida (1995)], configuration interaction [Koch and Harrison (1991)], coupled cluster [Koch and Jo/rgensen (1990)], and Jastrow-Slater wavefunctions [Mussard et al. (2017), Zhao and Neuscamman (2016)] amongst others, using analytic derivative techniques. Dynamical quantities can also be calculated in the time-domain. Here, quantum chemical methods typically formulate the equation of motion for the wavefunction from the Dirac-Frenkel (time-dependent) 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)].

Time-propagation has also been investigated in conjunction with DMRG wavefunctions. Although a wide variety of time-propagation 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 time-evolving block decimation [Vidal (2003)], are specialized to Hamiltonians with short-range interactions on a 1D lattice. For quantum chemistry, it is necessary to work with long-range interactions, and one of the early time-dependent DMRG (td-DMRG) algorithms that supported such Hamiltonians was the time-step targeting time-dependent DMRG method [Feiguin and White (2005)]. There have also been many other important developments in time-dependent DMRG which we do not discuss here, including translating time-propagation 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 time-propagation using the time-dependent variational principle [Kinder et al. (2011), Haegeman et al. (2011)], and matrix product operator representations of the time-evolution operator with improved global time-step error [Zaletel et al. (2015)].

In the current work, we explore frequency-dependent and time-dependent 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 Head-Gordon (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), Olivares-Amaya 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 RPA-like correlation energy contributions for some small molecules [Nakatani et al. (2014)]. To our knowledge, time-dependent 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 time-step targeting time-dependent DMRG (td-DMRG) methods. We concentrate on these techniques rather than the analytic DMRG response or other time-dependent formulations for two reasons. First, DDMRG and td-DMRG 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 ground-state, thus excited states with quite different entanglement structure to the ground-state 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 td-DMRG methods, which treat the response wavefunction or time-evolved state on an equal footing with the ground-state or initial state. In particular, we will introduce two small improvements to the techniques, that we call DDMRG and td-DMRG. Although the change to the algorithms is small and easy to implement within existing DDMRG and td-DMRG 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 frequency-dependent and time-dependent algorithms to compute Green’s functions. We subsequently give some background on DMRG and MPS, before discussing the detailed theory of the DDMRG and td-DMRG algorithms, as well as their DDMRG and td-DMRG improvements. In section 3 we benchmark DDMRG and td-DMRG 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 td-DMRG to obtain the complex polarization function to characterize the metallicity of the ground-state as a function of bond-length. 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 time-dependent perturbation , where is the field strength. The linear response of the observable is given by


where and the Kubo formula[Kubo (1957)] for the generalized susceptibility is:


The frequency dependent susceptibility is:


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)],


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,


where respectively, are creation/annihilation operators, and is the anticommutator. Its Lehmann representation reads


The spectral function or density of states (LDOS) becomes


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 time-domain calculations of Green’s functions

We can obtain equivalent information on the linear response in the frequency and in the time-domain. 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,


It is convenient to compute the Green’s function from the response equation:


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


Using real arithmetic, we solve for the real () and imaginary parts () of the correction vector separately. To compute the imaginary part from the equation,


we can in general minimize the Hylleraas-like functional [Jeckelmann (2002)],


From the imaginary part, the real part can be obtained as:


In the time () domain the IP part of the Green’s function (5) is written as:


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,


Eq. (14) can be evaluated by a real-time propagation of an initial state (). There are many methods to carry out the time-propagation [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 Runge-Kutta (RK4) algorithm, which requires calculating four vectors:


where is the wavefunction at the initial time-step and is the time-step. From these four vectors the state at time can then be obtained as:


The total accumulated time-step 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 time-dependent 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 non-linear parametrization for the wave function of the form:


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,


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 so-called mixed canonical gauge around site :


where the tensors to the left and right of satisfy the orthogonality conditions respectively:


Because of the orthogonality conditions, the and tensors collectively define orthogonal sets of many-particle renormalized bases, recursively,


and the MPS wavefunction may be equivalently written in the space of these renormalized states as:


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 ground-state DMRG algorithm, at each site , we build the renormalized site Hamiltonian () by projecting the Hamiltonian () into the renormalized basis of the site:


where projects into the basis . Then, Eq. (19) becomes a quadratic function in . We then solve for the ground-state of through:


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:


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 site-wavefunction at site can be obtained by transforming Chan and Head-Gordon (2002):


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 state-averaged algorithms. To construct the common renormalized bases at each site, one transforms bases from site to site via the “state-averaged” density matrix:


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 non-zero 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 td-DMRG(G) algorithm below, or if one needs to use a smaller bond dimension in the DDMRG calculation than in the ground-state DMRG calculation. The reduction in bond dimension can in general be achieved via a variational compression by constructing the “least-squares” functional,


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 frequency-dependent 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 state-averaged 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 ground-state 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 state-averaged formalism, which means that some accuracy in the representation of each state is lost for a given bond dimension . For example, the ground-state wavefunction in DDMRG for a given is less accurate than that obtained in the standard ground-state DMRG algorithm. A related side-effect is that even after completing a ground-state DMRG calculation, it is necessary to re-optimize the (worse) ground-state 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 re-use the solution of a ground-state DMRG sweep as and there is no loss of accuracy in the ground-state representation during the DDMRG sweeps. The modified DDMRG scheme can be summarized as follows:

  • A ground-state 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 td-DMRG

The time-dependent DMRG (td-DMRG) algorithm that we will discuss was introduced by Feiguin and White and belongs to the family of adaptive time-dependent DMRG (td-DMRG) methods. It is based on the 4 order Runge Kutta (RK4) algorithm described in Sec. 2.2. The advantage of this td-DMRG algorithm is that it is quite simple to implement for Hamiltonians with non-local interactions (as relevant for quantum chemistry) within a standard DMRG program. We first describe Feiguin and White’s td-DMRG algorithm and then describe an improvement to this algorithm that we will call td-DMRG.

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 state-averaged formalism, where all MPS appearing in the equations share the same renormalized basis at site . Thus the four Runge-Kutta 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 Runge-Kutta vectors. Similarly to as in DDMRG, we invoke the approximation:


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 td-DMRG, this is done by first computing wavefunctions at the intermediate times and using linear combinations of the vectors:


The density matrix used for the renormalization is the weighted average of all the (site) wavefunctions at different times:


Feiguin and White Feiguin and White (2005) found by experimentation that the choice of weights


gave the best convergence with bond dimension during the time-propagation.

The accuracy of a td-DMRG simulation is controlled by the bond dimension as well as the time-step 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 time-step decreases the Runge-Kutta 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 time-step should be chosen to balance the intrinsic time-integration error with the error due to DMRG compressions.

Similarly to DDMRG, the use of a state-averaged renormalized basis at each site introduces some undesirable errors into the td-DMRG algorithm. For example, the MPS at the beginning of a time-step, represented in the renormalized basis at time , becomes approximated by the renormalized basis at time at the end of the time-step, 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 Runge-Kutta scheme would be represented by its own independent MPS, i.e. , , and the Runge-Kutta vectors . Operations that increase the bond dimension of the MPS (e.g. when applying the Hamiltonian to construct the Runge-Kutta vectors, or adding the Runge-Kutta 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 Runge-Kutta context, td-DMRG(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 td-DMRG. In this case, we construct the four Runge-Kutta states as:


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 time-propagation using a first order time-step scheme (involving only the first Runge-Kutta vector ) then the above procedure is the same as td-DMRG(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 td-DMRG(G) are introduced into the higher Runge-Kutta vectors. However, additional compressions are avoided by reusing the projected Hamiltonian to construct the additional vectors. Importantly, the cost of the td-DMRG method is only a factor of two higher than the standard td-DMRG 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 td-DMRG algorithm consists of:

  • Carrying out ground-state DMRG to obtain and .

  • Computing the MPS for .

  • Propagating in real-time for a total time () as required for the desired accuracy in the spectrum. The propagation scheme consists of sweeps for each time-step. At each site , we compute the four Runge-Kutta 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 2-4 sweeps are sufficient).

  • If desired, is Fourier transformed using Eq. (15) to obtain the frequency-dependent Green’s function.

3 Results and Discussion

3.1 Benchmarking DDMRG and td-DMRG

The DDMRG and td-DMRG algorithms above have been implemented inside the Block DMRG code. Chan and Head-Gordon (2002); Chan (2004); Ghosh et al. (2008); Sharma and Chan (2012) We now examine the performance of the DDMRG and td-DMRG 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 STO-6G 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 spin-adapted implementations of the algorithms. We found that, similarly to ground-state simulations, spin-adaptation provides roughly a factor of two gain in the effective bond dimension (see Supplementary Material).

Here we first analyze the performance of DDMRG and td-DMRG 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 td-DMRG ( 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.

Figure 1: Dependence of the LDOS on bond dimension for a \ceH_10 chain at  . LDOSs at the central site using DDMRG (upper panel) and td-DMRG (lower panel). A broadening () of 0.005 a.u. has been used. For ease of visualization dots and lines are used to represent the same quantity (LDOS); different bond dimensions are represented by different colors.

From Fig. 1, we see that both DDMRG and td-DMRG approach the reference FCI result as is increased towards the maximum value (=100). However, DDMRG converges much more quickly than td-DMRG toward the exact result. In particular, DDMRG is indistinguishable from FCI already at =30, while td-DMRG requires =50-100 to reach the same accuracy. At =30, the td-DMRG 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 td-DMRG 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 td-DMRG 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 td-DMRG are evolutions of their parent algorithms because they do not restrict all MPS appearing in the equations to share the same state-averaged 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.

Figure 2: Comparison between DDMRG and DDMRG in the description of the spectral function of an equally spaced 10 atom hydrogen chain near the equilibrium bond distance ( ). A broadening () equal to 0.005 a.u. has been used.

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 ground-state wavefunction by an MPS in a state-averaged basis with only a small bond dimension. In contrast, even if we use an ground-state MPS in the DDMRG algorithm, it has a consistent converged energy across the sweep which gives rise to a much more stable spectrum.

Figure 3: Comparison between td-DMRG and td-DMRG in the description of the spectral function of an equally spaced 10 sites hydrogen chain at the equilibrium bond distance ( ). A broadening () equal to 0.005 a.u. has been used.

In Fig. 3 we compare the td-DMRG and td-DMRG algorithms for the 10 site hydrogen chain (\ceH_10). We see that the td-DMRG calculation is comparable in accuracy, if not better than, the td-DMRG calculation.

In both the DDMRG and td-DMRG cases, the cost of the calculations for fixed bond dimension is roughly twice the cost of the original DDMRG and td-DMRG algorithm. On the other hand, the effective bond dimension in DDMRG and td-DMRG appears to be close to twice that in DDMRG and td-DMRG respectively. Given that the scaling with bond dimension is like , we see that the DDMRG and td-DMRG algorithms offer significant savings in practice.

Additional understanding of the behaviour of td-DMRG can be obtained comparing the time-dependent Green’s function matrix elements ( in this case) calculated with td-DMRG, td-DMRG, and td-DMRG(G) using with both a linear propagator, and the 4th order Runge-Kutta propagator. Because of the cost of the td-DMRG(G) algorithm, which requires variational MPS compression at each time step, we performed comparisons for the simpler case of the 8-site Hubbard chain. Plots of the errors calculated against the exact FCI propagation are presented in Fig. 4.

Figure 4: Errors of td-DMRG, td-DMRG and td-DMRG(G) in the estimation of the matrix element of a 8-site Hubbard chain with respect to the exact FCI propagation. Results obtained using the linear propagator and the 4th order Runge-Kutta (RK4) scheme are presented in panel and respectively. In this plot different colors refer to different methods and solid and dashed lines are used to represent the real and imaginary parts of respectively.

The td-DMRG(G) calculations were carried out with a general purpose MPO/MPS library without spin-adaptation 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 td-DMRG and td-DMRG(G) significantly increases the accuracy of the propagation over the simple td-DMRG scheme. In particular, td-DMRG roughly allows for a doubling of the propagation time over td-DMRG before a noticeable buildup of error occurs, while td-DMRG(G) allows for a further doubling. In the case of the linear propagator, the only difference in principle between td-DMRG and td-DMRG(G) is the use of the real and imaginary averaged density matrix to determine the renormalized bases (compression) for the wavefunction at the next time-step (first case), rather than the exact MPS variational compression algorithm (latter method),1 and this is responsible for the difference in accuracy. In the case of the 4th order propagation scheme, td-DMRG(G) provides an accurate representation of all the Runge-Kutta vectors. This leads to an extremely stable propagation, but at the cost of a significantly larger number of compression steps (6 more compressions per time step). Based on this analysis, we can conclude that td-DMRG provides a good compromise between accuracy in the representation of the Runge-Kutta vectors, and efficiency in practice, when carrying out real-time propagation. Note that the error due to finite is smaller than the other errors analyzed in this section and thus we have not discussed it in detail. A more detailed analysis of the errors associated with the time-step is presented in the supplementary material.

3.2 Core-ionization potential of \ceH_2o

As a chemical application of the methods developed here we now consider the calculation of a core-ionization 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.

cc-pVDZ 543.34 543.27b 541.97 542.13 539.78
cc-pVTZ 540.68 540.66b 539.02 539.62
cc-pVQZ 539.73
cc-pCVDZ 542.69c 541.17 541.30
cc-pCVTZ 541.15 541.13c 540.03 540.10
cc-pVDZ-DK 542.53
cc-pVTZ-DK 539.96
cc-pVQZ-DK 540.16
  • Scalar relativistic effects have been introduced using the sf-X2C method.Liu (2010); Saue (2011); Peng and Reiher (2012)

  • Data from Ref. 63; 55

  • Data from Ref. 59

  • Data from Ref. 56

  • Data from Ref. 64

Table 1: \ceH_2O core ionization potentials (eV). Theoretical data have been calculated at the geometry of Ref. 59.

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 cc-pVQZ 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 sf-X2C Hamiltonian. Liu (2010); Saue (2011); Peng and Reiher (2012) The inclusion of scalar relativistic effects increases the IP by 0.35-0.4 eV. The final result in the largest cc-pVQZ basis including scalar relativistic effects is within 0.4 eV of the experimental value. The core-valence 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 EOM-CCSD and CVS-CCSD results, we find that the correlation effects beyond doubles amount to approximately 1 eV in the IP. Interestingly, the EOM-CCSD(2)* method recently developed by Dutta et alDutta et al. (2015) performs very well, with errors of roughly 0.1 eV. MRCC (UGA-SUMRCC) calculations, as performed by Sen et al in Ref. 59 also improve on the EOM-CCSD 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 metal-insulator 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 metal-insulator 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, ground-state 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 many-body complex polarization function Resta (1994, 1998); Stella et al. (2011); Hine and Foulkes (2007). Also, the algebraic decay of the off-diagonal elements of the single-particle 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 td-DMRG 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 STO-6G basis set Hehre et al. (1969) with open (OBC) and periodic boundary conditions (PBC). We also carry out ground-state DMRG and restricted and unrestricted Hartree-Fock 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.

Figure 5: LDOS for three equally spaced hydrogen chains (\ceH_10 - red, \ceH_30 - green, \ceH_50 - blue) at three different bond distances ( and  ). LDOS of \ceH_50 calculated at bond distances and   with PBC (blue dashed lines) have been also included. All the LDOSs have been calculated on the central site of the chain. A broadening () of 0.05 a.u. has been used. Vertical lines are used to indicate the position of the ionization potential (IP) and electron affinity (EA) of the systems.

As the chain-length 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 single-particle 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 ground-state 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 IP-EA 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 exchange-correlation 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 single-particle finite size effects discussed above.

Figure 6: Band gaps calculated for the \ceH_50 chain.

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 bond-order wave, as can be clearly seen from the off-diagonal bond-order matrix elements of the 1-particle density matrix (i.e. , ) see Fig. 7. The DMRG gaps are bounded by the RHF and UHF gaps for  .

Figure 7: Comparison of DDMRG (dots), RHF (solid line) and UHF (dashed line) density matrix off-diagonal elements for the equally spaced \ceH_50 chain as a function of the bond distance.

Another way to characterize the metallicity of the ground-state is from the complex polarization function. This quantity, denoted ,Resta (1998); Stella et al. (2011) is defined as:


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 ground-state and its modulus for metallic behaviour, while in an insulator. Although is a complicated many-body 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 ground-state. Here we compute using the td-DMRG 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.

Figure 8: DMRG and HF complex polarization functions. In panel a) complex polarization functions for \ceH_10, \ceH_30 and \ceH_50 using DMRG are presented. In panel b) complex polarization functions for \ceH_50 at the DMRG, RHF and UHF level of theory are presented. Periodic Boundary Conditions (PBC) have been used each case.

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 single-particle gap, the complex polarization function can vanish in a system even when single-particle 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 plane-wave 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 time-step targeting time-dependent DMRG (td-DMRG) algorithm. In particular, we proposed and implemented two improved variants of these algorithms, DDMRG and td-DMRG, in the context of computing Green’s functions and the density of states. DDMRG and td-DMRG yield improved dynamical quantities with respect to their parent DDMRG and td-DMRG 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 td-DMRG 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 td-DMRG will provide benchmarking capabilities for ab initio dynamical quantities similar to that provided by ground-state DMRG for ground-state 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 td-DMRG to compute a complicated measure of delocalization, the complex polarization function. Both these capabilities will be useful in establishing the physics of the correlated metal-insulator 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 td-DMRG may be fruitfully used to study the correlated density of states of more complex chemical systems, such as the multi-centre 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:CHE-1657286 and NSF:CHE-1650436. 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 spin-adaptation and time-step size on the accuracy of td-DMRG/td-DMRG, additional spectral functions of hydrogen chains.


  1. A single complex density matrix is used in td-DMRG(G) during the compression step.


  1. McWeeny, R. Methods of molecular quantum mechanics; Academic press, 1992.
  2. Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Recent advances in wave function-based methods of molecular-property calculations. Chem. Rev. 2012, 112, 543–631.
  3. Dalgarno, A.; Victor, G. The time-dependent coupled Hartree-Fock approximation. Proc. Roy. Soc. London A. 1966; pp 291–295.
  4. Casida, M. E. Response theory for molecules. Recent Advances in Density Functional Methods:(Part I) 1995, 1, 155.
  5. Koch, H.; Harrison, R. J. Analytical calculation of full configuration interaction response properties: Application to Be. J. Chem. Phys. 1991, 95, 7479–7485.
  6. Koch, H.; Jo/rgensen, P. Coupled cluster response functions. Journal Chem. Phys. 1990, 93, 3333–3344.
  7. Mussard, B.; Coccia, E.; Assaraf, R.; Otten, M.; Umrigar, C.; Toulouse, J. Time-dependent linear-response variational Monte Carlo. arXiv preprint arXiv:1705.09813 2017,
  8. 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.
  9. Krause, P.; Klamroth, T.; Saalfrank, P. Molecular response properties from explicitly time-dependent configuration interaction methods. J. Chem. Phys. 2007, 127, 034107.
  10. Li, X.; Smith, S. M.; Markevitch, A. N.; Romanov, D. A.; Levis, R. J.; Schlegel, H. B. A time-dependent Hartree–Fock approach for studying the electronic optical response of molecules in intense fields. Phys. Chem. Chem. Phys. 2005, 7, 233–239.
  11. Hoodbhoy, P.; Negele, J. Time-dependent coupled-cluster approximation to nuclear dynamics. I. Application to a solvable model. Phys. Rev. C 1978, 18, 2380.
  12. Hallberg, K. A. Density-matrix algorithm for the calculation of dynamical properties of low-dimensional systems. Phys. Rev. B 1995, 52, R9827.
  13. Jeckelmann, E. Dynamical density-matrix renormalization-group method. Phys. Rev. B 2002, 66, 045114.
  14. Ramasesha, S.; Pati, S. K.; Krishnamurthy, H.; Shuai, Z.; Brédas, J. Low-lying electronic excitations and nonlinear optic properties of polymers via symmetrized density matrix renormalization group method. Synth. Met. 1997, 85, 1019 – 1022.
  15. Pati, S. K.; Ramasesha, S.; Shuai, Z.; Brédas, J. L. Dynamical nonlinear optical coefficients from the symmetrized density-matrix renormalization-group method. Phys. Rev. B 1999, 59, 14827–14830.
  16. Kühner, T. D.; White, S. R. Dynamical correlation functions using the density matrix renormalization group. Phys. Rev. B 1999, 60, 335–343.
  17. Dorando, J. J.; Hachmann, J.; Chan, G. K.-L. Analytic response theory for the density matrix renormalization group. J. Chem. Phys. 2009, 130, 184111.
  18. 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.
  19. Haegeman, J.; Osborne, T. J.; Verstraete, F. Post-matrix product state methods: To tangent space and beyond. Phys. Rev. B 2013, 88, 075133.
  20. Cazalilla, M. A.; Marston, J. B. Time-Dependent Density-Matrix Renormalization Group: A Systematic Method for the Study of Quantum Many-Body Out-of-Equilibrium Systems. Phys. Rev. Lett. 2002, 88, 256403.
  21. Vidal, G. Efficient Classical Simulation of Slightly Entangled Quantum Computations. Phys. Rev. Lett. 2003, 91, 147902.
  22. Daley, A. J.; Kollath, C.; Schollwöck, U.; Vidal, G. Time-dependent density-matrix renormalization-group using adaptive effective Hilbert spaces. J. Stat. Mec. 2004, 2004, P04005.
  23. Feiguin, A. E.; White, S. R. Time-step targeting methods for real-time dynamics using the density matrix renormalization group. Phys. Rev. B 2005, 72, 020404.
  24. White, S. R.; Feiguin, A. E. Real-Time Evolution Using the Density Matrix Renormalization Group. Phys. Rev. Lett. 2004, 93, 076401.
  25. 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,
  26. Haegeman, J.; Cirac, J. I.; Osborne, T. J.; Pižorn, I.; Verschelde, H.; Verstraete, F. Time-dependent variational principle for quantum lattices. Phys. Rev. Lett. 2011, 107, 070601.
  27. Zaletel, M. P.; Mong, R. S. K.; Karrasch, C.; Moore, J. E.; Pollmann, F. Time-evolving a matrix product state with long-ranged interactions. Phys. Rev. B 2015, 91, 165112.
  28. 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.
  29. 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.
  30. 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.
  31. White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 1992, 69, 2863–2866.
  32. White, S. R. Density-matrix algorithms for quantum renormalization groups. Phys. Rev. B 1993, 48, 10345–10356.
  33. White, S. R.; Martin, R. L. Ab initio quantum chemistry using the density matrix renormalization group. J. Chem. Phys. 1999, 110, 4127–4130.
  34. Daul, S.; Ciofini, I.; Daul, C.; White, S. R. Full-CI quantum chemistry using the density matrix renormalization group. Int. J. Quantum Chem. 2000, 79, 331–342.
  35. 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.
  36. Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. J. Chem. Phys. 2002, 116, 4462–4476.
  37. Legeza, O.; Röder, J.; Hess, B. A. Controlling the accuracy of the density-matrix renormalization-group method: The dynamical block state selection approach. Phys. Rev. B 2003, 67, 125114.
  38. Chan, G. K.-L. An algorithm for large scale density matrix renormalization group calculations. J. Chem. Phys. 2004, 120, 3172–3178.
  39. Schollwöck, U. The density-matrix renormalization group. Rev. Mod. Phys. 2005, 77, 259–315.
  40. Moritz, G.; Reiher, M. Construction of environment states in quantum-chemical density-matrix renormalization group calculations. J. Chem. Phys. 2006, 124, 034103.
  41. 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.
  42. Marti, K. H.; Reiher, M. The Density Matrix Renormalization Group Algorithm in Quantum Chemistry. Z. Phys. Chem. 2010, 224, 583–599.
  43. Schollwöck, U. The density-matrix renormalization group in the age of matrix product states. Ann. Phys. 2011, 326, 96–192.
  44. Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithms for quantum chemistry. J. Chem. Phys. 2012, 136, 124121.
  45. Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantum chemistry. Eur. Phys. J. D 2014, 68, 272.
  46. Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. The ab-initio density matrix renormalization group in practice. J. Chem. Phys. 2015, 142, 034102.
  47. 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.
  48. Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. J. Phys. Soc. Jap. 1957, 12, 570–586.
  49. Coleman, P. Introduction to Many-Body Physics:; Cambridge University Press: Cambridge, 2015.
  50. 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.
  51. 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.
  52. 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.
  53. Hehre, W. J.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of Slater-Type Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657–2664.
  54. Li, Z.; Chan, G. K.-L. Spin-projected matrix product states (SP-MPS): a versatile tool for strongly correlated systems. J. Chem. Theory Comput. 2017, DOI:10.1021/acs.jctc.7b00270.
  55. Coriani, S.; Koch, H. Communication: X-ray absorption spectra and core-ionization potentials within a core-valence separated coupled cluster framework. J. Chem. Phys. 2015, 143, 181103.
  56. Dutta, A. K.; Vaval, N.; Pal, S. EOMIP-CCSD(2)*: An Efficient Method for the Calculation of Ionization Potentials. J. Chem. Theory Comput. 2015, 11, 2461–2472.
  57. Wenzel, J.; Wormit, M.; Dreuw, A. Calculating core-level excitations and x-ray absorption spectra of medium-sized closed-shell molecules with the algebraic-diagrammatic construction scheme for the polarization propagator. Journal of Computational Chemistry 2014, 35, 1900–1915.
  58. Brabec, J.-Y.; Bhaskaran-Nair, K.; Govind, N.; Pittner, J.-Y.; Kowalski, K. Communication: Application of state-specific multireference coupled cluster methods to core-level excitations. J. Chem. Phys. 2012, 137, 171101.
  59. 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.
  60. Liu, W. Ideas of relativistic quantum chemistry. Mol. Phys. 2010, 108, 1679–1706.
  61. Saue, T. Relativistic Hamiltonians for Chemistry: A Primer. ChemPhysChem 2011, 12, 3077–3094.
  62. Peng, D.; Reiher, M. Exact decoupling of the relativistic Fock operator. Theo. Chem. Acc. 2012, 131, 1081.
  63. Coriani, S.; Koch, H. Erratum: “Communication: X-ray absorption spectra and core-ionization potentials within a core-valence separated coupled cluster framework” [J. Chem. Phys. 143, 181103 (2015)]. J. Chem. Phys. 2016, 145, 149901.
  64. Ohtsuka, Y.; Nakatsuji, H. Inner-shell ionizations and satellites studied by the open-shell reference symmetry-adapted cluster/symmetry-adapted cluster configuration-interaction method. J. Chem. Phys. 2006, 124, 054110.
  65. Knecht, S.; Legeza, Ö.; Reiher, M. Communication: Four-component density matrix renormalization group. J. Chem. Phys. 2014, 140, 041101.
  66. 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.
  67. 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.
  68. Tsuchimochi, T.; Scuseria, G. E. Strong correlations via constrained-pairing mean-field theory. J. Chem. Phys. 2009, 131, 121102.
  69. Lin, N.; Marianetti, C. A.; Millis, A. J.; Reichman, D. R. Dynamical Mean-Field Theory for Quantum Chemistry. Phys. Rev. Lett. 2011, 106, 096402.
  70. Sinitskiy, A. V.; Greenman, L.; Mazziotti, D. A. Strong correlation in hydrogen chains and lattices using the variational two-electron reduced density matrix method. J. Chem. Phys. 2010, 133, 014104.
  71. Mazziotti, D. A. Large-Scale Semidefinite Programming for Many-Electron Quantum Mechanics. Phys. Rev. Lett. 2011, 106, 083001.
  72. Motta, M.; Ceperley, D. M.; Kin-Lic Chan, G.; Gomez, J. A.; Gull, E.; Guo, S.; Jimenez-Hoyos, 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 many-electron problem in real materials: equation of state of the hydrogen chain with state-of-the-art many-body methods. arXiv preprint arXiv:1705.01608 2017,
  73. Fertitta, E.; Paulus, B.; Barcza, G.; Legeza, O. Investigation of metal-insulator-like transition through the ab initio density matrix renormalization group approach. Phys. Rev. B 2014, 90, 245129.
  74. Wouters, S.; Jiménez-Hoyos, 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.
  75. Resta, R. Macroscopic polarization in crystalline dielectrics: the geometric phase approach. Rev. Mod. Phys. 1994, 66, 899–915.
  76. Resta, R. Quantum-Mechanical Position Operator in Extended Systems. Phys. Rev. Lett. 1998, 80, 1800–1803.
  77. Hine, N. D. M.; Foulkes, W. M. C. Localization lengths over metal to band insulator transitions. J. Phys.: Condensed. Matter 2007, 19, 506212.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

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