Details of the spin adaptation in block2

[bogdanoff]

Dear Huanchen, thank you very much for building up such a nice code and for providing very detailed documentation.

I am trying to make use of the CSF coefficients that can be calculating from the spin-adapted MPS state. Could you explain, what kind of spin-adaptation scheme is used, and how to read the CSF labels from the printout.

As a test system I chose N2 at some distance in using localized orbitals, see the attached input.
At the end I get the following printout:

DMRG energy = -107.438020323126224
LOC, reorder: [0 1 2 3 4 5]
dtrie finished 0.008127588080242276
Number of CSF =          3 (cutoff =      0.05)
Sum of weights of included CSF =    0.999999999727636

CSF          0 -++---  =   -0.745359113074040
CSF          1 +-+---  =   -0.527041871076544
CSF          2 ++----  =   -0.408248280360435

and

DMRG energy = -107.438020323076884
LOC, reorder: [0 3 1 4 2 5]
dtrie finished 0.008290691999718547
Number of CSF =          9 (cutoff =      0.05)
Sum of weights of included CSF =    0.999999999421752

CSF          0 -++---  =    0.527049368666756
CSF          1 +-+---  =    0.430331437547624
CSF          2 ++----  =    0.351364451030205
CSF          3 -+--+-  =   -0.304289567457688
CSF          4 --++--  =   -0.304289527965494
CSF          5 +--+--  =   -0.272164475050635
CSF          6 +---+-  =   -0.248451336091424
CSF          7 -+-+--  =   -0.248451238372456
CSF          8 ---++-  =    0.175678670854218

1. I can see that the CSF expansion depends on site ordering, as it should.
2. It looks like the position of - and + symbols correspond to orbitals after reordering, not to the indices in FCIDUMP.
3. I assume that + and - roughly corresponding to spin-up and spin-down coupling, but overall the strings are not of the UGA (unitary group approach) type, as there the CSF string can only start with u coupling, and along the string cannot go below 0. Such that duu and udd are not a valid UGA CSF labels.

from pyblock2._pyscf.ao2mo import integrals as itg
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
from pyscf import gto, scf

import numpy
from pyscf import lo, mcscf
from pyscf.tools import molden, mo_mapping

my_2ms = 2

mol = gto.M(atom="N1 0 0 -3; N2 0 0 3", basis="sto3g", symmetry="c1", verbose=0, spin=6)
mf = scf.RHF(mol).run(conv_tol=1E-14)

pm = lo.PM(mol, mf.mo_coeff, mf)
pm.pop_method = 'iao'
loc_mo = pm.kernel()

#sort                                                                                                                                                                                                              
tol = 0.9
s_idx = mo_mapping.mo_comps(["N1 2p","N2 2p"], mol, loc_mo) < tol
s_orbs = loc_mo[:, s_idx]

ord_labels = ["N{} 2p{}".format(i,x) for i in [1,2] for x in 'xyz']
mo_list = [s_orbs]
for label in ord_labels:
    mo_idx = mo_mapping.mo_comps([label], mol, loc_mo) > tol
    mo_list.append(loc_mo[:, mo_idx])
mo = numpy.hstack(mo_list)
molden.from_mo(mol, 'mo_'+str(my_2ms)+'.molden', mo)

ncas, n_elec, spin = 6, 6, my_2ms

mc = mcscf.CASCI(mf, ncas, n_elec)
mc.mo_coeff = mo
h1e, ecore = mc.get_h1eff()
g2e = mc.get_h2eff()


bond_dims = [250] * 4 + [500] * 5
noises = [1e-4] * 4 + [1e-5] * 5 + [0]
thrds = [1e-10] * 9

for reord in [numpy.array(range(6)), numpy.array([0,3,1,4,2,5])]:
    driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
    driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin)

    mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1, reorder=reord)

    ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
    energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
                         thrds=thrds, iprint=1)
    print('DMRG energy = %20.15f' % energy)

    print("LOC, reorder:",  reord)
    csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)

[hczhai]

Hi Nikolay,

Thanks for the good question. This part is indeed not clear in the documentation.

