The density matrix renormalization group for ab initio quantum chemistry

The density matrix renormalization group for ab initio quantum chemistry

Sebastian Wouters Center for Molecular Modelling, Ghent University, Technologiepark 903, 9052 Zwijnaarde, Belgium    Dimitri Van Neck Center for Molecular Modelling, Ghent University, Technologiepark 903, 9052 Zwijnaarde, Belgium
Received: date / Revised version: date
Abstract

During the past 15 years, the density matrix renormalization group (DMRG) has become increasingly important for ab initio quantum chemistry. Its underlying wavefunction ansatz, the matrix product state (MPS), is a low-rank decomposition of the full configuration interaction tensor. The virtual dimension of the MPS, the rank of the decomposition, controls the size of the corner of the many-body Hilbert space that can be reached with the ansatz. This parameter can be systematically increased until numerical convergence is reached. The MPS ansatz naturally captures exponentially decaying correlation functions. Therefore DMRG works extremely well for noncritical one-dimensional systems. The active orbital spaces in quantum chemistry are however often far from one-dimensional, and relatively large virtual dimensions are required to use DMRG for ab initio quantum chemistry (QC-DMRG). The QC-DMRG algorithm, its computational cost, and its properties are discussed. Two important aspects to reduce the computational cost are given special attention: the orbital choice and ordering, and the exploitation of the symmetry group of the Hamiltonian. With these considerations, the QC-DMRG algorithm allows to find numerically exact solutions in active spaces of up to 40 electrons in 40 orbitals.

pacs:
31.15.A-Ab initio calculations and 31.50.BcPotential energy surfaces of ground electronic states and 05.10.ccRenormalization in statistical physics

1 Introduction

At the basis of ab initio quantum chemistry lies Hartree-Fock (HF) theory Hartree (1928); Slater (1928); Fock (1926). In HF theory, a single Slater determinant (SD) is optimized by finding the set orbitals which minimize its energy expectation value. The occupancy of the HF orbitals is definite: occupied orbitals are filled with probability 1, and virtual orbitals are empty with probability 1. The exact ground state is a linear combination over all possible Slater determinants. The difference in energy between the HF solution and the exact ground state is the correlation energy. This energy is often (somewhat ambiguously) divided into two contributions: static and dynamic correlation Helgaker et al. (2000). When near-degeneracies between determinants occur, and more than one determinant is needed to describe the qualitative behaviour of a molecule, it is said to have static correlation. This type of correlation often arises in transition metal complexes or -conjugated systems, as well as for geometries far from equilibrium. It is typically resolved with only a few determinants. The Coulomb repulsion results in a small nonzero occupancy of many virtual HF orbitals in the true ground state. This effect is called dynamic correlation, and it constitutes the remainder of the energy gap.

All static and dynamic correlation can in principle be retrieved at HF cost with density functional theory (DFT). Hohenberg and Kohn have shown that the electron density provides sufficient information to determine all ground state properties, and that there exists a unique universal functional of the electron density which can be used to obtain the exact ground state density Hohenberg and Kohn (1964). Kohn and Sham rewrote the universal functional as the sum of the kinetic energy of a noninteracting system and an exchange-correlation functional Kohn and Sham (1965). This allows to represent the electron density by means of the Kohn-Sham Slater determinant, which immediately ensures correct N-representability. Unfortunately, the universal functional is unknown. Many approximate semi-empirical exchange-correlation functionals of various complexity have been proposed. Because the exact exchange-correlation functional is unknown, not all correlation is retrieved with DFT. For single-reference systems, for which the exact solution is dominated by a single SD, DFT is good in capturing dynamic correlation. For multireference (MR) systems, DFT fails to retrieve static correlation Dickson and Becke (2005).

Dynamic correlation can also be captured with ab initio post-HF methods. These start from the optimized HF orbitals and the corresponding SD, and build in dynamic correlation on top of the single SD reference. Commonly known are Møller-Plesset (Rayleigh-Schrödinger) perturbation theory Møller and Plesset (1934), the configuration interaction (CI) expansion Slater (1929); Condon (1930), and coupled cluster (CC) theory Coester (1958); Coester and Kümmel (1960); J. C̆íz̆ek (1966). These methods are truncated in their perturbation or expansion order. An important property of wavefunctions is size-consistency: the fact that for two noninteracting subsystems, the compound wavefunction should be multiplicatively separable and the total energy additively separable. CI with excitations is not size-consistent if there are more than electrons in the compound system, whereas CC is always size-consistent because of its exponential wavefunction ansatz Helgaker et al. (2000). Because these post-HF methods start from a single SD reference, they have difficulty building in static correlation. Mostly, very large expansion orders are required to retrieve static correlation.

It is therefore better to resort to MR methods for systems with pronounced static correlation. For such systems, the subset of important orbitals (the active space), in which the occupation changes over the dominant determinants, is often rather small. This allows for a particular MR solution method: the complete active space (CAS) self-consistent field (SCF) method Roos (1980); Roos et al. (1980); Siegbahn et al. (1981). From the HF solution, a subset of occupied and virtual orbitals is selected as active space. While the remaining occupied and virtual orbitals are kept frozen at HF level, the electronic structure in the active space is solved exactly (the CAS-part). Subsequently, the occupied, active, and virtual spaces are rotated to further minimize the energy. This two-step cycle, which is sometimes implemented together, is repeated until convergence is reached (the SCF-part). CASSCF resolves the static correlation in the system. Dynamic correlation can be built in on top of the CASSCF reference wavefunction by perturbation theory (CASPT2) Malmqvist et al. (1990); Andersson et al. (1992), a CI expansion (MRCI or CASCI) Buenker and Peyerimhoff (1974); Werner and Reinsch (1982); Siegbahn (1979, 1980); Brooks and Schaefer (1979), or CC theory (MRCC or CASCC) Oliphant and Adamowicz (1992); Stolarczyk (1994). For the latter, approximate schemes such as canonical transformation (CT) theory Yanai and Chan (2006) are often used.

