We present an efficient stochastic algorithm for the recently introduced perturbative density matrix renormalization group method for large active spaces. The stochastic implementation bypasses the computational bottleneck involved in solving the first order equation in the earlier deterministic algorithm. We demonstrate the efficiency and accuracy of the algorithm on the C2 and Cr2 molecular benchmark systems.

Molecular electronic structure with a mix of static correlation and dynamic correlation remains a challenging aspect of quantum chemistry. It is often treated within a complete active space (CAS) approach for the static correlation plus a second-order perturbative treatment of the dynamic correlation, such as with CASPT21,2 (complete-active-space second-order perturbation theory) and NEVPT23–5 (N-electron valence state perturbation theory). However, in many cases, the second-order treatment of correlations outside of the active space is not quantitatively accurate, and it becomes necessary to enlarge the active space with orbitals of intermediate correlation strength. For instance, in 3d transition metal complexes, the virtual 4d, semi-core (3s and 3p), or valence ligand orbitals often need to be included in the active space for good accuracy.6–8 The resulting very large number of active orbitals with a mix of different correlation character renders the density matrix renormalization group (DMRG),9–19 whose strength lies in treating large active spaces where all the orbitals are strongly interacting, relatively inefficient, and the subsequent dynamic correlation treatment also becomes very difficult.7,8,20–25

Recently, we developed a perturbative DMRG (p-DMRG)26 algorithm to efficiently target near-exact energies in large orbital spaces where there is a mix of strongly and intermediately correlated orbitals. We showed that the p-DMRG provides, in practice, benchmark quality energies as accurate as those obtained by variational DMRG calculations, but with a significantly reduced computational cost. The essential idea is to treat the strong correlation with a variational DMRG calculation with small bond dimension M0 (which can be systematically increased) within the full orbital space, and the residual correlation by second-order perturbation theory (PT2) within the same orbital space. The total energy obtained by p-DMRG converges very quickly with M0, which can thus be chosen significantly smaller than the bond dimension M in the fully variational calculation. The basic idea of p-DMRG is similar in spirit to that of selected configuration interaction (SCI) plus perturbation theory,27–35 and in particular, the recently developed Heat-bath CI+PT (HCI+PT),33,34 where the static correlation is treated by configuration interaction within a selected set of important determinants defined by a threshold, which can be tuned to systematically increase the accuracy much like M0 in p-DMRG, with the residual dynamic correlation in the same orbital space treated by standard PT2.

In this communication, we describe an efficient stochastic algorithm that further speeds up the p-DMRG algorithm by avoiding having to solve the first order equation deterministically. Such a stochastic extension to HCI+PT has already been described.34 However, in the selected CI setting, the stochastic extension is very natural due to a clear separation of the variational and first order determinant spaces. Since in p-DMRG, the zeroth order wavefunction is a matrix product state (MPS),19 the stochastic extension is less straightforward and this communication is thus concerned primarily with its correct and efficient formulation.

To begin with, we assume that a zeroth order wavefunction |Ψ(0)⟩, expressed as an MPS with a small bond dimension M0, has been optimized by the standard variational DMRG algorithm. Given a partitioning of the Hamiltonian, Ĥ=Ĥ0+V^ with Ĥ0(0)⟩ = E0(0)⟩, the first order wavefunction |Ψ(1)⟩ is defined by the first order equation,

Q(Ĥ0E0)Q|Ψ(1)=QV^|Ψ(0),
(1)

where Q = 1 − P and P = |Ψ(0)⟩⟨Ψ(0)| are projectors. Once |Ψ(1)⟩ is obtained from Eq. (1), which in our implementation is achieved by minimizing the Hylleraas functional using the MPSPT algorithm,20 the second order energy E2=Ψ(1)|QV^|Ψ(0) can be computed. For very large active spaces with 50–100 orbitals, the required bond dimension M1 for expanding |Ψ(1)⟩, which typically scales like O(K2M0) (K is the number of spatial orbitals), can be quite large (≫104). This creates a significant cost both in computation and in storage. In the previous deterministic p-DMRG algorithm,26 we used a sum of MPS representation |Ψ(1)=i|Ψi(1) as well as extrapolation, which helped alleviate the computational cost. By contrast, a stochastic algorithm can be expected to essentially eliminate this bottleneck, at the cost of introducing statistical errors.

To develop a stochastic variant of p-DMRG, we first rewrite E2 as

E2=Ψ(0)|V^Q[Q(Ĥ0E0)Q]1QV^|Ψ(0)
(2)

and then aim to find an explicit expression for E2 as a sum over terms, ideally of a similar form to i|Di|V^|Ψ(0)|2EiE0 as appears in HCI+PT with |Di⟩ being a determinant, which can then be sampled stochastically. In our previous paper,26 we chose Ĥ0 as

