How to deal with block getting trapped in too high states

[ppinski-hqs]

When performing calculations with Block (as with other FCI solvers), there is some risk of getting trapped in too high electronic states. One rather notorious case I encountered was for the toy example of a carbon atom in STO-3G basis.

My questions upfront:

• Is the behaviour demonstrated below "normal", or is it likely there are some problems with Block itself or with its PySCF interface?
• If this is normal behaviour, do you have suggestions on how to deal with the problem more systematically than by cranking up the number of roots as I did below? (Raising the question at what point to stop with more complicated problems.)

I'm looking forward to to any helpful advice.

---

The minimal example code accessing Block via its PySCF interface below demonstrates the problem:

from pyscf.gto import Mole
from pyscf.scf import ROHF
from pyscf.mcscf import CASCI
from pyscf.dmrgscf import DMRGCI

# Carbon atom with STO-3G basis and triplet spin.
mol = Mole(atom='C 0 0 0', basis='STO-3G', spin=2, verbose=4).build()
mf = ROHF(mol).run()

# (1) Run FCI via CASCI: triply degenerate ground state, as expected.
mc = CASCI(mf, 5, 6)
mc.fcisolver.nroots = 3
mc.kernel()

# (2) Run DMRG with three roots: degeneracy observed with normal CI solver not reproduced.
mc = CASCI(mf, 5, 6)
mc.fcisolver = DMRGCI(mol)
mc.fcisolver.nroots = 3
mc.kernel()

# (3) Increasing the number of roots all the way to 10: now I get a triply degenerate ground state.
mc = CASCI(mf, 5, 6)
mc.fcisolver = DMRGCI(mol)
mc.fcisolver.nroots = 10
mc.kernel()

I ran the calculations using PySCF 2.4.0 and block2 0.5.2 from pypi.org for python 3.11. dmrgscf was in its latest version (4ff57bf). The results:

(1) CASCI: triply degenerate energy, as expected.

CASCI state   0  E = -37.2187335506365  E(CI) = -37.2187335506365  S^2 = 2.0000000
CASCI state   1  E = -37.2187335506365  E(CI) = -37.2187335506365  S^2 = 2.0000000
CASCI state   2  E = -37.2187335506365  E(CI) = -37.2187335506365  S^2 = 2.0000000

(2) DMRG-CASCI with three roots: one of the energies corresponds to the excited state.

CASCI state   0  E = -37.2187335505787  E(CI) = -37.2187335505787  S^2 = 2.0000000
CASCI state   1  E = -37.2187335505787  E(CI) = -37.2187335505787  S^2 = 2.0000000
CASCI state   2  E = -36.8693631248328  E(CI) = -36.8693631248328  S^2 = 2.0000000

(3) DMRG-CASCI with ten roots (four to nine did not solve the problem): first three states are degenerate.

CASCI state   0  E = -37.2187335506362  E(CI) = -37.2187335506362  S^2 = 2.0000000
CASCI state   1  E = -37.2187335504622  E(CI) = -37.2187335504622  S^2 = 2.0000000
CASCI state   2  E = -37.2187335497707  E(CI) = -37.2187335497707  S^2 = 2.0000000
CASCI state   3  E = -36.8693631250086  E(CI) = -36.8693631250086  S^2 = 2.0000000
CASCI state   4  E = -36.8693631250084  E(CI) = -36.8693631250084  S^2 = 2.0000000
CASCI state   5  E = -36.8693631250034  E(CI) = -36.8693631250034  S^2 = 2.0000000
CASCI state   6  E = -36.8693631249764  E(CI) = -36.8693631249764  S^2 = 2.0000000
CASCI state   7  E = -36.8693631248167  E(CI) = -36.8693631248167  S^2 = 2.0000000
CASCI state   8  E = -36.7970258229902  E(CI) = -36.7970258229902  S^2 = 2.0000000
CASCI state   9  E = -36.7970258229130  E(CI) = -36.7970258229130  S^2 = 2.0000000

[hczhai]

Hi,

Thanks for reporting the problem and providing the script. I can reproduce the result you observed. I think the behavior is normal considering the following facts:

1. DMRG is not a black-box method, and it has many adjustable parameters.
2. The "best" parameters for DMRG should be system dependent. Namely, they will not be the same for large and small systems.
3. In the example script you provide, you implicitly used the default parameters provided by dmrgscf. And it is reasonable that the default parameters is only good for some medium-sized systems, but not for tiny systems.

In the block2 documentation, you can find a concrete set of recommended parameters for small systems:

# set very tight thresholds for small system
mc.fcisolver.scheduleSweeps = [0, 4, 8, 12, 16]
mc.fcisolver.scheduleMaxMs = [250, 500, 500, 500, 500]
mc.fcisolver.scheduleTols = [1e-08, 1e-10, 1e-12, 1e-12, 1e-12]
mc.fcisolver.scheduleNoises = [0.0001, 0.0001, 5e-05, 5e-05, 0.0]
mc.fcisolver.maxIter = 30
mc.fcisolver.twodot_to_onedot = 20
mc.fcisolver.block_extra_keyword = ['singlet_embedding', 'full_fci_space', 'fp_cps_cutoff 0', 'cutoff 0']
mc.fcisolver.conv_tol = 1e-14
mc.conv_tol = 1e-11

With this set of parameters, the discrepancy between FCI and DMRG is below 1E-14 for all three states. For larger system I do not know the best parameters and you have to do the test for your own case. But for medium-to-large systems it is less likely to get trapped in high energy local minima. You can also read the following paper to get a better understanding of the DMRG parameters:

Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. The ab-initio density matrix renormalization group in practice. The Journal of Chemical Physics 2015, 142, 034102. doi: 10.1063/1.4905329

[ppinski-hqs]

Thank you for the quick reply and for the advice.

Indeed, we typically use the default parameters, with the exception of very few settings such as the bond dimension. The goal is to calculate entropies to help with active space selection, so the results only need to be "good enough" without requiring quantitative accuracy. Most of the time it works fine, except for those cases when the calculation gets trapped in solutions as described above. Possibly this thinking was too simplistic, and we may need to invest a bit more effort into suitable DMRG settings, after all.

Thank you also for providing block2 as a pip package - that has made installing and using it quite convenient.