Improve CI vector sign convention consistency in csf_fci.

pyscf/csf_fci/csfstring.py

@@ -609,8 +609,14 @@ def _transform_det2csf (inparr, norb, neleca, nelecb, smult, reverse=False, csd_
             inparr[:,det_addrs] = 0
             continue
 
+        # sign convention of PySCF is A' B' |vac>
+        # sign convention of get_spin_evecs is A'(unpaired) B'(unpaired) C'(pairs) |vac>
+        ncomm = npair * (npair-1) // 2 # A'(paired) B'(paired) -> C'(pairs)
+        ncomm += npair * max (0, nelecb-npair) # A'B' -> AB'(unpaired) AB'(paired)
+        sgn = (-1) ** (ncomm % 2)
+
         t_ref = lib.logger.perf_counter ()
-        umat = np.asarray_chkfinite (get_spin_evecs (nspin, neleca, nelecb, smult))
+        umat = np.asarray_chkfinite (sgn * get_spin_evecs (nspin, neleca, nelecb, smult))
         size_umat = max (size_umat, umat.nbytes)
         ncsf_blk = ncsf # later on I can use this variable to implement a generator form of get_spin_evecs to save
         #                 memory when there are too many csfs
@@ -740,8 +746,14 @@ def ax_b (mat):
             dj = di + nconf*ndet
             mat_ij = mat[:,di:dj].reshape (nrow, nconf, ndet)
 
+            # sign convention of PySCF is A' B' |vac>
+            # sign convention of get_spin_evecs is A'(unpaired) B'(unpaired) C'(pairs) |vac>
+            ncomm = npair * (npair-1) // 2 # A'(paired) B'(paired) -> C'(pairs)
+            ncomm += npair * max (0, nelecb-npair) # A'B' -> AB'(unpaired) AB'(paired)
+            sgn = (-1) ** (ncomm % 2)
+
             nspin = neleca + nelecb - 2*npair
-            umat = np.asarray_chkfinite (get_spin_evecs (nspin, neleca, nelecb, smult))
+            umat = np.asarray_chkfinite (sgn * get_spin_evecs (nspin, neleca, nelecb, smult))
 
             outmat[:,ci:cj] = np.tensordot (mat_ij, umat, axes=1).reshape (nrow, ncsf*nconf, order='C')
 

pyscf/csf_fci/spin_op.py

@@ -0,0 +1,192 @@
+# The code in this file actually has nothing to do with CSFs
+# Ideally, this file would be a part of pyscf/fci/spin_op.py
+# However, I am sticking it here to avoid messy conflicts
+
+from functools import reduce
+import numpy as np
+from scipy import linalg
+from pyscf import lib
+from pyscf.fci import cistring
+from pyscf.fci.addons import _unpack_nelec
+
+def norm_sdown (smult, nelec):
+    s2 = smult-1
+    m2 = nelec[0] - nelec[1]
+    return np.sqrt ((s2-m2+2)*(s2+m2)/4)
+
+def norm_sup (smult, nelec):
+    s2 = smult-1
+    m2 = nelec[0] - nelec[1]
+    return np.sqrt ((s2-m2)*(s2+m2+2)/4)
+
+def contract_sladder(fcivec, norb, nelec, op=-1):
+    ''' Contract spin ladder operator S+ or S- with fcivec.
+        Changes neleca - nelecb without altering <S2>
+        Obtained by modifying pyscf.fci.spin_op.contract_ss
+    '''
+    neleca, nelecb = _unpack_nelec(nelec)
+    na = cistring.num_strings(norb, neleca)
+    nb = cistring.num_strings(norb, nelecb)
+    fcivec = fcivec.reshape(na,nb)
+    assert (op in (-1,1)), 'op = -1 or 1'
+    if ((op==-1 and (neleca==0 or nelecb==norb)) or
+        (op==1 and (neleca==norb or nelecb==0))): return np.zeros ((0,0))
+    # ^ Annihilate vacuum state ^
+
+    def gen_map(fstr_index, nelec, des=True):
+        a_index = fstr_index(range(norb), nelec)
+        amap = np.zeros((a_index.shape[0],norb,2), dtype=np.int32)
+        if des:
+            for k, tab in enumerate(a_index):
+                amap[k,tab[:,1]] = tab[:,2:]
+        else:
+            for k, tab in enumerate(a_index):
+                amap[k,tab[:,0]] = tab[:,2:]
+        return amap
+
+    if op==-1:
+        aindex = gen_map(cistring.gen_des_str_index, neleca)
+        bindex = gen_map(cistring.gen_cre_str_index, nelecb, False)
+    else:
+        aindex = gen_map(cistring.gen_cre_str_index, neleca, False)
+        bindex = gen_map(cistring.gen_des_str_index, nelecb)
+
+    ci1 = np.zeros((cistring.num_strings(norb,neleca+op),
+                   cistring.num_strings(norb,nelecb-op)))
+    nspin_comm = neleca-1 if op==-1 else neleca
+    spin_comm_fac = (-1) ** nspin_comm
+    for i in range(norb):
+        signa = aindex[:,i,1]
+        signb = bindex[:,i,1]
+        maska = np.where(signa!=0)[0]
+        maskb = np.where(signb!=0)[0]
+        addra = aindex[maska,i,0]
+        addrb = bindex[maskb,i,0]
+        citmp = lib.take_2d(fcivec, maska, maskb)
+        citmp *= signa[maska].reshape(-1,1)
+        citmp *= signb[maskb]
+        citmp *= spin_comm_fac
+        #: ci1[addra.reshape(-1,1),addrb] += citmp
+        lib.takebak_2d(ci1, citmp, addra, addrb)
+    return ci1
+
+def contract_sdown (ci, norb, nelec):
+    return contract_sladder (ci, norb, nelec, op=-1)
+
+def contract_sup (ci, norb, nelec):
+    return contract_sladder (ci, norb, nelec, op=1)
+
+def as_list (ci0):
+    if isinstance (ci0, np.ndarray) and ci0.ndim == 2:
+        ci1 = [ci0,]
+    else:
+        ci1 = list (ci0)
+    return ci1
+
+def like_ci0 (ci1, ci0):
+    if isinstance (ci0, np.ndarray) and ci0.ndim == 2:
+        ci1 = ci1[0]
+    else:
+        ci1 = np.stack (ci1, axis=0)
+    return ci1
+
+def mdown (ci0, norb, nelec, smult):
+    '''For a high-spin CI vector, rotate the spin axis away from the Z axis.
+
+    Args:
+        ci0 : ndarray or list of ndarrays
+            Contains one or multiple CI vectors in the Hilbert space of
+            (sum (nelec) + smult - 1) // 2 spin-up electrons and
+            (sum (nelec) - smult + 1) // 2 spin-down electrons (i.e., m=s)
+            with spin multiplicity given by smult
+        norb : integer
+            Number of orbitals
+        nelec : (integer, integer)
+            Number of spin-up and spin-down electrons in ci1
+        smult : integer
+            Spin multiplicity of ci0 and ci1
+
+    Returns:
+        ci1 : ndarray
+            Contains one or multiple CI vectors in the Hilbert space of nelec[0] spin-up electrons
+            and nelec[1] spin-down electrons occupying norb orbitals. If ci0 is spin-pure, then ci1
+            is the same state with a different Z axis (i.e., a different m microstate).
+    '''
+    ci1 = as_list (ci0)
+    neleca, nelecb = _unpack_nelec (nelec)
+    for i in range (len (ci1)):
+        for j in range (((smult-1)-(neleca-nelecb))//2,0,-1):
+            ci1[i] = contract_sdown (ci1[i], norb, (neleca+j,nelecb-j))
+            ci1[i] /= norm_sdown (smult, (neleca+j,nelecb-j))
+    return like_ci0 (ci1, ci0)
+
+def mup (ci0, norb, nelec, smult):
+    '''For a high-spin CI vector, rotate the spin axis towards the Z axis.
+
+    Args:
+        ci0 : ndarray or list of ndarrays
+            Contains one or multiple CI vectors in the Hilbert space of nelec[0] spin-up electrons
+            and nelec[1] spin-down electrons occupying norb orbitals.
+        norb : integer
+            Number of orbitals
+        nelec : (integer, integer)
+            Number of spin-up and spin-down electrons in ci0
+        smult : integer
+            Spin multiplicity of ci0 and ci1
+
+    Returns:
+        ci1 : ndarray
+            Contains one or multiple CI vectors in the Hilbert space of
+            (sum (nelec) + smult - 1) // 2 spin-up electrons and
+            (sum (nelec) - smult + 1) // 2 spin-down electrons (i.e., m=s).
+            If ci0 is spin-pure, then ci1 is the same state with a different Z axis (i.e., a
+            different m microstate).
+    '''
+    ci1 = as_list (ci0)
+    neleca, nelecb = _unpack_nelec (nelec)
+    for i in range (len (ci1)):
+        for j in range (((smult-1)-(neleca-nelecb))//2):
+            ci1[i] = contract_sup (ci1[i], norb, (neleca+j,nelecb-j))
+            ci1[i] /= norm_sup (smult, (neleca+j,nelecb-j))
+    return like_ci0 (ci1, ci0)
+
+if __name__ == '__main__':
+    import sys
+    import time
+    from pyscf.fci.direct_spin1 import contract_2e
+    from pyscf.fci.spin_op import spin_square0
+    t0, w0 = lib.logger.process_clock (), lib.logger.perf_counter ()
+    nelec, norb = (int (argv) for argv in sys.argv[1:])
+    nelec = (min (norb, nelec), nelec - min(norb, nelec))
+    smult = nelec[0]-nelec[1]+1
+    print ("Testing the spin ladder operators for a {}e, {}o s={} space".format (sum (nelec), norb, (smult-1)/2))
+    cishape = [cistring.num_strings (norb, ne) for ne in nelec]
+    np.random.seed(1)
+    ci = np.random.rand (*cishape)
+    eri = np.random.rand (norb,norb,norb,norb)
+    ci /= linalg.norm (ci)
+    print (" neleca nelecb ndeta ndetb {:>5s} {:>13s} {:>5s} {:>5s}".format ("cc","chc","ss","2s+1"))
+    def print_line (c, ne):
+        try:
+            ss, smult = spin_square0 (c, norb, ne)
+            hc = contract_2e (eri, c, norb, ne)
+            chc = c.conj ().ravel ().dot (hc.ravel ())
+        except AssertionError:
+            assert (any ([n<0 for n in ne]))
+            ss, smult, chc = (0.0, 1.0, 0.0)
+        cc = linalg.norm (c)
+        ndeta, ndetb = c.shape
+        print (" {:>6d} {:>6d} {:>5d} {:>5d} {:5.3f} {:13.9f} {:5.3f} {:5.3f}".format (
+            ne[0], ne[1], ndeta, ndetb, cc, chc, ss, smult))
+    print_line (ci, nelec)
+    for ndown in range (smult-1):
+        ci = contract_sdown (ci, norb, nelec)
+        nelec = (nelec[0]-1, nelec[1]+1)
+        print_line (ci, nelec)
+    for nup in range (smult):
+        ci = contract_sup (ci, norb, nelec)
+        nelec = (nelec[0]+1, nelec[1]-1)
+        print_line (ci, nelec)
+    print ("Time elapsed {} clock ; {} wall".format (lib.logger.process_clock () - t0,
+                                                     lib.logger.perf_counter () - w0))
+

pyscf/csf_fci/test/test_spin_op.py

@@ -0,0 +1,112 @@
+import numpy as np
+import unittest
+from scipy import linalg
+from pyscf.csf_fci.csfstring import CSFTransformer
+from pyscf.csf_fci.spin_op import contract_sdown, contract_sup, mdown, mup
+from pyscf.csf_fci.spin_op import norm_sdown, norm_sup
+from pyscf.fci.direct_spin1 import contract_2e
+from pyscf.fci.spin_op import spin_square0
+from pyscf.fci import cistring
+from itertools import product
+
+def get_cc_chc_smult (eri, c, norb, ne):
+    cc = c.conj ().ravel ().dot (c.ravel ())
+    chc = c.conj ().ravel ().dot (contract_2e (eri, c, norb, ne).ravel ())
+    ss, smult = spin_square0 (c, norb, ne)
+    return cc, chc, smult
+
+def opstring (op, s, ms):
+    parity = int(round(2*s))%2
+    if parity:
+        fmt = '{:.1f}'
+        s = round (s,1)
+        ms = round (ms,1)
+    else:
+        fmt = '{:d}'
+        s = int (round (s))
+        ms = int (round (ms))
+    fmt = 'S{}|' + fmt + ',' + fmt + '>'
+    return fmt.format (('-','z','+')[op+1], s, ms)
+
+class KnownValues(unittest.TestCase):
+
+    def test_ladder (self):
+        for norb, nelec_tot in product ((6,7), repeat=2):
+            nelec = (min (norb, nelec_tot), nelec_tot - min(norb, nelec_tot))
+            smult = nelec[0]-nelec[1]+1
+            cishape = [cistring.num_strings (norb, ne) for ne in nelec]
+            ci = np.random.rand (*cishape)
+            eri = np.random.rand (norb,norb,norb,norb)
+            ci /= linalg.norm (ci)
+            cc_ref, chc_ref, smult_ref = get_cc_chc_smult (eri, ci, norb, nelec)
+            self.assertAlmostEqual (cc_ref, 1.0, 9)
+            self.assertAlmostEqual (smult_ref, smult, 9)
+            s = (smult-1)/2
+            with self.subTest(norb=norb, nelec=nelec_tot, op=opstring (1,s,s)):
+                self.assertEqual (contract_sup (ci, norb, nelec).size, 0)
+            for ms in np.arange (s,-s,-1):
+                ci = contract_sdown (ci, norb, nelec)
+                ci /= norm_sdown (smult, nelec)
+                nelec = (nelec[0]-1, nelec[1]+1)
+                cc_test, chc_test, smult_test = get_cc_chc_smult (eri, ci, norb, nelec)
+                with self.subTest(norb=norb, nelec=nelec_tot, op=opstring(-1,s,ms)):
+                    self.assertAlmostEqual (cc_test, 1.0, 9)
+                    self.assertAlmostEqual (smult_test, smult, 9)
+                    self.assertAlmostEqual (chc_test, chc_ref, 9)
+            with self.subTest(norb=norb, nelec=nelec_tot, op=opstring (-1,s,-s)):
+                self.assertEqual (contract_sdown (ci, norb, nelec).size, 0)
+            for ms in np.arange (-s,s,1):
+                ci = contract_sup (ci, norb, nelec)
+                ci /= norm_sup (smult, nelec)
+                nelec = (nelec[0]+1, nelec[1]-1)
+                cc_test, chc_test, smult_test = get_cc_chc_smult (eri, ci, norb, nelec)
+                with self.subTest(norb=norb, nelec=nelec_tot, op=opstring(1,s,ms)):
+                    self.assertAlmostEqual (cc_test, 1.0, 9)
+                    self.assertAlmostEqual (smult_test, smult, 9)
+                    self.assertAlmostEqual (chc_test, chc_ref, 9)
+
+    def test_mdown_mup (self):
+        for norb, nelec_tot in product ((6,7), repeat=2):
+            nelec = (min (norb, nelec_tot), nelec_tot - min(norb, nelec_tot))
+            smult = nelec[0]-nelec[1]+1
+            cishape = [cistring.num_strings (norb, ne) for ne in nelec]
+            ci = np.random.rand (*cishape)
+            ci /= linalg.norm (ci)
+            s = (smult-1)/2
+            for itype, ci0 in enumerate ((ci, [ci,], ci[None,:,:])):
+                nelec1 = nelec
+                for ms in np.arange (s,-s,-1):
+                    with self.subTest (itype=itype, norb=norb, nelec_tot=nelec_tot, ms=ms):
+                        ci1 = mup (mdown (ci0, norb, nelec1, smult), norb, nelec1, smult)
+                        ovlp = np.ravel (ci0).dot (np.ravel (ci1))
+                        self.assertAlmostEqual (ovlp, 1.0, 9)
+                        nelec1 = (nelec1[0]-1,nelec1[1]+1)
+
+    def test_csf_sign (self):
+        rng = np.random.default_rng ()
+        for smult, ndocc, nvirt in product (range (1,8), range(3), range(3)):
+            if smult==1 and ndocc==0: continue
+            nelec = ((smult-1) + ndocc, ndocc)
+            norb = (smult-1) + ndocc + nvirt
+            trans0 = CSFTransformer (norb, nelec[0], nelec[1], smult)
+            ci_csf = 2 * rng.random (trans0.ncsf) - 1
+            ci_csf /= linalg.norm (ci_csf)
+            ci0 = trans0.vec_csf2det (ci_csf)
+            for spin in range (1-smult, smult-1, 2):
+                nelec = sum (nelec)
+                nelec = ((nelec + spin) // 2, (nelec - spin) // 2)
+                trans1 = CSFTransformer (norb, nelec[0], nelec[1], smult)
+                ci1 = mup (trans1.vec_csf2det (ci_csf).reshape (trans1.ndeta,trans1.ndetb),
+                           norb,
+                           nelec,
+                           smult)
+                with self.subTest (norb=norb, nelec=sum(nelec), smult=smult, spin=spin):
+                    self.assertAlmostEqual (ci1.ravel ().dot (ci0.ravel ()), 1.0, 8)
+
+
+if __name__ == "__main__":
+    print("Full Tests for fci.spin_op")
+    unittest.main()
+
+
+

pyscf/lib/csf/csfstring.c

@@ -105,11 +105,12 @@ void FCICSFmakecsf (double * umat, uint64_t * detstr, uint64_t * coupstr, int ns
         numerator = 1;
         denominator = 1;
         sgn = 1;
-        osgn = 1;
+        // Commute each spin-down electron past each spin-up electron
+        osgn = -track2MS;
         for (ispin = 1; ispin < nspin; ispin++){
             sup = (1ULL << ispin) & coupstr[icoup];
             msup = (1ULL << ispin) & detstr[idet];
-            if (msup){ track2MS++; osgn *= -1; } else { track2MS--; osgn *= -1; sgn *= osgn; }
+            if (msup){ track2MS++; osgn *= -1; } else { track2MS--; sgn *= osgn; }
             /* Clebsch-Gordan coefficient <j1,j2,m1,m2|J,M=m1+m2> (j2 = 1/2, m2 = +-1/2)
  *              = sgn * sqrt (num / denom)
  *              sgn = sgn(J2-j1)^delta(m2,+1/2)