Because the many-body Hilbert space grows exponentially with the number of single-particle states, only small active spaces, of up to 18 electrons in 18 orbitals, can be treated in the CAS-part. In 1999, the density matrix renormalization group (DMRG) was introduced in ab initio quantum chemistry (QC) White and Martin (1999). This MR method allows to find numerically exact solutions in significantly larger active spaces, of up to 40 electrons in 40 orbitals.

2 Matrix product states

The electronic Hamiltonian can be written in second quantization as

(1)

The Latin letters denote spatial orbitals and the Greek letters electron spin projections. The and are the one- and two-electron integrals, respectively. In the occupation number representation, the basis states of the many-body Hilbert space are

(2)

The symmetry group of the Hamiltonian (1) is , or total electronic spin, particle-number, and molecular point group symmetry. By defining the operators

(3)
(4)
(5)
(6)
(7)

it can be easily checked that , , , and form a set of commuting observables. This constitutes the total electronic spin and particle-number symmetries. For fixed particle number , Eq. (1) can also be written as

(8)

with

(9)

The molecular point group symmetry consists of the rotations, reflections, and inversions which leave the external potential due to the nuclei invariant. These symmetry operations map nuclei with equal charges onto each other. The point group symmetry has implications for the spatial orbitals. Linear combinations of the single-particle basis functions can be constructed which transform according to a particular row of a particular irreducible representation (irrep) of Cornwell (1984). As the Hamiltonian transforms according to the trivial irrep of , can only be nonzero if the reductions of and have at least one irrep in common. Most molecular electronic structure programs make use of the abelian point groups with real-valued character tables.

An eigenstate of the Hamiltonian (8) can be written as

(10)

The size of the full CI (FCI) tensor grows as , exponentially fast with . This tensor can be exactly decomposed by a singular value decomposition (SVD) as follows:

(11)

Define

(12)

and decompose the right unitary again with an SVD as follows:

(13)

Define

(14)

Continue by successively decomposing the right unitaries . In this way, the FCI tensor can be exactly rewritten as the following contracted matrix product:

(15)

which is graphically represented in Fig. 1. Except for the first and last orbital (or site), Eq. (15) introduces a rank-3 tensor per site. One of its indices corresponds to the physical index , the other two to the virtual or bond indices and . In Fig. 1, tensors are represented by circles, physical indices by open lines, and virtual indices by connected lines. The graph hence represents how the contracted matrix product decomposes the FCI tensor. Since no assumptions are made about the FCI tensor, the dimension of the indices has to grow exponentially towards the middle of this contracted product:

(16)

This is solely due to the increasing matrix dimensions in the successive SVDs. Instead of variationally optimizing over the FCI tensor, one may as well optimize over the tensors of its decomposition (15). To make Eq. (15) of practical use, its dimensions can be truncated:

(17)

The corresponding ansatz is called a matrix product state (MPS) with open boundary conditions. The truncation dimension is called the bond or virtual dimension. The MPS ansatz can be optimized by the DMRG algorithm White (1992, 1993); White and Martin (1999), yielding a variational upper bound for the ground state energy.

Figure 1: Tensors are represented by circles, physical indices by open lines, and virtual indices by connected lines. The MPS graph hence represents how the contracted matrix product decomposes the FCI tensor.

DMRG was invented in 1992 by White in the field of condensed matter theory White (1992). Östlund and Rommer discovered in 1995 its underlying variational ansatz, the MPS Östlund and Rommer (1995); Rommer and Östlund (1997). The discovery of the MPS ansatz allowed to understand DMRG by means of quantum information theory. The area law for one-dimensional quantum systems, see section 3, was proven by Hastings in 2007 Hastings (2007), and constitutes a hard proof that an MPS is very efficient in representing the ground state of noncritical one-dimensional quantum systems.

The MPS ansatz was in fact discovered earlier, under various names. Nishino found that they were used in statistical physics as a variational optimization technique Nishino (August 23 - September 3, 2010): in 1941 by Kramers and Wannier Kramers and Wannier (1941) and in 1968 by Baxter Baxter (1968). Nightingale and Blöte recycled Baxter’s ansatz in 1986 to approximate quantum eigenstates Nightingale and Blöte (1986). In 1987, Affleck, Kennedy, Lieb and Tasaki constructed the exact valence-bond ground state of a particular next-nearest-neighbour spin chain Affleck et al. (1987). They obtained an MPS with bond dimension 2. In mathematics, the translationally invariant valence-bond state is known as a finitely correlated state Fannes et al. (1989, 1992), and in the context of information compression, an MPS is known as a tensor train Oseledets (2011); Savostyanov et al. (2014).