1. Since the CSF notation depends on the orbital ordering and reordering of the orbitals will generate CG factors, it is recommended not to use orbital reordering when you need the CSF strings.
2. By default, in the python interface, singlet embedding is used. The CSF string starts with a fictitious spin which is equal to the target spin. So the string "-++---" should be interpreted as "((((((1 - 1/2) + 1/2) + 1/2) - 1/2) - 1/2) - 1/2) = 0" where the first "1" is your target spin, and the right hand side will always be 0 (singlet).
3. If you take the reverse string and exchange "+/-" you may get a normal non-singlet-embedding string (the "UGA" type). So "-++---" becomes "+++--+". But I have not carefully check whether this gives the correct fermion sign.

[bogdanoff]

Thanks for the quick reply! Your explanation is clear, I need to think a bit how to deal with singlet embedding. I forgot about it and also thought that those orbitals are added at the end.

I think as the order matters, we can't look backwards to the string.

Strangely, the very last CSF in my example: ---++- that goes to uu|ddduud (uu| represent the singlet embedding here) still gives a bad label, as after third real electron S13=(((1-1/2)-1/2)-1/2) would be -1/2.

[hczhai]

Strangely, the very last CSF in my example: ---++- that goes to uu|ddduud (uu| represent the singlet embedding here) still gives a bad label, as after third real electron S13=(((1-1/2)-1/2)-1/2) would be -1/2.

The CSF output from block2 will not make sense if you use reorder=np.array([0,3,1,4,2,5]) in get_qc_mpo. This orbital reordering feature is not compatible with the CSF because as you know adjusting the CSF orbital order is a non-trivial thing. If you really need this orbital reordering, you can do that on your integrals before you put them into block2, as the following:

idx = np.array([0,3,1,4,2,5])
h1e = h1e[idx][:, idx]
g2e = g2e[idx][:, idx][:, :, idx][:, :, :, idx]

I think as the order matters, we can't look backwards to the string.

I agree that an arbitrary reordering of orbitals in the CSF is hard to implement. But if you simply reverse the orbital order and flip the coupling signs, it will not generate any CG factors. The resulting CSF string will match that from UGA-FCI, as long as you also use a reversed orbital order in your UGA-FCI inputs.

The following is a script (tested on your example system) to show that the CSF strings generated by the "reverse and exchange +/-" trick for CSFs from the singlet embedding MPS, can produce the correct energy expectation when interpreted as CSF strings in DRT/UGA. So it is ok to do this.

from pyblock2.driver.core import DMRGDriver, SymmetryTypes
from pyblock2._pyscf import mcscf as b2mcscf
from pyscf import gto, scf, mcscf, ao2mo
import numpy as np

mol = gto.M(atom="N 0 0 -3; N 0 0 3", basis="sto3g", symmetry="c1", verbose=0, spin=6)
mf = scf.RHF(mol).run(conv_tol=1E-14)

# get localized orbitals
lo_coeff, lo_occ, lo_energy = b2mcscf.get_uno(mf.to_uhf())
selected = b2mcscf.select_active_space(mol, lo_coeff, lo_occ, ao_labels=["N-2p"], atom_order=[0, 1])
lo_coeff[:, np.array(sorted(selected))] = lo_coeff[:, np.array(selected)]

# get active space integrals
ncas, n_elec, spin = 6, 6, 2
mc = mcscf.CASCI(mf, ncas, n_elec)
mc.mo_coeff = lo_coeff
h1e, ecore = mc.get_h1eff()
g2e = ao2mo.restore(1, mc.get_h2eff(), ncas)

g2e[np.abs(g2e) < 1E-12] = 0
h1e[np.abs(h1e) < 1E-12] = 0

# manual orbital reordering (optional)
# idx = np.argsort(np.array([0,2,5,1,3,4]))
# h1e = h1e[idx][:, idx]
# g2e = g2e[idx][:, idx][:, :, idx][:, :, :, idx]

bond_dims = [250] * 4 + [500] * 5
noises = [1e-4] * 4 + [1e-5] * 5 + [0]
thrds = [1e-20] * 9

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
                     thrds=thrds, cutoff=0, iprint=1)
print('DMRG energy = %20.15f' % energy)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)