Ĥ0=PE0P+QĤdQ,
(3)

where Ĥd contains those operators in Ĥ which do not change the occupation numbers of spatial orbitals, i.e.,

Ĥd=phppÊpp+12pq(pp|qq)êpqqp+12pq(pq|qp)êpqpq,
(4)

with Êpq=σapσaqσ and êpqrs=σ,τapσaqτarτasσ=EpsEqrδqsEpr. This Hamiltonian is spin-free; however, it is only block-diagonal in the determinant space since it contains additional couplings for determinants with the same spatial occupations.26 Thus, in this work, we instead used the Epstein-Nesbet (EN) form, which is equivalent to neglecting the above off-diagonal couplings in the determinant space. In the following, we will use Ĥd to denote the EN form, which satisfies ⟨Di|Ĥd|Dj⟩ = δijDi|Ĥ|Dj⟩ ≜ δijEi. This choice will result in a slight difference in the computed E2 compared with our previous p-DMRG, usually found to be less than 1mEh (vide post), and the difference will gradually decrease as M0 increases.

The EN partition is, of course, commonly used in SCI+PT.27–35 In this case, Q[Q(Ĥ0E0)Q]1Q in (2) can be simplified to Q(ĤdE0)1Q for Q = i|Di⟩⟨Di| and [Q, Ĥd] = 0, giving E2=i|Di|V^|Ψ(0)|2EiE0 as a sum over determinants in the Q-space, as mentioned above. However, in our case, Q(ĤdE0)Q and hence Q[Q(Ĥ0E0)Q]1Q are not diagonal in the determinant space due to the more complicated structure of Q arising from the fact that the zeroth order MPS potentially can have nonzero overlap with every determinant. To derive a similar expression for E2, we first rewrite

Q(Ĥ0E0)Q=Q(ĤdE0)QXY,X=(ĤdE0),Y=PX+XPPXP.
(5)

In general, the operator X = ĤdE0 can be assumed to be invertible by properly choosing E0, and this important issue will be discussed later. Then, using the relation Q(XY)1Q=QX1(1YX1)1Q=QX1Q+QX1YX1Q+QX1YX1YX1Q+ and directly evaluating each term, we can find that the series can be greatly simplified as P is one-dimensional such that Q(XY)−1Q becomes

Q(Q(Ĥ0E0)Q)1Q=QX1QQX1PX1QΨ(0)|X1|Ψ(0),
(6)

which can be simply verified by multiplying the right-hand side with Q(Ĥ0E0)Q to get the identity operator Q in the Q-space. The last term is the correction for the fact that Q(ĤdE0)1P is in general nonzero, which clearly vanishes in the SCI+PT case as [Q, Ĥd] = 0. The result is also easily derived by using Löwdin partitioning for X = (ĤdE0) and blockwise inversion to fold in the self-energy contribution of the single P space Hamiltonian element on the Q space denominators. Substituting Eq. (6) into Eq. (2) for E2 and invoking the resolution of identity i|Di⟩⟨Di| = 1, each term in E2 can finally be formulated as a sum over determinants, viz.,

E2A+|B|2C,
(7)
AΨ(0)|V^QX1QV^|Ψ(0)=i|Ψ(0)|V^Q|Di|2EiE0,
(8)
BΨ(0)|V^QX1|Ψ(0)=iΨ(0)|V^Q|DiDi|Ψ(0)EiE0,
(9)
CΨ(0)|X1|Ψ(0)=i|Ψ(0)|Di|2EiE0.
(10)

These formulae constitute the final working equations for stochastic p-DMRG. Note that unlike in the case of SCI+PT, here the summation over |Di⟩ goes over all the determinants.

In the practical evaluation of A, B, and C, instead of sampling |Di⟩ uniformly, we use importance sampling to improve the efficiency. For simplicity, we first discuss the evaluation of the term C(10). It can be evaluated as C=1EiE0Pi=|Ψ(0)|Di|2, where the subscript indicates that the average of 1EiE0 is taken with respect to the population of the determinants generated with the probability Pi = |⟨Ψ(0)|Di⟩|2. To achieve this, we use an algorithm similar to that in Refs. 36 and 37. Specifically, suppose |Ψ(0)⟩ is in the right canonical form,

|Ψ(0)=n1n2nKCn1[1]Rn2[2]RnK[K]|n1n2nK,
(11)