The concept of a renormalization group was first used in quantum electrodynamics. The coarse-grained view of a point-like electron breaks down at small distance scales (or large energy scales). The electron itself consists of electrons, positrons, and photons. The mass and charge contributions from this fine structure lead to infinities. These were successfully resolved by Tomonaga, Schwinger, and Feynman Tomonaga (1946); Schwinger (1948a, b); Feynman (1949a, b). Later, Wilson used a numerical renormalization group (NRG) to solve the long-standing Kondo problem Wilson (1975). He turned the coupling of the impurity to the conduction band into a half-infinite lattice problem by discretizing the conduction band in momentum space. For increasing lattice sizes, only the lowest energy states are kept at each renormalization step. These are sufficient to study the low-temperature thermodynamics of the impurity system. Although very successful for impurity systems, NRG fails for real-space lattice systems such as the discretized particle-in-a-box, spin-lattice, and Hubbard models. For these systems, the low energy states of a small subsystem are often irrelevant for the ground state of the total system White and Noack (1992). Consider for example the ground state of the particle-in-a-box problem. By concatenating the solution of two smaller sized boxes, an unphysical node is introduced in the approximation of the ground state of the larger problem. It was White who pointed out this problem and resolved it with his DMRG method White (1992). Instead of selecting the degrees of freedom with lowest energy, the most relevant degrees of freedom should be selected.

3 Entanglement and the von Neumann entropy

This section attempts to clarify the broader context of DMRG. A brief introduction to quantum entanglement, the von Neumann entropy, and the area law is given.

Figure 2: Bipartition of the single-particle states.

Consider the bipartition of orthonormal single-particle states into two subsystems and in Fig. 2. Suppose and are the orthonormal basis states of the many-body Hilbert spaces of resp. subsystem and . The Hilbert space of the composite system is spanned by the product space , and a general quantum many-body state of the composite system can be written as

(18)

The Schmidt decomposition of is obtained by performing an SVD on and by rotating the orthonormal bases and with the unitary matrices:

(19)

For normalized :

(20)

For the given bipartition, one is sometimes interested in the optimal approximation of in a least squares sense . It can be shown that the optimal approximation, with a smaller number of terms in the summation (18), is obtained by keeping the states with the largest Schmidt numbers in Eq. (19). This fact will be of key importance for the DMRG algorithm (see section 4.3).

In classical theories, the sum over can contain only one nonzero value . A measurement in subsystem then does not influence the outcome in subsystem , and the two subsystems are not entangled. In quantum theories, the sum over can contain many nonzero values . State in subsytem occurs with probability , as can be observed from the reduced density matrix (RDM) of subsystem :

(21)

Analogously the RDM of subsystem can be constructed:

(22)

From (19), it follows that the measurement of in subsystem implies the measurement of in subsystem with probability 1. Measurements in and are hence not independent, and the two subsystems are said to be entangled.

Consider for example two singly occupied orbitals and in the spin-0 singlet state:

(23)

The measurements of the spin projections of the electrons are not independent. Each possible spin projection of the electron in can be measured with probability , but the simultaneous measurement of both spin projections will always yield

(24)

with probability 1.

The RDMs and allow to define the von Neumann entanglement entropy von Neumann (1927):

(25)

This quantum analogue of the Shannon entropy is a measure of how entangled subsystems and are. If they are not entangled, and for , which implies . If they are maximally entangled, for all and , which implies , with the minimum of the sizes of the many-body Hilbert spaces of and .

A Hamiltonian which acts on a -dimensional quantum lattice system in the thermodynamic limit is called local if there exists a distance cutoff beyond which the interaction terms decay at least exponentially. Consider the ground state of a gapped -dimensional quantum system in the thermodynamic limit, and select as subsystem a hypercube with side and volume . The von Neumann entropy is believed to obey an area law Plenio et al. (2005); Eisert et al. (2010); Van Acoleyen et al. (2013):

(26)

This is the result of a finite correlation length, as only lattice sites in the immediate vicinity of the hypercube’s boundary are then correlated with lattice sites on the other side of the boundary. This is a theorem for one-dimensional systems Hastings (2007) and a conjecture in higher dimensions Eisert et al. (2010), supported by numerical examples and theoretical arguments Van Acoleyen et al. (2013). For critical quantum systems, with a closed excitation gap, there can be logarithmic corrections to the area law Vidal et al. (2003); Eisert et al. (2010).

For gapped one-dimensional systems, consider as subsystem a line segment of length . Its boundary consists of two points. Due to the finite correlation length in the ground state, the entanglement of the subsystem does not increase with , if is significantly larger than the correlation length. The von Neumann entropy is then a constant independent of , and the ground state can be well represented by retaining only a finite number of states in the Schmidt decomposition of any bipartition of the lattice in two semi-infinite line segments. This is the reason why the MPS ansatz and the corresponding DMRG algorithm work very well to study the ground states of gapped one-dimensional systems.

Figure 3: Several tensor network states. Tensors are represented by circles, physical indices by open lines, and virtual indices by connected lines. The graph hence represents how the ansatz decomposes the FCI tensor.

The MPS ansatz

(27)

is shown graphically in Fig. 3. Except for the first and last orbital (or site), the MPS ansatz introduces a rank-3 tensor per site. One of its indices corresponds to the physical index , the other two to the virtual indices and . Similar to Fig. 1, tensors are represented by circles, physical indices by open lines, and virtual indices by connected lines in Fig. 3. The graph hence represents how the ansatz decomposes the FCI tensor. The finite size of the virtual indices can capture finite-length correlations along the one-dimensional chain. Stated more rigorously: for a system in the thermodynamic limit, all correlation functions measured in an MPS ansatz with finite decay exponentially with increasing site distance Fannes et al. (1992); Evenbly and Vidal (2011):

(28)

Unless the lattice size is reasonably small Stoudenmire and White (2012), an MPS is not efficient to represent the ground state of higher dimensional or critical systems. Fortunately, efficient tensor network states (TNS) for higher dimensional and critical lattice systems, which do obey the correct entanglement scaling laws, have been developed Evenbly and Vidal (2011). There even exists a continuous MPS ansatz for one-dimensional quantum field theories Verstraete and Cirac (2010).

