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 C_{2} and Cr_{2} 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 CASPT2^{1,2} (complete-active-space second-order perturbation theory) and NEVPT2^{3–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 3*d* transition metal complexes, the virtual 4*d*, semi-core (3*s* and 3*p*), 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 *M*_{0} (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 *M*_{0}, 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 *M*_{0} 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 *M*_{0}, has been optimized by the standard variational DMRG algorithm. Given a partitioning of the Hamiltonian, $\u0124=\u01240+V^$ with *Ĥ*_{0}|Ψ^{(0)}⟩ = *E*_{0}|Ψ^{(0)}⟩, the first order wavefunction |Ψ^{(1)}⟩ is defined by the first order equation,

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=\u27e8\Psi (1)|QV^|\Psi (0)\u27e9$ can be computed. For very large active spaces with 50–100 orbitals, the required bond dimension *M*_{1} for expanding |Ψ^{(1)}⟩, which typically scales like *O*(*K*^{2}*M*_{0}) (*K* is the number of spatial orbitals), can be quite large (≫10^{4}). 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 $|\Psi (1)=\u2211i|\Psi 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 *E*_{2} as

and then aim to find an explicit expression for *E*_{2} as a sum over terms, ideally of a similar form to $\u2211i\u2212|\u27e8Di|V^|\Psi (0)\u27e9|2Ei\u2212E0$ as appears in HCI+PT with |*D*_{i}⟩ being a determinant, which can then be sampled stochastically. In our previous paper,^{26} we chose *Ĥ*_{0} as

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

with $\xcapq=\u2211\sigma ap\sigma \u2020aq\sigma $ and $\xeapqrs=\u2211\sigma ,\tau ap\sigma \u2020aq\tau \u2020ar\tau as\sigma =EpsEqr\u2212\delta 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 ⟨*D*_{i}|*Ĥ*_{d}|*D*_{j}⟩ = *δ*_{ij}⟨*D*_{i}|*Ĥ*|*D*_{j}⟩ ≜ *δ*_{ij}*E*_{i}. This choice will result in a slight difference in the computed *E*_{2} compared with our previous p-DMRG, usually found to be less than 1m*E*_{h} (vide post), and the difference will gradually decrease as *M*_{0} increases.

The EN partition is, of course, commonly used in SCI+PT.^{27–35} In this case, $Q[Q(\u01240\u2212E0)Q]\u22121Q$ in (2) can be simplified to $Q(\u0124d\u2212E0)\u22121Q$ for *Q* = *∑*_{i}|*D*_{i}⟩⟨*D*_{i}| and [*Q*, *Ĥ*_{d}] = 0, giving $E2=\u2211i\u2212|\u27e8Di|V^|\Psi (0)\u27e9|2Ei\u2212E0$ as a sum over determinants in the *Q*-space, as mentioned above. However, in our case, *Q*(*Ĥ*_{d} − *E*_{0})*Q* and hence $Q[Q(\u01240\u2212E0)Q]\u22121Q$ 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 *E*_{2}, we first rewrite

In general, the operator *X* = *Ĥ*_{d} − *E*_{0} can be assumed to be invertible by properly choosing *E*_{0}, and this important issue will be discussed later. Then, using the relation $Q(X\u2212Y)\u22121Q\u2009=\u2009QX\u22121(1\u2212YX\u22121)\u22121Q\u2009=\u2009QX\u22121Q\u2009+\u2009QX\u22121YX\u22121Q+\u2009QX\u22121YX\u22121YX\u22121Q\u2009+\u2009\cdots \u2009$ and directly evaluating each term, we can find that the series can be greatly simplified as *P* is one-dimensional such that *Q*(*X* − *Y*)^{−1}*Q* becomes

which can be simply verified by multiplying the right-hand side with *Q*(*Ĥ*_{0} − *E*_{0})*Q* to get the identity operator *Q* in the *Q*-space. The last term is the correction for the fact that $Q(\u0124d\u2212E0)\u22121P$ 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* = (*Ĥ*_{d} − *E*_{0}) 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 *E*_{2} and invoking the resolution of identity *∑*_{i}|*D*_{i}⟩⟨*D*_{i}| = 1, each term in *E*_{2} can finally be formulated as a sum over determinants, viz.,

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

In the practical evaluation of *A*, *B*, and *C*, instead of sampling |*D*_{i}⟩ 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=\u27e81Ei\u2212E0\u27e9Pi=|\u27e8\Psi (0)|Di\u27e9|2$, where the subscript indicates that the average of $1Ei\u2212E0$ is taken with respect to the population of the determinants generated with the probability *P*_{i} = |⟨Ψ^{(0)}|*D*_{i}⟩|^{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,

where $C\alpha 1n1[1]$ and $R\alpha K\u22121nK[K]$ are matrices and $R\alpha k\u22121\alpha knk[k]$ are rank-3 tensors satisfying the right canonical condition $\u2211nk\alpha kR\alpha k\u22121\u2032\alpha knk[k]R\alpha k\u22121\alpha knk[k]=\delta \alpha k\u22121\u2032\alpha k\u22121$, then a determinant |*D*_{i}⟩ = |*n*_{1}*n*_{2}…*n*_{K}⟩ can be sampled according to *P*_{i} by a single sweep from left to right as follows: The first occupation number *n*_{1} ∈{|$v$*ac*⟩, |1_{β}⟩, |1_{α}⟩, |1_{α}1_{β}⟩} is generated according to $p1(n1)\u225c\u2211\alpha 1|C\alpha 1n1|2$, which satisfies $\u2211n1p1=1$ as |Ψ^{(0)}⟩ is normalized. Given *n*_{1}, at the second site, the generation probability for *n*_{2} is defined as $p2(n2|n1)\u225c\u2211\alpha 2|(Cn1Rn2)\alpha 2|2/N2$. Importantly, the normalization constant *N*_{2} can be found as $N2=\u2211n2p2(n2|n1)=p1(n1)$ due to the right canonical property. Repeating this procedure recursively, the generation probability for *n*_{k} given *n*_{1}⋯*n*_{k−1} can be defined as $pk(nk|n1\cdots nk\u22121)\u225c\u2211\alpha 1|(Cn1Rn2\cdots Rnk)\alpha k|2/Nk$ with *N*_{k} = *p*_{k−1}(*n*_{k−1}|*n*_{1}⋯*n*_{k−2}). Thus, after the occupation of the last site is chosen, the total generation probability is $P(n1\cdots nK)=p1(n1)p2(n2|n1)\cdots p(nk|n1\cdots nK\u22121)=|(Cn1Rn2\cdots RnK)|2=|\u27e8\Psi (0)|Di\u27e9|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 $CSk\u22121Mk\u22121skmkSkMk$ 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^|\Psi (0)$ can be expressed as an MPS faithfully. However, since converting $QV^|\Psi (0)$ into an MPS exactly would incur a bond dimension of *O*(*K*^{2}*M*_{0}), this becomes prohibitive for large active spaces. Thus, we first use a bond dimension *M*_{1} [that is small compared with *O*(*K*^{2}*M*_{0})] to compress $QV^|\Psi (0)$ variationally as an MPS, i.e., $|\Phi \u2248QV^|\Psi (0)$. This approximate state |Φ⟩ is then used to define a generation probability for |*D*_{i}⟩, *P*_{i} = |⟨Φ|*D*_{i}⟩|^{2}, for importance sampling, such that

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

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 $|\Phi =QV^|\Psi (0)$. It deserves to be noted that in Eq. (12), there is a subtlety. For the equality to hold, the set of determinants $S\u225c{|Di:\u27e8\Phi |Di\u27e9\u22600}$ interacting with |Ψ^{(0)}⟩ must include the set that interacts with $\u27e8\Psi (0)|V^Q|Di\u27e9\u22600$, otherwise some contributions to *E*_{2} will be missing. This is guaranteed for sufficiently large *M*_{1}, and fortunately, we found that to converge to chemical accuracy ca. 1m*E*_{h}, the required *M*_{1} can be much smaller than *O*(*K*^{2}*M*_{0}) (vide post). For the term *B*(9), we simply evaluate it as $B=\u27e8\u27e8\Psi (0)|V^Q|Di\u27e9\u27e8Di|\Psi (0)\u27e9Pi(Ei\u2212E0)\u27e9Pi=|\u27e8\Phi |Di\u27e9|2$ since it turns out to be less important in most cases due to its small size and the quadratic dependence of *E*_{2} on it.

The above algorithm constitutes a stochastic algorithm for p-DMRG. Its computational cost depends on two parts: the compression for $QV^|\Psi (0)$ on the right-hand side of the first order equation (1), which scales as $O(K3M12M0)$ assuming *M*_{1} ≫ *M*_{0}, and the cost for the stochastic evaluation of *E*_{2}. The latter is dominated by the number of samples *N*_{s} times the cost for evaluating the matrix elements $\u27e8\Psi (0)|V^Q|Di\u27e9$. We found the following scheme to be efficient at least for the systems investigated in this work. By using the identity $\u27e8\Psi (0)|V^Q|Di\u27e9=\u27e8\Psi (0)|\u0124Q|Di\u27e9=\u27e8\Psi (0)|(H\u2212EDMRG(0))|Di\u27e9=\u2211j\u27e8\Psi (0)|Dj\u27e9\u27e8Dj|(H\u2212EDMRG(0))|Di\u27e9$, the evaluation is converted into *O*(*K*^{2}*N*^{2}) evaluations of overlaps between determinants and the zeroth order wavefunction ⟨Ψ^{(0)}|*D*_{j}⟩, 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 $\u27e8\Psi (0)|(H\u2212EDMRG(0))|Di\u27e9$ is directly computed by considering |*D*_{i}⟩ as an MPS with bond dimension one. However, the efficiency of our choice in practice may stem from the fact that in computing ⟨Ψ^{(0)}|*D*_{j}⟩, the sparsity in the MPS tensors can be utilized such that the actual computational time is less than the direct evaluation of $\u27e8\Psi (0)|(H\u2212EDMRG(0))|Di\u27e9$. 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 *M*_{1} 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: C_{2} at its equilibrium bond length 1.242 53 Å^{38} in the cc-pVDZ basis set^{39} and Cr_{2} at its equilibrium bond length 1.68 Å^{40} in the Ahlrichs’ SV basis^{41} and the cc-pVDZ-DK basis.^{42} For C_{2} and Cr_{2} in the SV basis, all electrons were correlated, resulting in orbital spaces with 12 electrons in 28 orbitals, (12e, 28o) for C_{2} and (48e,42o) for Cr_{2}, respectively. For Cr_{2} in the cc-pVDZ-DK basis, the same orbital space as in our previous work^{26} was employed, viz., (28e, 76o) with the 1*s*, 2*s*, and 2*p* frozen. As in our previous work on deterministic p-DMRG,^{26} the zeroth order energy is defined by an interpolation of two limits,

For simplicity, in this work, we only explored two cases, *λ* = 0 and $\lambda =12$, where in the former case $E0(\lambda =0)=EDMRG(0)$, and in the latter case $E0(\lambda =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 C_{2} obtained with stochastic p-DMRG with *N*_{s} = 36 000 are listed in Table I together with the p-DMRG results. (Note that we report *N*_{s} 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 *M*_{0}, and converge toward the exact value computed by variational DMRG with *M* = 4000 very quickly. The $\lambda =12$ energies are better than those with *λ* = 0 due to the reasons discussed earlier. The bond dimension *M*_{1} 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 *QĤ*_{d}*Q* these couplings will make the singlet perturbers have higher energies in the original p-DMRG compared with *E*_{i} 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 *M*_{0} = 400 is about 1 min for stochastic p-DMRG, which is 4 times faster than the deterministic one.

. | p-DMRG . | Stochastic p-DMRG^{a}
. | |||
---|---|---|---|---|---|

M_{0}
. | DMRG . | λ = 0
. | $\lambda =12$ . | λ = 0
. | $\lambda =12$ . |

Total energies of C_{2} | |||||

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 C_{2} | |||||

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-DMRG . | Stochastic p-DMRG^{a}
. | |||
---|---|---|---|---|---|

M_{0}
. | DMRG . | λ = 0
. | $\lambda =12$ . | λ = 0
. | $\lambda =12$ . |

Total energies of C_{2} | |||||

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 C_{2} | |||||

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 *E*_{2}, 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 *E*_{2} 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.

M_{0}
. | A
. | B
. | C
. | E_{2}
. |
---|---|---|---|---|

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

M_{0}
. | A
. | B
. | C
. | E_{2}
. |
---|---|---|---|---|

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 Cr_{2} dimer is a more difficult problem than C_{2} 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 *M*_{1} = 8000, and the number of samples is *N*_{s} = 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 1m*E*_{h} compared with the extrapolated result. By contrast, the convergence to chemical accuracy using p-DMRG and its stochastic variant is extremely fast using $\lambda =12$. As shown in Table III, for p-DMRG, chemical accuracy is achieved at *M*_{0} = 300, while for stochastic p-DMRG it is achieved at *M*_{0} = 200 (as its energy is usually slightly lower due to the same reasons as discussed for C_{2}). The convergence with *λ* = 0 is considerably slower compared with the $\lambda =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 *M*_{0}. All stochastic evaluations took less than 15 min (wall time on a 28-core node) for a statistical error of less than 1m*E*_{h}, while the variational compression for $QV^|\Psi (0)$ took about 30 min with *M*_{1} = 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.

DMRG energy (E in E_{h}) and discarded weights ($w$) with a reversed sweep schedule.
. | ||||||
---|---|---|---|---|---|---|

M
. | 500 . | 1000 . | 2000 . | 4000 . | 8000 . | ∞
. |

$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 E_{h}) and discarded weights ($w$) with a reversed sweep schedule.
. | ||||||
---|---|---|---|---|---|---|

M
. | 500 . | 1000 . | 2000 . | 4000 . | 8000 . | ∞
. |

$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 (N_{s} = 28 000) obtained with 28 000 samples.^{b}
. | ||||||
---|---|---|---|---|---|---|

. | . | p-DMRG . | Stochastic p-DMRG^{c}
. | |||

M_{0}
. | DMRG . | λ = 0
. | $\lambda =12$ . | λ = 0
. | $\lambda =12$ . | $\lambda =12$^{d}
. |

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 (N_{s} = 28 000) obtained with 28 000 samples.^{b}
. | ||||||
---|---|---|---|---|---|---|

. | . | p-DMRG . | Stochastic p-DMRG^{c}
. | |||

M_{0}
. | DMRG . | λ = 0
. | $\lambda =12$ . | λ = 0
. | $\lambda =12$ . | $\lambda =12$^{d}
. |

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 estimated^{16} 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.

We found, in practice, that the bond dimension *M*_{1} 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 *M*_{1} for |Φ⟩ are shown for *M*_{0} = 500, *λ* = 0, and *N*_{s} = 28 000. Using Eq. (12), *M*_{1} = 4000 is sufficient to converge to 1m*E*_{h} accuracy. For too small *M*_{1}, 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 *M*_{1}. 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 *N*_{s}. Therefore, in practice, this approximation is not very useful.

M_{1}
. | E_{2} from Eq. (12)
. | E_{2} 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, Cr_{2} 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 *M*_{1} = 7500 for the first-order wavefunction, $|\Psi (1)=\u2211i=15|\Psi i(1)$, the deterministic p-DMRG energies differ from the extrapolated energy (*M*_{1} = *∞*) by 5–10m*E*_{h} depending on *M*_{0}. However, the energies provided by the stochastic p-DMRG with a given *M*_{1} = 10 000 and *N*_{s} = 80 000 are already very close to the extrapolated *M*_{1} = *∞* 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 *M*_{0} and *M*_{1} 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.

M_{0}
. | 1000 . | 2000 . | 3000 . | 4000 . |
---|---|---|---|---|

DMRG | −0.8346 | −0.8617 | −0.8743 | −0.8818 |

p-DMRGa (M_{1} = 5 × 7500) | −0.9036 | −0.9035 | −0.9035 | −0.9037 |

p-DMRGb (M_{1} = ∞) | −0.9080 | −0.9109 | −0.9129 | −0.9141 |

stochastic p-DMRG | −0.905(2) | −0.909(2) | −0.909(1) | −0.911(1) |

M_{0}
. | 1000 . | 2000 . | 3000 . | 4000 . |
---|---|---|---|---|

DMRG | −0.8346 | −0.8617 | −0.8743 | −0.8818 |

p-DMRGa (M_{1} = 5 × 7500) | −0.9036 | −0.9035 | −0.9035 | −0.9037 |

p-DMRGb (M_{1} = ∞) | −0.9080 | −0.9109 | −0.9129 | −0.9141 |

stochastic p-DMRG | −0.905(2) | −0.909(2) | −0.909(1) | −0.911(1) |

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 E*_{2} *can be approximated by the term A in Eq.* (8) *and used a different way to compute E*_{2} *where* |Ψ^{(0)}⟩ *was represented stochastically along the lines of the original stochastic HCI+PT*.^{34} *The 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.*

## ACKNOWLEDGMENTS

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.