where Cα1n1[1] and RαK1nK[K] are matrices and Rαk1αknk[k] are rank-3 tensors satisfying the right canonical condition nkαkRαk1αknk[k]Rαk1αknk[k]=δαk1αk1, then a determinant |Di⟩ = |n1n2nK⟩ can be sampled according to Pi by a single sweep from left to right as follows: The first occupation number n1 ∈{|vac⟩, |1β⟩, |1α⟩, |1α1β⟩} is generated according to p1(n1)α1|Cα1n1|2, which satisfies n1p1=1 as |Ψ(0)⟩ is normalized. Given n1, at the second site, the generation probability for n2 is defined as p2(n2|n1)α2|(Cn1Rn2)α2|2/N2. Importantly, the normalization constant N2 can be found as N2=n2p2(n2|n1)=p1(n1) due to the right canonical property. Repeating this procedure recursively, the generation probability for nk given n1nk−1 can be defined as pk(nk|n1nk1)α1|(Cn1Rn2Rnk)αk|2/Nk with Nk = pk−1(nk−1|n1nk−2). Thus, after the occupation of the last site is chosen, the total generation probability is P(n1nK)=p1(n1)p2(n2|n1)p(nk|n1nK1)=|(Cn1Rn2RnK)|2=|Ψ(0)|Di|2=Pi, equal to our target distribution. For a spin-adapted DMRG implementation,15 the determinants can be generated similarly, but at each step the Clebsch-Gordon coefficient CSk1Mk1skmkSkMk needs to be incorporated to map each reduced tensor to the full one.

For term A(8), it would be possible to employ the same strategy, if QV^|Ψ(0) can be expressed as an MPS faithfully. However, since converting QV^|Ψ(0) into an MPS exactly would incur a bond dimension of O(K2M0), this becomes prohibitive for large active spaces. Thus, we first use a bond dimension M1 [that is small compared with O(K2M0)] to compress QV^|Ψ(0) variationally as an MPS, i.e., |ΦQV^|Ψ(0). This approximate state |Φ⟩ is then used to define a generation probability for |Di⟩, Pi = |⟨Φ|Di⟩|2, for importance sampling, such that

A=|Ψ(0)|V^Q|Di|2Pi(EiE0)Pi=|Φ|Di|2.
(12)

In the case that |Φ⟩ is a good approximation to QV^|Ψ(0), an approximate estimator for A is just

A1EiE0Pi=|Φ|Di|2,
(13)

which becomes similar to the expression for term C but with respect to a different distribution. This approximation becomes exact in the limiting case that |Φ=QV^|Ψ(0). It deserves to be noted that in Eq. (12), there is a subtlety. For the equality to hold, the set of determinants S{|Di:Φ|Di0} interacting with |Ψ(0)⟩ must include the set that interacts with Ψ(0)|V^Q|Di0, otherwise some contributions to E2 will be missing. This is guaranteed for sufficiently large M1, and fortunately, we found that to converge to chemical accuracy ca. 1mEh, the required M1 can be much smaller than O(K2M0) (vide post). For the term B(9), we simply evaluate it as B=Ψ(0)|V^Q|DiDi|Ψ(0)Pi(EiE0)Pi=|Φ|Di|2 since it turns out to be less important in most cases due to its small size and the quadratic dependence of E2 on it.

The above algorithm constitutes a stochastic algorithm for p-DMRG. Its computational cost depends on two parts: the compression for QV^|Ψ(0) on the right-hand side of the first order equation (1), which scales as O(K3M12M0) assuming M1M0, and the cost for the stochastic evaluation of E2. The latter is dominated by the number of samples Ns times the cost for evaluating the matrix elements Ψ(0)|V^Q|Di. We found the following scheme to be efficient at least for the systems investigated in this work. By using the identity Ψ(0)|V^Q|Di=Ψ(0)|ĤQ|Di=Ψ(0)|(HEDMRG(0))|Di=jΨ(0)|DjDj|(HEDMRG(0))|Di, the evaluation is converted into O(K2N2) evaluations of overlaps between determinants and the zeroth order wavefunction ⟨Ψ(0)|Dj⟩, each of which scales as O(KM02). Thus, the total cost for the stochastic evaluation is O(NsN2K3M02), which is formally higher than O(NsK3M02) if the matrix element Ψ(0)|(HEDMRG(0))|Di is directly computed by considering |Di⟩ as an MPS with bond dimension one. However, the efficiency of our choice in practice may stem from the fact that in computing ⟨Ψ(0)|Dj⟩, the sparsity in the MPS tensors can be utilized such that the actual computational time is less than the direct evaluation of Ψ(0)|(HEDMRG(0))|Di. A detailed comparison of these two choices for large orbital spaces will be presented in future. The overall time of the stochastic step is usually found to be much smaller than the compression step. Thus, the present stochastic algorithm is more efficient than the previous p-DMRG algorithm, which required the iterative solution of the first order equation Eq. (1), which scales as O(K2M13+K3M12M0), and with a much larger M1 than required in the present algorithm.