The ansatz for two-dimensional systems is called the projected entangled pair state (PEPS) Verstraete and Cirac (2004), see Fig. 3. Instead of two virtual indices, each tensor now has four virtual indices, which allows to arrange the sites in a square lattice. A finite virtual dimension still introduces a finite correlation length, but due to the topology of the PEPS, this is sufficient for two-dimensional systems, even in the thermodynamic limit. Analogous extensions exist for other lattice topologies.

The ansatz for critical one-dimensional systems is called the multi-scale entanglement renormalization ansatz (MERA) Vidal (2007), see Fig. 3. This ansatz has two axes: along the physical one-dimensional lattice and along the renormalization direction. Consider two sites separated by along . The number of virtual bonds between those sites is only of order . With finite , all correlation functions measured in a MERA decay exponentially with increasing renormalization distance :

(29)

and therefore only algebraically with increasing lattice distance Vidal (2007); Evenbly and Vidal (2011).

An inconvenient property of the PEPS, MERA, and MPS with periodic boundary conditions Verstraete et al. (2004), is the introduction of loops in the network. This results in the inability to exploit the TNS gauge invariance to work with orthonormal renormalized environment states, see sections 4.2 and 4.3. One particular network which avoids such loops, but which is still able to capture polynomially decaying correlation functions, is the tree TNS (TTNS) Shi et al. (2006); Ferris (2013), see Fig. 3. From a central tensor with virtual bonds, consecutive onion-like layers are built of tensors with also virtual bonds. The last layer consists of tensors with only 1 virtual bond. An MPS is hence a TTNS with . The number of sites increases as Murg et al. (2010, 2014):

(30)

Hence for . The maximum number of virtual bonds between any two sites is . The correlation functions in a TTNS with finite and decrease exponentially with increasing separation :

(31)

and therefore only algebraically with increasing number of sites Shi et al. (2006); Ferris (2013).

For higher-dimensional or critical systems, DMRG can still be useful Stoudenmire and White (2012). The virtual dimension then has to be increased to a rather large size to obtain numerical convergence. In the case of multiple dimensions, the question arises if one should work in real or momentum space, and how the corresponding single-particle degrees of freedom should be mapped to the one-dimensional lattice Xiang (1996). Ab initio quantum chemistry can be considered as a higher-dimensional system, due to the full-rank two-body interaction in the Hamiltonian (8), and the often compact spatial extent of molecules. Nevertheless, DMRG turned out to be very useful for ab initio quantum chemistry (QC-DMRG) White and Martin (1999); Daul et al. (2000); Mitrushenkov et al. (2001); Chan and Head-Gordon (2002); Legeza et al. (2003a); Chan and Head-Gordon (2003); Legeza et al. (2003b); Mitrushenkov et al. (2003); Legeza and Sólyom (2003); Chan (2004); Chan et al. (2004); Legeza and Sólyom (2004); Moritz et al. (2005a); Chan and Van Voorhis (2005); Moritz et al. (2005b); Moritz and Reiher (2006); Hachmann et al. (2006); Rissler et al. (2006); Moritz and Reiher (2007); Dorando et al. (2007); Hachmann et al. (2007); Marti et al. (2008); Zgid and Nooijen (2008a, b, c); Ghosh et al. (2008); Chan (2008); Yanai et al. (2009); Dorando et al. (2009); Kurashige and Yanai (2009); Yanai et al. (2010); Neuscamman et al. (2010); Marti and Reiher (2010a); Luo et al. (2010); Mizukami et al. (2010); Marti et al. (2010); Murg et al. (2010); Marti and Reiher (2011); Barcza et al. (2011); Boguslawski et al. (2011); Kurashige and Yanai (2011); Mitrushchenkov et al. (2012); Sharma and Chan (2012a); Wouters et al. (2012); Boguslawski et al. (2012a); Yanai et al. (2012); Boguslawski et al. (2012b); Mizukami et al. (2013); Nakatani and Chan (2013); Boguslawski et al. (2013); Kurashige et al. (2013); Ma and Ma (2013); Saitow et al. (2013); Liu et al. (2013); Tecmer et al. (2014); Nakatani et al. (2014); Knecht et al. (2014); Wouters et al. (2014); Harris et al. (2014); Mottet et al. (2014); Lan et al. (2014); Murg et al. (2014); Sharma et al. (2014); Kurashige et al. (2014); Wouters et al. (2014); Fertitta et al. (2014); Chan et al. (2008); Chan and Zgid (2009); Marti and Reiher (2010b); Chan and Sharma (2011); Chan (2012); Kurashige (2014); Keller and Reiher (2014).

An excellent description of QC-DMRG in terms of renormalization transformations is given in Chan and Head-Gordon Chan and Head-Gordon (2002). Section 4 contains a description in terms of the underlying MPS ansatz, because this approach will be used in section 9 to introduce symmetry in the DMRG algorithm. The properties of the DMRG algorithm are discussed in section 5. Several convergence strategies are listed in section 6. An overview of the strategies to choose and order orbitals is given in section 7. A converged DMRG calculation can be the starting point of other methods. These methods are summarized in section 8. Section 10 gives an overview of the currently existing QC-DMRG codes, and the systems which have been studied with them.

4 The QC-DMRG algorithm

4.1 The MPS ansatz

DMRG can be formulated as the variational optimization of an MPS ansatz Östlund and Rommer (1995); Rommer and Östlund (1997). The MPS ansatz (27) has open boundary conditions, because sites 1 and L only have one virtual index. The sites are assumed to be orbitals, which have 4 possible occupancies , , , and . Henceforth will be used as a shorthand for . To be of practical use, the virtual dimensions are truncated to : . With increasing , the MPS ansatz spans a larger region of the full Hilbert space, but it is of course not useful to make larger than as the MPS ansatz then spans the whole Hilbert space.

