bug report of the multi-quantum number code

[jiangtong1000]

146 added the multi-quantum number support.

Recently I am using that part and found something weird, I will call the package starting from that commit as new and the version before as old. And here is what I found for electronic structure calculations. Apparently there is a bug.

H10 (3 Bohr) old: -4.974258044176154 (M=30) new: -4.967491721857599 (M=30), -4.967422012064133 (M=50) FCI: -4.974243429294

H20 (3 Bohr) old: -9.952918380537831 (M=30) new: -9.945641599914602 (M=30)

And If I use M=50 for H20, there will be quantum number error reported when doing svd_qn.

[jiangtong1000]

from renormalizer.model import h_qc, Model
from renormalizer.mps import Mpo, Mps, optimize_mps
from pyscf import lo, tools, gto
import tempfile
import numpy as np

nocca = 5
noccb = 5
nelec = nocca + noccb
r0 = 3.0
mol = gto.M(
    atom=[("H", i * r0, 0, 0) for i in range(nelec)],
    basis='sto-6g',
    unit='Bohr',
    verbose=5
)
mol_nelec = [nocca, noccb]
screen_tol = 1.e-5

bonddim = 30

s1e = mol.intor("int1e_ovlp_sph")
ao_coeff = lo.orth.lowdin(s1e)
spatial_norbs = ao_coeff.shape[0]
ftmp = tempfile.NamedTemporaryFile()
tools.fcidump.from_mo(mol, ftmp.name, ao_coeff)
h1e, h2e, nuc = h_qc.read_fcidump(ftmp.name, spatial_norbs)
h1e[np.abs(h1e) < screen_tol] = 0
h2e[np.abs(h2e) < screen_tol] = 0
basis1, ham_terms1 = h_qc.qc_model(h1e, h2e)
model = Model(basis1, ham_terms1)
mpo = Mpo(model)
nelec = sum(mol.nelec)
np.random.seed(7)
energy_list = []