We now demonstrate the accuracy and efficiency of the stochastic algorithm in comparison with the previous deterministic p-DMRG algorithm for two prototypical molecules: C2 at its equilibrium bond length 1.242 53 Å38 in the cc-pVDZ basis set39 and Cr2 at its equilibrium bond length 1.68 Å40 in the Ahlrichs’ SV basis41 and the cc-pVDZ-DK basis.42 For C2 and Cr2 in the SV basis, all electrons were correlated, resulting in orbital spaces with 12 electrons in 28 orbitals, (12e, 28o) for C2 and (48e,42o) for Cr2, respectively. For Cr2 in the cc-pVDZ-DK basis, the same orbital space as in our previous work26 was employed, viz., (28e, 76o) with the 1s, 2s, and 2p frozen. As in our previous work on deterministic p-DMRG,26 the zeroth order energy is defined by an interpolation of two limits,

E0(λ)=(1λ)EDMRG(0)+λEd(0),EDMRG(0)=Ψ(0)|Ĥ|Ψ(0),Ed(0)=Ψ(0)|Ĥd|Ψ(0).
(14)

For simplicity, in this work, we only explored two cases, λ = 0 and λ=12, where in the former case E0(λ=0)=EDMRG(0), and in the latter case E0(λ=12) is an average (EDMRG(0)+Ed(0))/2, which in practice we found to deliver significantly better energies for challenging systems due to the more balanced treatment of the zeroth order state and the determinant-like perturbers.26 The limit λ = 1 was previously found to be prone to the intruder state problem and not considered here.

The results for C2 obtained with stochastic p-DMRG with Ns = 36 000 are listed in Table I together with the p-DMRG results. (Note that we report Ns as the total number of samples across all computational cores in the calculation, not the number of samples per core as is sometimes reported.) The orbitals were obtained by a complete active space self-consistent field (CASSCF) calculation with an active space CAS(6e,6o). It is seen that both p-DMRG and its stochastic variant improve on the variational DMRG results at each M0, and converge toward the exact value computed by variational DMRG with M = 4000 very quickly. The λ=12 energies are better than those with λ = 0 due to the reasons discussed earlier. The bond dimension M1 used for compression is chosen as 2000, which is sufficient for Eq. (12) to hold. The total energies of stochastic p-DMRG are seen to be very close to those from its deterministic counterpart for each λ. Neglecting off-diagonal couplings in Eq. (4), which is equivalent to employing the EN partition, in the stochastic variant, leads to slightly lower energies. This can be rationalized by the fact that after diagonalization of dQ these couplings will make the singlet perturbers have higher energies in the original p-DMRG compared with Ei for determinants, and hence the original p-DMRG will have smaller second-order energies for singlet states. It is seen that in all cases, the statistical errors are very small with a moderate number of samples. The wall time at M0 = 400 is about 1 min for stochastic p-DMRG, which is 4 times faster than the deterministic one.

TABLE I.

Total energy and bond dissociation energy De of C2 (in Eh) in the cc-pVDZ basis set (12 electrons in 28 orbitals) obtained by DMRG and stochastic p-DMRG (Ns = 36 000) with various M0. Values in parentheses are the statistical uncertainties.

p-DMRGStochastic p-DMRGa
M0DMRGλ = 0λ=12λ = 0λ=12
Total energies of C2 
50 −75.708 599 −75.729 676 −75.731 740 −75.730 32(8) −75.732 53(9) 
100 −75.724 500 −75.731 540 −75.732 029 −75.731 73(3) −75.732 27(3) 
200 −75.729 502 −75.731 845 −75.731 990 −75.731 919(8) −75.732 065(9) 
400 −75.731 380 −75.731 939 −75.731 966 −75.731 947(3) −75.731 976(3) 
4000 −75.731 960     
Bond dissociation energies of C2 
50 0.184 82 0.205 87 0.207 93 0.206 51 0.208 72 
100 0.200 69 0.207 73 0.208 22 0.207 92 0.208 46 
200 0.205 69 0.208 04 0.208 18 0.208 11 0.208 26 
400 0.207 57 0.208 13 0.208 16 0.208 14 0.208 17 
4000 0.208 15     
p-DMRGStochastic p-DMRGa
M0DMRGλ = 0λ=12λ = 0λ=12
Total energies of C2 
50 −75.708 599 −75.729 676 −75.731 740 −75.730 32(8) −75.732 53(9) 
100 −75.724 500 −75.731 540 −75.732 029 −75.731 73(3) −75.732 27(3) 
200 −75.729 502 −75.731 845 −75.731 990 −75.731 919(8) −75.732 065(9) 
400 −75.731 380 −75.731 939 −75.731 966 −75.731 947(3) −75.731 976(3) 
4000 −75.731 960     
Bond dissociation energies of C2 
50 0.184 82 0.205 87 0.207 93 0.206 51 0.208 72 
100 0.200 69 0.207 73 0.208 22 0.207 92 0.208 46 
200 0.205 69 0.208 04 0.208 18 0.208 11 0.208 26 
400 0.207 57 0.208 13 0.208 16 0.208 14 0.208 17 
4000 0.208 15     
a