A Slater determinant has gauge freedom: a rotation in the occupied orbital space alone, or a rotation in the virtual orbital space alone, does not change the physical wavefunction. Only occupied-virtual rotations change the wavefunction. An MPS has gauge freedom as well. If for two neighbouring sites and , the left MPS tensors are right-multiplied with the non-singular matrix :

(32)

and the right MPS tensors are left-multiplied with the inverse of :

(33)

the wavefunction does not change, i.e. :

(34)

4.2 Canonical forms

The two-site DMRG algorithm consists of consecutive sweeps or macro-iterations, where at each sweep step the MPS tensors of two neighbouring sites are optimized in the micro-iteration. Suppose these sites are and . The gauge freedom of the MPS is used to bring it in a particular canonical form. For all sites to the left of , the MPS tensors are left-normalized:

(35)

and for all sites to the right of , the MPS tensors are right-normalized:

(36)

Left-normalization can be performed with consecutive QR-decompositions:

(37)

The MPS tensor is now left-normalized. The -matrix is multiplied into . From site 1 to , the MPS tensors are left-normalized this way, without changing the wavefunction. Right-normalization occurs analogously with LQ-decompositions. In section 4.4, it will become clear that this normalization procedure only needs to occur at the start of the DMRG algorithm.

At this point, it is instructive to make the analogy to the renormalization group formulation of the DMRG algorithm. Define the following vectors:

(38)
(39)

Due to the left- and right-normalization described above, these vectors are orthonormal:

(40)
(41)

and are renormalized bases of the many-body Hilbert spaces spanned by resp. orbitals 1 to and orbitals to . Consider for example the left side. For site from 1 to , the many-body basis is augmented by one orbital and subsequently truncated again to at most renormalized basis states:

(42)

DMRG is hence a renormalization group for increasing many-body Hilbert spaces. The next section addresses how this renormalization transformation is chosen.

4.3 Micro-iterations

Combine the MPS tensors of the two sites under consideration into a single two-site tensor:

(43)

At the current micro-iteration of the DMRG algorithm, (the flattened column form of the tensor ) is used as an initial guess for the effective Hamiltonian equation. This equation is obtained by variation of the Lagrangian Chan (2008)

(44)

with respect to the complex conjugate of :

(45)

The canonical form in Eqs. (35)-(36) ensured that no overlap matrix is present in this effective Hamiltonian equation. In the DMRG language, this equation can be interpreted as the approximate diagonalization of the exact Hamiltonian in the orthonormal basis , see Fig. 4. Because of the underlying MPS ansatz, DMRG is variational: is always an upper bound to the energy of the true ground state.

Figure 4: Optimization of the MPS tensors at sites and in the two-site DMRG algorithm. The effective Hamiltonian equation (45), obtained by variation of the Lagrangian (44), can be interpreted as the approximate diagonalization of the exact Hamiltonian in the orthonormal basis .

The lowest eigenvalue and corresponding eigenvector of the effective Hamiltonian are searched with iterative sparse eigensolvers. Typical choices are the Lanczos or Davidson algorithms Lanczos (1950); Davidson (1975). Once is found, it is decomposed with an SVD:

(46)

Note that is hence left-normalized and right-normalized. The sum over is truncated if there are more than nonzero Schmidt values , thereby keeping the largest ones. This is the optimal approximation for the bipartition of into and . In the original DMRG algorithm, and were obtained as the eigenvectors of resp.  and .

A discarded weight can be associated with the truncation of the sum over :

(47)

This is the probability to measure one of the discarded states in the subsystems or . The approximation introduced by the truncation becomes better with increasingly small discarded weight. Instead of working with a fixed , one could also choose dynamically in order to keep below a preset threshold, as is done in Legeza’s dynamic block state selection approach Legeza et al. (2003a).

4.4 Macro-iterations or sweeps

So far, we have looked at a micro-iteration of the DMRG algorithm. This micro-iteration happens during left or right sweeps. During a left sweep, is constructed, the corresponding effective Hamiltonian equation solved, the solution decomposed, the Schmidt spectrum truncated, is contracted into , is set to this contraction , is set to , and is decreased by 1. Note that is right-normalized for the next micro-iteration as required. This stepping to the left occurs until , and then the sweep direction is reversed from left to right. Based on energy differences, or wavefunction overlaps, between consecutive sweeps, a convergence criterium is triggered, and the sweeping stops.

DMRG can be regarded as a self-consistent field method: at convergence the neighbours of an MPS tensor generate the field which yields the local solution, and this local solution generates the field for its neighbours Chan and Head-Gordon (2002); Hachmann et al. (2006); Chan (2008).

4.5 Renormalized operators and their complements

The effective Hamiltonian in Eq. (45) is too large to be fully constructed as a matrix. Only its action on a particular guess is available as a function. In order to construct efficiently for general quantum chemistry Hamiltonians, several tricks are used. Suppose that a right sweep is performed and that the MPS tensors of sites and are about to be optimized.

Renormalized operators such as with are constructed and stored on disk White and Martin (1999); Chan and Head-Gordon (2002); Kurashige and Yanai (2009). The renormalized operators needed for the previous micro-iteration can be recycled to this end. Suppose :

(48)

Note that no phases appear because an even number of second-quantized operators was transformed. For an odd number, there should be an additional phase at the right-hand side (RHS) due to the Jordan-Wigner transformation Jordan and Wigner (1928). Renormalized operators to the right of can be loaded from disk, as they have been saved during the previous left sweep.