# DRT-FCI
fd = driver.bw.b.FCIDUMP()
fd.initialize_su2(ncas, n_elec, spin, 0, ecore, h1e, g2e)
iqs = driver.bw.VectorSX([driver.bw.SX(n_elec, spin, 0)])
# Here 'True' indicates a reserved ordering of orbitals inside DRT
big = driver.bw.b.su2.DRTBigSite(iqs, True, driver.n_sites, driver.orb_sym, fd, 0)
big.drt = driver.bw.b.su2.DRT(big.drt.n_sites, big.drt.get_init_qs(), big.drt.orb_sym)
print('\n%r' % big.drt)

hamil_op = driver.bw.b.su2.OpElement(driver.bw.b.OpNames.H, driver.bw.b.SiteIndex(), driver.vacuum)
hmat = big.get_site_op(0, hamil_op)[0]
dci = np.zeros((hmat.n, ), dtype=float)
hmat.diag(driver.bw.b.Matrix(dci.reshape(-1, 1)))

def hop(kci):
    bci = np.zeros((hmat.n, ), dtype=kci.dtype)
    driver.bw.b.CSRMatrixFunctions.multiply(hmat, 0, driver.bw.b.Matrix(kci.reshape(-1, 1)), 0,
                                            driver.bw.b.Matrix(bci.reshape(-1, 1)), 1.0, 0.0)
    return bci

kci = np.zeros((hmat.m, ))
kci[big.drt >> (big.drt ^ 0)] = 1
eners = driver.bw.b.IterativeMatrixFunctions.davidson(hop, dci, [kci], True, conv_thrd=1E-20)
print('E(%10s) = %15.10f' % ("DRT-FCI", eners[0] + ecore))

for ix, (csf, val) in enumerate([(c, v) for c, v in sorted(zip(big.drt, kci), key=lambda x: -abs(x[1])) if abs(v) >= 0.05]):
    print('CSF', "%10d" % ix, csf, " = %20.15f" % val)

# Energy expectation from CSFs generated from MPS
mpsci = np.zeros((hmat.m, ))
for occ, val in zip(csfs, coeffs):
    csf = ''.join("0-+2"[x] for x in occ[::-1])
    mpsci[big.drt.index(csf)] = val

print('E(%10s) = %15.10f' % ("MPS-CSF", np.dot(mpsci, hop(mpsci)) + ecore))

A few notes:

1. Your code for orbital localization does not work for me and it cannot generate consistent orbitals when running the script multiple times. I simplified the orbital localization part. A small bug fix is required: 5c5b543
2. driver.bw.b.su2.DRTBigSite(iqs, True, ...) implicitly reverses the orbital ordering in the integrals for DRT. If you use False here then it will not reverse the orbital ordering in DRT but then the E[MPS-CSF] will not match the correct energy.
3. csf = ''.join("0-+2"[x] for x in occ[::-1]) implements the "reverse and exchange +/-" trick.

[bogdanoff]

Thank you for the answer! I would say then, that the UG spin coupling is actually done from the end of the list, such that $$(S_n), (S_n + S_{n-1}), (S_n + S_{n-1}+S_{n-2}), ... , (S_n + ... + S_0)$$ are good quantum numbers. It can be seen from a calculation of N3 molecule atom="N 0 0 -3; N 0 0 3; N 0 0 0", where 012 and 210 orders give rise to different WF's.

Maybe as a related question, is there a way to get coefficients of CSF's not with the highest weights, but from a finite list. Something like get_csf_coefficients(ket, list=my_csf_strings, iprint=1)?

[hczhai]

I would say then, that the UG spin coupling is actually done from the end of the list

Yes I agree.

is there a way to get coefficients of CSF's not with the highest weights, but from a finite list.

This is enabled from the new commit a31e999. You can do, for example,

csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.0, given_dets=['2-+0--', '0+-2--'], iprint=1)

to get only the coefficients from a given CSF set. The CSF notation should match the block2 notation, namely, with singlet embedding.

If cutoff != 0.0 and (given_dets == [] or given_dets is None), then it will fall back to the old behavior.
If cutoff != 0.0 and given_dets == [...], then the two constraints will apply simultaneously. This means, if the abs of the coefficient of a given CSF such as 2-+0-- is below the given cutoff, you will get the zero coefficient for this CSF as a result (but the actual coefficient may not be zero).

The complexity of the algorithm is $\mathcal{O}(KLD^2)$, depending on the MPS bond dimension $D$, the length of given_dets $L$, and the number of orbitals $K$.