Stochastic p-DMRG does not converge to p-DMRG with increasing sample sizes due to the use of different Ĥ0 (see text).

Next, the different components of E2, viz., A, B, and C in Eqs. (8)–(10), are given in Table II. The B term is found to be very small such that the |B|2/C contribution is negligible for E2 in this case. However, since the calculation of B and C is relatively cheap compared with the evaluation of A, we always include them in all our calculations.

TABLE II.

Values of different components in stochastic p-DMRG with λ = 0 for C2 (in Eh) in the cc-pVDZ basis set: E2, A, B, and C in Eqs. (7)–(10), respectively.

M0ABCE2
50 −0.021 82(8) −0.000 9(1) 2.521(5) −0.021 82(8) 
100 −0.007 23(3) 0.001 5(1) 2.407(4) −0.007 23(3) 
200 −0.002 418(8) 0.000 03(5) 2.353(4) −0.002 418(8) 
400 −0.000 567(3) −0.000 007(9) 2.343(4) −0.000 567(3) 
M0ABCE2
50 −0.021 82(8) −0.000 9(1) 2.521(5) −0.021 82(8) 
100 −0.007 23(3) 0.001 5(1) 2.407(4) −0.007 23(3) 
200 −0.002 418(8) 0.000 03(5) 2.353(4) −0.002 418(8) 
400 −0.000 567(3) −0.000 007(9) 2.343(4) −0.000 567(3) 

The Cr2 dimer is a more difficult problem than C2 and can be taken as a prototype of a challenging molecule in quantum chemistry. The corresponding results are shown in Table III and Fig. 1, with orbitals determined by a CASSCF calculation with CAS(12e,12o). In the stochastic calculations, the bond dimension for compression is M1 = 8000, and the number of samples is Ns = 28 000 except for one column where 10 times more samples were used to illustrate the convergence of the stochastic error. From Table III, we see that the usual variational DMRG converges very slowly, and even with M = 8000 the error is about 1mEh compared with the extrapolated result. By contrast, the convergence to chemical accuracy using p-DMRG and its stochastic variant is extremely fast using λ=12. As shown in Table III, for p-DMRG, chemical accuracy is achieved at M0 = 300, while for stochastic p-DMRG it is achieved at M0 = 200 (as its energy is usually slightly lower due to the same reasons as discussed for C2). The convergence with λ = 0 is considerably slower compared with the λ=12 case, which we found in our previous study also, due to the unbalanced treatment of the zeroth order state and determinant-like perturbers.26 However, it is still significantly cheaper than the variational DMRG for the same level of accuracy due to the smaller bond dimension M0. All stochastic evaluations took less than 15 min (wall time on a 28-core node) for a statistical error of less than 1mEh, while the variational compression for QV^|Ψ(0) took about 30 min with M1 = 8000. Thus, in total, the stochastic algorithm is about 5 times faster than the deterministic algorithm and further parallelizes much better than the deterministic one and uses much less memory and disk.

TABLE III.

Total energy (E + 2086 in Eh) of Cr2 in the Ahlrichs’ SV basis (48 electrons in 42 orbitals) obtained by DMRG and stochastic p-DMRG.