Once three second-quantized operators are on one side of , they are multiplied with the matrix elements , and a summation is performed over the common indices to construct complementary renormalized operators Xiang (1996); White and Martin (1999); Chan and Head-Gordon (2002); Kurashige and Yanai (2009):

(49)

For two, three, and four second-quantized operators on one side of , these complementary renormalized operators are constructed. A bare renormalized operator (without matrix elements) is only constructed for one or two second-quantized operators.

Hermitian conjugation and commutation relations:

(50)

are also used to further limit the storage requirements for the (complementary) renormalized operators. Examples of renormalized operators and the fermion sign handling can be found in, for example, Refs. Wouters et al. (2012); Wouters (2014a).

4.6 Computational cost

This section describes the cost of the QC-DMRG algorithm per sweep in terms of memory, disk, and computational time White and Martin (1999); Chan and Head-Gordon (2002); Kurashige and Yanai (2009). To analyze this cost, let us first look at the cost per micro-iteration. A micro-iteration consists of three steps: solving the effective Hamiltonian equation (45), performing an SVD of the solution (46), and constructing the (complementary) renormalized operators for the next micro-iteration.

To solve the effective Hamiltonian equation with the Lanczos or Davidson algorithms, a set of trial vectors are kept in memory, as well as . To construct , (complementary) renormalized operators should also be stored in memory. The latter have at most two site indices. The total memory cost is hence .

The action of on is divided into several contributions. Each contribution consists of the joint action of a renormalized operator and the corresponding complementary renormalized operator. For each contribution, two matrix-matrix multiplications need to be performed, of computational cost . In total there are contributions, because complementary renormalized operators have at most two site indices. The total computational cost is hence for the multiplications, and for the summation of the different contributions.

The SVD of the solution takes computational time and memory.

The construction of one particular renormalized operator takes computational time and memory, and there are such operators. The most tedious part to analyze is the construction of the two-site complementary renormalized operators, e.g.

(51)

which takes at first sight computational time and memory per operator. There are such operators, and a naive implementation would hence result in a computational cost of per micro-iteration. However, this summation needs to be performed only once for each operator, at the moment when the second second-quantized operator is added:

(52)

From then on, this operator can be transformed as in Eq. (48). The total computational cost per micro-iteration is hence reduced to for the summation (there are three variable site indices in Eq. (52)), and for the transformation (there are operators to be transformed). The one-site complementary renormalized operator (the complement of three second-quantized operators) can be constructed from the two-site complementary renormalized operators at the moment when the third second-quantized operator is added. From then on, this operator can also be transformed as in Eq. (48).

As mentioned earlier, the (complementary) renormalized operators are stored to disk, as well as the MPS site tensors, in order to be recycled when the sweep direction is reversed. An overview of the resulting total cost per macro-iteration is given in Tab. 1. For a given virtual dimension , the DMRG algorithm is of polynomial cost in . The computational requirements in Tab. 1 are upper bounds if the symmetry group of the Hamiltonian is exploited, see section 9. Then the MPS tensors and corresponding (complementary) renormalized operators become block-sparse, and is not full rank. An example of the scaling of the computational time per DMRG sweep with the number of orbitals is shown in Fig. 5. Due to the imposed symmetry, CheMPS2 Wouters et al. (2014) achieves a scaling below .

time memory disk
-
SVD and basis truncation
Renormalized operators
Complementary renormalized operators
Total

The memory for the (complementary) renormalized operators is mentioned separately.

Table 1: Computational requirements per macro-iteration or sweep of the QC-DMRG algorithm.
Figure 5: The geometries of all-trans polyenes CH were optimized at the B3LYP/6-31G** level of theory for , 14, 16, 18, 20, 22 and 24. The -orbitals were kept frozen at the RHF/6-31G level of theory. The -orbitals in the 6-31G basis were localized by means of the Edmiston-Ruedenberg localization procedure Edmiston and Ruedenberg (1963), which maximizes . The localized -orbitals belong to the irrep of the point group, and were ordered according to the one-dimensional topology of the polyene. For all polyenes, the average CPU time per DMRG sweep was determined with CheMPS2 Wouters et al. (2014), for two reduced virtual dimensions . For the values of shown here, the energies are converged to accuracy due to the one-dimensional topology of the localized and ordered -orbitals. Due to the imposed symmetry, all tensors become block-sparse, see section 9, which causes the scaling to be below .

5 Properties

5.1 DMRG is variational

The DMRG algorithm is variational, because it can be formulated as the optimization of an MPS ansatz. All energies obtained during all micro-iterations are therefore upper bounds to the true ground state energy. These energies do not go down monotonically however, because the basis in which is diagonalized changes between different micro-iterations due to the truncation of the Schmidt spectrum Chan and Head-Gordon (2002).

5.2 Energy extrapolation

With increasing virtual dimension , the MPS ansatz spans an increasing part of the many-body Hilbert space. In the following, denotes the minimum energy encountered in Eq. (45) during the micro-iterations for a given virtual dimension . Several calculations with increasing can be performed, in order to assess the convergence. This even allows to make an extrapolation of the energy to the FCI limit. Several extrapolation schemes have been suggested. Note that and below are parameters to be fitted. The maximum discarded weight encountered during the last sweep before convergence is abbreviated as:

(53)

The initial assumption of exponential convergence White and Martin (1999)

(54)

was rapidly abandoned for the relation Legeza and Fáth (1996); Chan and Head-Gordon (2002); Legeza et al. (2003a)

(55)