mps = Mps.random(model, [nelec//2, nelec//2], bonddim, percent=1.0)
mps.optimize_config.method = '2site'
mps.optimize_config.procedure = [[bonddim, 0.4], [bonddim, 0.2], [bonddim, 0.1], [bonddim, 0],
                                 [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0],
                                 [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0],
                                 [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0], [bonddim, 0]]
energies, mps = optimize_mps(mps.copy(), mpo)
print(mps.expectation(mpo) + nuc)

[jiangtong1000]

If you like to try, the above is the code. Please check between the commit bac8bed3d1567 (the older one) and commit ed424bd5854d32e2 (the newer one). ps: use the trick #158 to speed up the mpo construction, especially for H20.

[jiangtong1000]

Also I remembered there is initial guess sensitivity worries mentioned at https://github.com/shuaigroup/Renormalizer/pull/146#issuecomment-1369393017 Maybe it is because the bug too

[jiangtong1000]

At least two bugs should be fixed, otherwise there will be qn_alpha exceeds qntot_alpha and negative qn_beta number

1. https://github.com/shuaigroup/Renormalizer/blob/1576c3c2bc16ffb635d751fd4edcce5d5752343b/renormalizer/mps/mps.py#L140-L142 should be

   if qntot[0] < iblock[0] or qntot[1] < iblock[1]:
       continue

2. https://github.com/shuaigroup/Renormalizer/blob/1576c3c2bc16ffb635d751fd4edcce5d5752343b/renormalizer/mps/svd_qn.py#L176-L178 should be

     if qntot[0] <nl[0] or qntot[1] < nl[1]:
              continue

[jiangtong1000]

The results are getting better, but still not totally fixed.

H-10 (r=3 Bohr) bugfix: 10: -4.960871569710858 30: -4.974115514369827 50: -4.974278368650478 70: -4.97427932154667 90: -4.974279544655352

old: 10: -4.9718418455855025 30: -4.974258806240129 50: -4.97427852927401 70: -4.974279439932787 90: -4.974279543996439

H20 (r=3 Bohr) bugfix: 10: -9.599322484758492 30: -9.358948827767495 50: -9.95273240869022

old: 10: -9.94530042750198 30: -9.952918182626156 50: -9.953052459821894

[jiangtong1000]

I went through the code, and wasn't able to find another bug again. Here are some numerical tests. With more careful design of sweep procedure (by using larger pencent value for more initial sweeps, H20's energy for M=50 can agree with the old code up to 1e-6. While for M=30, even with more initial sweeps with larger percent value, it can't agree to the old one (0.5 Hartree away), unless using the M=50's converged MPS as the initial guess for it, then it can agree with the old code's result up to 1e-6. I am still concerning about the possible bugs, but the numerical results make sense to me, since restricting the alpha and beta spins actually reduced the variational space, as compared to only restricting the total spins, then it may become easier to be trapped in local minimum. If this is true, I think retaining the old feature becomes more important, + https://github.com/shuaigroup/Renormalizer/discussions/160#discussion-5525300.

[liwt31]

Thanks for the report and the investigation.

And If I use M=50 for H20, there will be quantum number error reported when doing svd_qn.

Can you paste the whole error message (if still accessible)?

At least two bugs should be fixed, otherwise there will be qn_alpha exceeds qntot_alpha and negative qn_beta number

This is actually a known issue that does not have a good fix. In the context of electronic structure, negative quantum numbers do not make any sense. But in a more general context, it is possible to have negative quantum numbers. Consider, for example, a system that has conserved $S_z=0$.

So I guess maybe it's a better solution to provide another specific initialization function for electronic structure problems. Another issue is the basis set used for the calculation. In our code we use Spin-1/2 directly. However the more common practice is to use a site with {0, 1-alpha, 1-beta, 2} states. Thus our variational space is smaller than the common practice, which makes the algorithm less robust.

[jiangtong1000]

Can you paste the whole error message (if still accessible)?

That error has been removed if fixing the random function as mentioned.

This is actually a known issue that does not have a good fix. In the context of electronic structure, negative quantum numbers do not make any sense. But in a more general context, it is possible to have negative quantum numbers. Consider, for example, a system that has conserved Sz=0.

Thanks for letting me know.

However the more common practice is to use a site with {0, 1-alpha, 1-beta, 2} states. Thus our variational space is smaller than the common practice, which makes the algorithm less robust.

I thought fixing qntot has larger variational space then both the DoF=4/2 choices? And the extra components will not contribute to the final ground state but will actually act as noises, to some extent preventing from the local minimum.

So I guess maybe it's a better solution to provide another specific initialization function for electronic structure problems.

Sounds like a good compromise for now.

[jjren]

1. The random initialization is rather crude. If you know the approximate initial state, e.g. the HF state or antiferromagnetic state, you can use it and expand the bond dimension with expand_bond_dimension.
2. Of course, bad initialization will bring difficulty in optimization. But, In my opinion, the bugs Tong reported are not bugs, they will not lead to wrong results because I guess the wrong qn sectors will be removed automatically in the DMRG sweeping.
3. There are more robust algorithms than the one we use to select basis (with argument percent), for example, Phys. Rev. B 72, 180403. If getting stuck in a local minimum really matters to you, I think you can try to implement this algorithm.

[liwt31]

Maybe a first step towards solving the issue is to try different initialization methods - If by using a proper initialization, such as the method suggested by Jiajun, multi-quantum number implementation is able to produce similar results to single-quantum number implementation, then it is likely not a bug but an initialization issue. If with a decent initial guess the algorithm still does not work, then we can dig deeper and try to discover the bug.

[jiangtong1000]

I see. There is an example which used hf initial guess, but I think the mps orbitals are MO, right? I think I can introduce orho ao way at sometime later.

[jjren]

With Lowdin orth ao, I guess the initial state can be an antiferromagnetic state |up down up down> + |down up down up>.

[jiangtong1000]

I am closing this issue since all the questions have been resolved. Indeed, what I reported is not a bug. Instead, if we really restrict the quantum number sector as I suggested above. The result will be wrong for some other cases (ground state calls are still good) since we also use the svd_qn for canonicalization, therefore, the direct canonicalise/compression of a general MPO or a MPS that is not ground state might have errors.