DMRG energy (E in Eh) and discarded weights (w) with a reversed sweep schedule.
M5001000200040008000
w 5.4 × 10−4 2.1 × 10−4 1.1 × 10−4 5.8 × 10−5 2.8 × 10−5  
E −0.4955 −0.5052 −0.5108 −0.5140 −0.5158 −0.5174(3)a 
DMRG energy (E in Eh) and discarded weights (w) with a reversed sweep schedule.
M5001000200040008000
w 5.4 × 10−4 2.1 × 10−4 1.1 × 10−4 5.8 × 10−5 2.8 × 10−5  
E −0.4955 −0.5052 −0.5108 −0.5140 −0.5158 −0.5174(3)a 
Stochastic p-DMRG energy (Ns = 28 000) obtained with 28 000 samples.b
p-DMRGStochastic p-DMRGc
M0DMRGλ = 0λ=12λ = 0λ=12λ=12d
100 −0.4236 −0.5000 −0.5147 −0.5025(3) −0.5177(4) −0.5182(1) 
200 −0.4568 −0.5068 −0.5154 −0.5072(2) −0.5167(3) −0.5167(1) 
300 −0.4785 −0.5107 −0.5161 −0.5118(2) −0.5170(2) −0.51726(8) 
400 −0.4861 −0.5123 −0.5165 −0.5130(2) −0.5171(2) −0.51728(7) 
500 −0.4921 −0.5136 −0.5169 −0.5143(2) −0.5174(1) −0.51751(6) 
Stochastic p-DMRG energy (Ns = 28 000) obtained with 28 000 samples.b
p-DMRGStochastic p-DMRGc
M0DMRGλ = 0λ=12λ = 0λ=12λ=12d
100 −0.4236 −0.5000 −0.5147 −0.5025(3) −0.5177(4) −0.5182(1) 
200 −0.4568 −0.5068 −0.5154 −0.5072(2) −0.5167(3) −0.5167(1) 
300 −0.4785 −0.5107 −0.5161 −0.5118(2) −0.5170(2) −0.51726(8) 
400 −0.4861 −0.5123 −0.5165 −0.5130(2) −0.5171(2) −0.51728(7) 
500 −0.4921 −0.5136 −0.5169 −0.5143(2) −0.5174(1) −0.51751(6) 
a

The error bar for the extrapolated energy () is estimated16 as 1/5 of the difference between the extrapolation energy and M = 8000 energy.

b

All calculations were performed on a 28-core node. Values in parentheses are the statistical uncertainties.

c

Stochastic p-DMRG does not converge to p-DMRG with increasing sample sizes due to the use of different Ĥ0 (see text).

d

Results obtained with 10 times more samples.

FIG. 1.

Convergence of variational DMRG, p-DMRG, and stochastic p-DMRG energies as a function of M0 for Cr2 in Ahlrichs’ SV basis. Both p-DMRG and stochastic p-DMRG improve the convergence significantly, and the choice λ=12 gives better results than λ = 0.

FIG. 1.

Convergence of variational DMRG, p-DMRG, and stochastic p-DMRG energies as a function of M0 for Cr2 in Ahlrichs’ SV basis. Both p-DMRG and stochastic p-DMRG improve the convergence significantly, and the choice λ=12 gives better results than λ = 0.

Close modal

We found, in practice, that the bond dimension M1 for compression does not need to be too large to achieve chemical accuracy. In Table IV, the dependence of the stochastic p-DMRG results on M1 for |Φ⟩ are shown for M0 = 500, λ = 0, and Ns = 28 000. Using Eq. (12), M1 = 4000 is sufficient to converge to 1mEh accuracy. For too small M1, due to the problem of a poor importance sampling function, which limits the set of the interacting determinants to be explored, the energy is slightly higher than the converged ones with larger M1. It is interesting to see that the approximate estimate (13) requires much larger bond dimension, approximately by a factor of 4, to achieve the same level of accuracy as Eq. (12), although its statistical error is smaller with the same Ns. Therefore, in practice, this approximation is not very useful.

TABLE IV.

Dependence of the total energy (E + 2086 in Eh) obtained by stochastic p-DMRG on the bond dimension M1 used to compress QV^|Ψ(0) as an MPS |Φ⟩. The results are illustrated with M0 = 500 and λ = 0.

M1E2 from Eq. (12)E2 from Eq. (13)
1000 −0.5100(15) −0.501 44(3) 
2000 −0.5122(14) −0.507 35(6) 
4000 −0.5143(7) −0.511 51(7) 
6000 −0.5148(7) −0.512 72(8) 
8000 −0.5143(2) −0.513 30(8) 
M1E2 from Eq. (12)E2 from Eq. (13)
1000 −0.5100(15) −0.501 44(3) 
2000 −0.5122(14) −0.507 35(6) 
4000 −0.5143(7) −0.511 51(7) 
6000 −0.5148(7) −0.512 72(8) 
8000 −0.5143(2) −0.513 30(8) 

Finally, we compare the performance of the stochastic p-DMRG against the deterministic p-DMRG for a larger problem, Cr2 in an orbital space of (28e, 76o), in Table V. For this system, the convergence of the previous deterministic p-DMRG is very slow. Even with a sum of five MPS with bond dimension M1 = 7500 for the first-order wavefunction, |Ψ(1)=i=15|Ψi(1), the deterministic p-DMRG energies differ from the extrapolated energy (M1 = ) by 5–10mEh depending on M0. However, the energies provided by the stochastic p-DMRG with a given M1 = 10 000 and Ns = 80 000 are already very close to the extrapolated M1 = p-DMRG energies, demonstrating the efficiency of the stochastic p-DMRG. As we showed in our deterministic p-DMRG calculations,26 one can further extrapolate in M0 and M1 to obtain an estimate of the exact energy without any truncation error. In future, we will explore the possibility of such extrapolations with stochastic p-DMRG also.