because the energy is a linear function of the RDM Chan and Head-Gordon (2002). An example of an extrapolation with Eq. (55) is shown in Fig. 6. The tail of the distribution of RDM eigenvalues scales as Chan et al. (2002); Chan and Head-Gordon (2002)

(56)

Substituting this relation in Eq. (55) yields an improved version of Eq. (54) Chan and Head-Gordon (2002):

(57)

An example of an extrapolation with Eq. (57) is shown in Fig. 10. Eqs. (55) and (57) are the most widely used extrapolation schemes in QC-DMRG. Three other relations have been proposed as well. A relation for incremental energies has been suggested Mitrushenkov et al. (2003):

(58)

but the extrapolated often violates the variational principle. An alternative relation based on the discarded weight has also been proposed Mitrushenkov et al. (2003):

(59)

as well as a Richardson-type extrapolation scheme, based on the assumption that the energy is an analytic function of Marti and Reiher (2010a):

(60)
Figure 6: Extrapolation of the variational DMRG ground-state energy with the discarded weight , for N in the cc-pVDZ basis near equilibrium (nuclear separation 2.118 a.u.). The calculation was performed with CheMPS2 Wouters et al. (2014) with symmetry, see section 9. denotes the number of reduced virtual basis states. The irrep ordering in the DMRG calculation was [ABBBBBBA] in order to place bonding and antibonding orbitals close to each other on the one-dimensional DMRG lattice, see section 7.3 and Fig. 8.

5.3 The CI content of the wavefunction

To analyze the MPS wavefunction (27), suppose that the orthonormal orbitals are the HF orbitals. An important difference with traditional post-HF methods such as CI expansions, is that no FCI coefficients are a priori zero. An MPS hence captures CI coefficients of any particle-excitation rank relative to HF Chan et al. (2004); Hachmann et al. (2006). A small virtual dimension implies little information content in the FCI coefficient tensor, or equivalently that the many nonzero FCI coefficients are in fact highly correlated. This has to be contrasted with CI expansions, which are truncated in their particle-excitation rank and therefore set many FCI coefficients a priori to zero. The nonzero FCI coefficients are however not a priori correlated in a CI expansion: they are entirely free to be variationally optimized.

5.4 Size-consistency

For a method to be size-consistent, the compound wavefunction should be multiplicatively separable and the energy additively separable for noninteracting subsystems and . From the discussion of the Schmidt decomposition above, it follows immediately that an MPS is size-consistent if the orbitals of subsystems and do not overlap, and if they are separated into two groups on the one-dimensional DMRG lattice Chan and Head-Gordon (2002); Chan and Sharma (2011). The latter is for example realized if orbitals to correspond to subsystem and orbitals to correspond to subsystem . DMRG will then automatically retrieve a product wavefunction, in which only one Schmidt value is nonzero at the corresponding boundary.

5.5 DMRG is not FCI

A good variational energy does not necessarily imply that the wavefunction is accurate. Suppose we have an orthonormal MPS with virtual dimension which has been variationally optimized to approximate the true ground state with energy . Suppose that

(61)

with . Then

(62)

and

(63)

The energy converges quadratically in the wavefunction error. Most DMRG convergence criteria rely on energy convergence (), see Fig. 6. An important implication is that, except for tremendously large virtual dimensions where , the MPS wavefunction is not invariant to orbital rotations. The orbital choice and their ordering on a one-dimensional lattice also influence the convergence rate with . Strategies to choose and order orbitals are discussed in section 7. Sparse iterative FCI eigensolvers converge the FCI tensor to a predefined threshold instead of the energy. An FCI solution can therefore be considered invariant to orbital rotations.

6 Convergence strategies

The DMRG algorithm can get stuck in a local minimum or a limit cycle, if is insufficiently large Chan and Head-Gordon (2002). The chance of occurrence is larger for inconvenient orbital choices and orderings. Because the virtual dimension cannot be increased indefinitely in practice, it is important to choose the set of orbitals and their ordering well, see section 7. Additional considerations to enhance convergence are described here.

6.1 The number of sites to be optimized in a micro-iteration

It is better to use the two-site DMRG algorithm than the one-site version White (2005). In the one-site version, the Hamiltonian is diagonalized during the micro-iterations in the basis instead of . Because of the larger variational freedom in the two-site DMRG algorithm, lower energy solutions are obtained, and the algorithm is less likely to get stuck Zgid and Nooijen (2008b). It might therefore be worthwhile to optimize three or more MPS tensors simultaneously in a micro-iteration, or to group several orbitals into a single DMRG lattice site White and Martin (1999).

The two-site algorithm has another important advantage, when the symmetry group of the Hamiltonian is exploited. The virtual dimension is then distributed over several symmetry sectors, see section 9. In the one-site algorithm, the virtual dimension of a symmetry sector has to be changed manually during the sweeps Zgid and Nooijen (2008b), while the SVD (46) in the two-site algorithm automatically picks the best distribution.

6.2 Perturbative corrections and noise

White suggested to add perturbative corrections to the RDM in order to enhance convergence White (2005). Instead of using perturbative corrections, one can also add noise to the RDM prior to diagonalization or to prior to SVD Chan and Head-Gordon (2002). The corrections or noise help to reintroduce lost symmetry sectors (lost quantum numbers) in the renormalized basis, which are important for the true ground state. Instead of adding noise or perturbative corrections, one can also reserve a certain percentage of the virtual dimension to be distributed equally over all symmetry sectors Chan (2004).

6.3 Getting started

The wavefunction from which the QC-DMRG algorithm starts has an influence on the converged energy (by getting stuck in a local minimum) and on the rate of convergence Legeza et al. (2003a); Legeza and Sólyom (2003); Moritz and Reiher (2006). The effect of the starting guess is estimated to be an order of magnitude smaller than the effect of the choice and ordering of the orbitals Moritz and Reiher (2006). Nevertheless, it deserves attention.

One possibility is to choose a small active space to start from, and subsequently augment this active space stepwise with previously frozen orbitals Mitrushenkov et al. (2001), in analogy to the infinite-system DMRG algorithm White (1992). Natural orbitals from a small CASSCF calculation or HF orbitals can be used to this end Moritz and Reiher (2006). An alternative is to make an a priori guess of how correlated the orbitals are. This can be done with a DMRG calculation with small virtual dimension , from which the approximate single-orbital entropies can be obtained. The subsystem is then chosen to be a single orbital in Eq. (25). The larger the single-orbital entropy, the more it is correlated. The active space can then be chosen and dynamically extended based on the single-orbital entropies Barcza et al. (2011).

One can also decompose the wavefunction from a cheap CI calculation with single and double excitations into an MPS to start from Chan and Head-Gordon (2002); Moritz and Reiher (2006). Another possibility is to distribute equally over the symmetry sectors, and to fill the MPS with noise. This retrieves energies below the HF energy well within the first macro-iteration Wouters et al. (2014); Wouters (2014a).

To achieve a very accurate MPS quickly, it is also best to start from calculations with relatively small virtual dimension , and to enlarge it stepwise Chan and Head-Gordon (2002); Moritz and Reiher (2006); White (1996).

7 Orbital choice and ordering

There are many ways to set up a renormalization group flow, and the specific setup influences the outcome. One consideration of key importance in QC-DMRG is the choice and ordering of orbitals. Most molecules or active spaces are far from one-dimensional. By placing the orbitals on a one-dimensional lattice, and by assuming an MPS ansatz with modest , an artifical correlation length is introduced in the system, which can be a bad approximation. Over time, several rules of thumb have been established to choose and order the orbitals.

7.1 Elongated molecules

Quantum information theory learns that locality is an important concept, see section 3. The Coulomb interaction is however long-ranged. On the other hand, the mutual screening of electrons and nuclei can result in an effectively local interaction. For elongated molecules such as hydrogen chains Chan and Head-Gordon (2002); Hachmann et al. (2006); Mitrushchenkov et al. (2012); Wouters et al. (2012); Nakatani and Chan (2013); Ma and Ma (2013), polyenes Chan and Head-Gordon (2002); Chan and Van Voorhis (2005); Hachmann et al. (2006); Ghosh et al. (2008); Yanai et al. (2010), or acenes Hachmann et al. (2007); Dorando et al. (2007); Mizukami et al. (2013), which are more or less one-dimensional, choosing a spatially local basis has turned out to be very beneficial. There are roughly three ways to choose a local basis: symmetric orthogonalization as it lies closest to the original gaussian basis functions Hachmann et al. (2007); Dorando et al. (2007); Mitrushchenkov et al. (2012); Wouters et al. (2012); Ma and Ma (2013); Carlson and Keller (1957), explicit localization procedures such as Pipek-Mezey or Edmiston-Ruedenberg Ghosh et al. (2008); Mizukami et al. (2013); Pipek and Mezey (1989); Edmiston and Ruedenberg (1963), and working in a biorthogonal basis Chan and Van Voorhis (2005); Mitrushchenkov et al. (2012). For the latter, the effective Hamiltonian is not hermitian anymore. The DMRG algorithm should then be correspondingly adapted Mitrushenkov et al. (2003); Chan and Van Voorhis (2005); Mitrushchenkov et al. (2012). The adapted algorithm is slower and prone to convergence issues, and it is therefore better to use one of the other two localized bases Chan and Van Voorhis (2005); Mitrushchenkov et al. (2012). Fig. 7 illustrates the speed-up in energy convergence by using a localized basis for all-trans polyenes.

Figure 7: The computational details were discussed in the caption of Fig. 5. The active space of CH, which consists of 28 -orbitals, is studied both with ordered localized orbitals (Edmiston-Ruedenberg) and canonical orbitals (restricted HF). The energy converges significantly faster with the number of reduced virtual basis states when ordered localized orbitals are used.

7.2 Hamiltonian measures

If the topology of the molecule does not provide hints for choosing and ordering orbitals, it was investigated whether the Hamiltonian (1) can be of use. Several integral measures have been proposed, for which a minimal bandwidth is believed to yield a good orbital order. Chan and Head-Gordon proposed to minimize the bandwidth of the one-electron integral matrix of the HF orbitals Chan and Head-Gordon (2002). In quantum chemistry, it is often stated that the one-electron integrals are an order of magnitude larger than the two-electron integrals, and that quantum chemistry therefore corresponds to the small- limit of the Hubbard model Legeza et al. (2003a); Barcza et al. (2011); Hubbard (1963). On the other hand, there are many two-electron integrals, and they may become important due to their number. When other orbitals than the HF orbitals are used, it may therefore be interesting to minimize the bandwidth of the Fock matrix Legeza et al. (2003b):

(64)

Other proposed integral measures are the MP2-inspired matrix Mitrushenkov et al. (2003):

(65)

where are the HF single-particle energies, as well as several measures in Ref. Moritz et al. (2005a). These are the Coulomb matrix , the exchange matrix , the mean-field matrix , and two derived quantities:

(66)
(67)

While the one-electron integrals vanish when orbitals and belong to different molecular point group irreps, and do not. Ref. Moritz et al. (2005a) used a genetic algorithm to find the optimal HF orbital ordering, in order to assess the proposed integral measures. This genetic algorithm was expensive, which limited its usage to small test systems. It favoured bandwidth minimization, although no de