TABLE V.

Total energies (E + 2099 in Eh) of Cr2 in the cc-pVDZ-DK basis set (28 electrons in 76 orbitals) computed by deterministic and stochastic p-DMRG with λ=12.

M01000200030004000
DMRG −0.8346 −0.8617 −0.8743 −0.8818 
p-DMRGa (M1 = 5 × 7500) −0.9036 −0.9035 −0.9035 −0.9037 
p-DMRGb (M1 = −0.9080 −0.9109 −0.9129 −0.9141 
stochastic p-DMRG −0.905(2) −0.909(2) −0.909(1) −0.911(1) 
M01000200030004000
DMRG −0.8346 −0.8617 −0.8743 −0.8818 
p-DMRGa (M1 = 5 × 7500) −0.9036 −0.9035 −0.9035 −0.9037 
p-DMRGb (M1 = −0.9080 −0.9109 −0.9129 −0.9141 
stochastic p-DMRG −0.905(2) −0.909(2) −0.909(1) −0.911(1) 
a

Reference 26: energy for |Ψ(1)=i=15|Ψi(1) with M1 = 7500.

b

Reference 26: extrapolated energy for M1 = .

In summary, we have presented an efficient stochastic algorithm to overcome the previous bottleneck in p-DMRG,26 developed recently for problems with large active spaces. We demonstrated that, in combination with a good choice of zeroth order Hamiltonian, the stochastic p-DMRG algorithm can provide highly accurate total energies for challenging systems with significantly reduced computational resources as compared with the deterministic p-DMRG, and both of these are much cheaper than the original variational DMRG, in large orbital spaces with a mix of correlation strengths.

Note:During the finalization of this work, we became aware of a related study in Ref. 43 . This work assumed that E2can be approximated by the term A in Eq.(8)and used a different way to compute E2where(0)was represented stochastically along the lines of the original stochastic HCI+PT.34The systematic accuracy of the final energies from that algorithm should be comparable to the p-DMRG energies with the choice λ = 0, although the statistical properties will not be the same.

This work was supported by the US National Science Foundation through No. CHE 1665333. Additional support was provided by No. OAC 1657286. Z.L. is supported by the Simons Collaboration on the Many-Electron Problem. G.K.-L.C. is a Simons Investigator in Physics.

1.
K.
Andersson
,
P. A.
Malmqvist
,
B. O.
Roos
,
A. J.
Sadlej
, and
K.
Wolinski
,
J. Phys. Chem.
94
,
5483
(
1990
).
2.
B. O.
Roos
,
K.
Andersson
,
M. P.
Flscher
,
P.-Å.
Malmqvist
,
L.
Serrano-Andrs
,
K.
Pierloot
, and
M.
Merchn
, in
Advances in Chemical Physics
, edited by
I.
Prigogine
and
S. A.
Rice
(
John Wiley & Sons, Inc.
,
1996
), pp.
219
331
.
3.
C.
Angeli
,
R.
Cimiraglia
,
S.
Evangelisti
,
T.
Leininger
, and
J.-P.
Malrieu
,
J. Chem. Phys.
114
,
10252
(
2001
).
4.
C.
Angeli
,
R.
Cimiraglia
, and
J.-P.
Malrieu
,
Chem. Phys. Lett.
350
,
297
(
2001
).
5.
C.
Angeli
,
R.
Cimiraglia
, and
J.-P.
Malrieu
,
J. Chem. Phys.
117
,
9138
(
2002
).
6.
K.
Andersson
and
B. O.
Roos
,
Chem. Phys. Lett.
191
,
507
(
1992
).
7.
Y.
Kurashige
and
T.
Yanai
,
J. Chem. Phys.
135
,
094104
(
2011
).
8.
S.
Guo
,
M. A.
Watson
,
W.
Hu
,
Q.
Sun
, and
G. K.-L.
Chan
,
J. Chem. Theory Comput.
12
,
1583
(
2016
).
9.
10.
S. R.
White
,
Phys. Rev. B
48
,
10345
(
1993
).
11.
S. R.
White
and
R. L.
Martin
,
J. Chem. Phys.
110
,
4127
(
1999
).
12.
A. O.
Mitrushenkov
,
G.
Fano
,
F.
Ortolani
,
R.
Linguerri
, and
P.
Palmieri
,
J. Chem. Phys.
115
,
6815
(
2001
).
13.
G. K.-L.
Chan
and
M.
Head-Gordon
,
J. Chem. Phys.
116
,
4462
(
2002
).
14.
Ö.
Legeza
,
J.
Röder
, and
B. A.
Hess
,
Phys. Rev. B
67
,
125114
(
2003
).
15.
S.
Sharma
and
G. K.-L.
Chan
,
J. Chem. Phys.
136
,
124121
(
2012
).
16.
R.
Olivares-Amaya
,
W.
Hu
,
N.
Nakatani
,
S.
Sharma
,
J.
Yang
, and
G. K.-L.
Chan
,
J. Chem. Phys.
142
,
034102
(
2015
).
17.
S.
Keller
,
M.
Dolfi
,
M.
Troyer
, and
M.
Reiher
,
J. Chem. Phys.
143
,
244118
(
2015
).
18.
T.
Yanai
,
Y.
Kurashige
,
W.
Mizukami
,
J.
Chalupský
,
T. N.
Lan
, and
M.
Saitow
,
Int. J. Quantum Chem.
115
,
283
(
2015
).
19.
G. K.-L.
Chan
,
A.
Keselman
,
N.
Nakatani
,
Z.
Li
, and
S. R.
White
,
J. Chem. Phys.
145
,
014102
(
2016
).
20.
S.
Sharma
and
G. K.-L.
Chan
,
J. Chem. Phys.
141
,
111101
(
2014
).
21.
N.
Nakatani
and
S.
Guo
,
J. Chem. Phys.
146
,
094102
(
2017
).
22.
S.
Sharma
,
G.
Knizia
,
S.
Guo
, and
A.
Alavi
,
J. Chem. Theory Comput.
13
,
488
(
2017
).
23.
A. Y.
Sokolov
and
G. K.-L.
Chan
,
J. Chem. Phys.
144
,
064102
(
2016
).
24.
L.
Freitag
,
S.
Knecht
,
C.
Angeli
, and
M.
Reiher
,
J. Chem. Theory Comput.
13
,
451
(
2017
).
25.
A. Y.
Sokolov
,
S.
Guo
,
E.
Ronca
, and
G. K.-L.
Chan
,
J. Chem. Phys.
146
,
244102
(
2017
).
26.
S.
Guo
,
Z.
Li
, and
G. K.-L.
Chan
, e-print arXiv:1803.07150.
27.
B.
Huron
,
J.
Malrieu
, and
P.
Rancurel
,
J. Chem. Phys.
58
,
5745
(
1973
).
28.
R. J.
Buenker
and
S. D.
Peyerimhoff
,
Theor. Chem. Acc.
35
,
33
(
1974
).
29.
R. J.
Harrison
,
J. Chem. Phys.
94
,
5021
(
1991
).
30.
J. B.
Schriber
and
F. A.
Evangelista
,
J. Chem. Phys.
144
,
161106
(
2016
).
31.
N. M.
Tubman
,
J.
Lee
,
T. Y.
Takeshita
,
M.
Head-Gordon
, and
K. B.
Whaley
,
J. Chem. Phys.
145
,
044112
(
2016
).
32.
W.
Liu
and
M. R.
Hoffmann
,
J. Chem. Theory Comput.
12
,
1169
(
2016
).
33.
A. A.
Holmes
,
N. M.
Tubman
, and
C. J.
Umrigar
,
J. Chem. Theory Comput.
12
,
3674
(
2016
).
34.
S.
Sharma
,
A. A.
Holmes
,
G.
Jeanmairet
,
A.
Alavi
, and
C. J.
Umrigar
,
J. Chem. Theory Comput.
13
,
1595
(
2017
).
35.
Y.
Garniron
,
A.
Scemama
,
P.-F.
Loos
, and
M.
Caffarel
,
J. Chem. Phys.
147
,
034101
(
2017
).
36.
37.
E. M.
Stoudenmire
and
S. R.
White
,
New J. Phys.
12
,
055026
(
2010
).
38.
M.
Douay
,
R.
Nietmann
, and
P.
Bernath
,
J. Mol. Spectrosc.
131
,
250
(
1988
).
39.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
40.
V. E.
Bondybey
and
J. H.
English
,
Chem. Phys. Lett.
94
,
443
(
1983
).
41.
A.
Schäfer
,
H.
Horn
, and
R.
Ahlrichs
,
J. Chem. Phys.
97
,
2571
(
1992
).
42.
N. B.
Balabanov
and
K. A.
Peterson
,
J. Chem. Phys.
123
,
064107
(
2005
).
43.
S.
Sharma
, preprint arXiv:1803.04341 (
2018
).