From 79aa4e9a43d3e99f985613cf2b4a449b97e25431 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 11:57:00 -0500 Subject: [PATCH 1/7] add sign testing and spin_op.py --- pyscf/csf_fci/spin_op.py | 178 ++++++++++++++++++++++++++++++++ pyscf/csf_fci/test/test_sign.py | 34 ++++++ 2 files changed, 212 insertions(+) create mode 100644 pyscf/csf_fci/spin_op.py create mode 100644 pyscf/csf_fci/test/test_sign.py diff --git a/pyscf/csf_fci/spin_op.py b/pyscf/csf_fci/spin_op.py new file mode 100644 index 000000000..d23d6b37e --- /dev/null +++ b/pyscf/csf_fci/spin_op.py @@ -0,0 +1,178 @@ +# 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 contract_sladder(fcivec, norb, nelec, op=-1): + ''' Contract spin ladder operator S+ or S- with fcivec. + Changes neleca - nelecb without altering + Obtained by modifying pyscf.fci.spin_op.contract_ss + ''' + norm_ci0 = linalg.norm (fcivec) + 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) + norm_ci1 = linalg.norm (ci1) + np.finfo (float).tiny + ci1 *= norm_ci0 / norm_ci1 # ??? + 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)) + 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)) + 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 as e: + 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)) + diff --git a/pyscf/csf_fci/test/test_sign.py b/pyscf/csf_fci/test/test_sign.py new file mode 100644 index 000000000..a4fa59543 --- /dev/null +++ b/pyscf/csf_fci/test/test_sign.py @@ -0,0 +1,34 @@ +import unittest +import numpy as np +from scipy import linalg +from pyscf.csf_fci import spin_op +from pyscf.csf_fci.csfstring import CSFTransformer +import itertools + +class KnownValues(unittest.TestCase): + + def test_sign (self): + rng = np.random.default_rng () + for smult, ndocc, nvirt in itertools.product (range (1,6), 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 = spin_op.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 csfstring sign") + unittest.main () + From 323d3cc078899cdd9a224ae553ec602b7fc91f64 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 14:57:06 -0500 Subject: [PATCH 2/7] csf <-> det sign factor fix Apply the correct sign factor corresponding to arranging the electrons like A'(unpaired) B'(unpaired) C'(pairs) where ci' = ai' bi', compared to the PySCF convention of A' B'. --- pyscf/csf_fci/csfstring.py | 16 ++++++++++++++-- pyscf/csf_fci/test/test_sign.py | 5 +++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/pyscf/csf_fci/csfstring.py b/pyscf/csf_fci/csfstring.py index ef51238b0..bf7932067 100644 --- a/pyscf/csf_fci/csfstring.py +++ b/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') diff --git a/pyscf/csf_fci/test/test_sign.py b/pyscf/csf_fci/test/test_sign.py index a4fa59543..746bc15f2 100644 --- a/pyscf/csf_fci/test/test_sign.py +++ b/pyscf/csf_fci/test/test_sign.py @@ -9,7 +9,8 @@ class KnownValues(unittest.TestCase): def test_sign (self): rng = np.random.default_rng () - for smult, ndocc, nvirt in itertools.product (range (1,6), range(3), range(3)): + for smult, ndocc, nvirt in itertools.product (range (1,6), range(4), range(3)): + if nvirt > 0: continue if smult==1 and ndocc==0: continue nelec = ((smult-1) + ndocc, ndocc) norb = (smult-1) + ndocc + nvirt @@ -27,7 +28,7 @@ def test_sign (self): 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 csfstring sign") unittest.main () From afb1508f45190fe320948dc80e70438f772e478d Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 15:17:58 -0500 Subject: [PATCH 3/7] get_spin_evecs operator order sign fix Correct the sign convention of csfstring.c:FCICSFmakecsf. It had been like ai' a(i+1)' a(i+2)' ... b(j)' b(j+1)' b(j+2)'. It is now correctly ai' a(i+1)' a(i+2)' ... b(j+2)' b(j+1)' b(j)'. Also add test_spin_op.py from mrh. --- pyscf/csf_fci/test/test_sign.py | 35 --------- pyscf/csf_fci/test/test_spin_op.py | 110 +++++++++++++++++++++++++++++ pyscf/lib/csf/csfstring.c | 5 +- 3 files changed, 113 insertions(+), 37 deletions(-) delete mode 100644 pyscf/csf_fci/test/test_sign.py create mode 100644 pyscf/csf_fci/test/test_spin_op.py diff --git a/pyscf/csf_fci/test/test_sign.py b/pyscf/csf_fci/test/test_sign.py deleted file mode 100644 index 746bc15f2..000000000 --- a/pyscf/csf_fci/test/test_sign.py +++ /dev/null @@ -1,35 +0,0 @@ -import unittest -import numpy as np -from scipy import linalg -from pyscf.csf_fci import spin_op -from pyscf.csf_fci.csfstring import CSFTransformer -import itertools - -class KnownValues(unittest.TestCase): - - def test_sign (self): - rng = np.random.default_rng () - for smult, ndocc, nvirt in itertools.product (range (1,6), range(4), range(3)): - if nvirt > 0: continue - 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 = spin_op.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 csfstring sign") - unittest.main () - diff --git a/pyscf/csf_fci/test/test_spin_op.py b/pyscf/csf_fci/test/test_spin_op.py new file mode 100644 index 000000000..55221feb5 --- /dev/null +++ b/pyscf/csf_fci/test/test_spin_op.py @@ -0,0 +1,110 @@ +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.fci.direct_spin1 import contract_2e +from pyscf.fci.spin_op import spin_square0 +from pyscf.fci import cistring +from itertools import product + +np.random.seed(1) +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) + 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) + 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() + + + diff --git a/pyscf/lib/csf/csfstring.c b/pyscf/lib/csf/csfstring.c index 4b2e3df9c..a5e61f45c 100644 --- a/pyscf/lib/csf/csfstring.c +++ b/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 (j2 = 1/2, m2 = +-1/2) * = sgn * sqrt (num / denom) * sgn = sgn(J2-j1)^delta(m2,+1/2) From 14268c5785e58a4abdfb8e9c6c9aa14ae3d28b8e Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 17:27:41 -0500 Subject: [PATCH 4/7] Norm handling in spin_op.py Instead of just fixing output to the same norm as input, use the theoretically derived normalization factors. --- pyscf/csf_fci/spin_op.py | 14 ++++++++++++-- pyscf/csf_fci/test/test_spin_op.py | 3 +++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/pyscf/csf_fci/spin_op.py b/pyscf/csf_fci/spin_op.py index d23d6b37e..bf7119220 100644 --- a/pyscf/csf_fci/spin_op.py +++ b/pyscf/csf_fci/spin_op.py @@ -9,6 +9,16 @@ 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 @@ -59,8 +69,6 @@ def gen_map(fstr_index, nelec, des=True): citmp *= spin_comm_fac #: ci1[addra.reshape(-1,1),addrb] += citmp lib.takebak_2d(ci1, citmp, addra, addrb) - norm_ci1 = linalg.norm (ci1) + np.finfo (float).tiny - ci1 *= norm_ci0 / norm_ci1 # ??? return ci1 def contract_sdown (ci, norb, nelec): return contract_sladder (ci, norb, nelec, op=-1) @@ -107,6 +115,7 @@ def mdown (ci0, norb, nelec, smult): 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): @@ -136,6 +145,7 @@ def mup (ci0, norb, nelec, smult): 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__': diff --git a/pyscf/csf_fci/test/test_spin_op.py b/pyscf/csf_fci/test/test_spin_op.py index 55221feb5..06f7f50df 100644 --- a/pyscf/csf_fci/test/test_spin_op.py +++ b/pyscf/csf_fci/test/test_spin_op.py @@ -3,6 +3,7 @@ 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 @@ -46,6 +47,7 @@ def test_ladder (self): 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)): @@ -56,6 +58,7 @@ def test_ladder (self): 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)): From a7f78b51667de0306916a726582b1efe1b6dcb8d Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 17:35:42 -0500 Subject: [PATCH 5/7] linting --- pyscf/csf_fci/spin_op.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/pyscf/csf_fci/spin_op.py b/pyscf/csf_fci/spin_op.py index bf7119220..90577c79f 100644 --- a/pyscf/csf_fci/spin_op.py +++ b/pyscf/csf_fci/spin_op.py @@ -24,7 +24,6 @@ def contract_sladder(fcivec, norb, nelec, op=-1): Changes neleca - nelecb without altering Obtained by modifying pyscf.fci.spin_op.contract_ss ''' - norm_ci0 = linalg.norm (fcivec) neleca, nelecb = _unpack_nelec(nelec) na = cistring.num_strings(norb, neleca) nb = cistring.num_strings(norb, nelecb) @@ -71,8 +70,11 @@ def gen_map(fstr_index, nelec, des=True): 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 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: @@ -169,12 +171,13 @@ def print_line (c, ne): ss, smult = spin_square0 (c, norb, ne) hc = contract_2e (eri, c, norb, ne) chc = c.conj ().ravel ().dot (hc.ravel ()) - except AssertionError as e: + 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 (" {:>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) @@ -184,5 +187,6 @@ def print_line (c, ne): 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)) + print ("Time elapsed {} clock ; {} wall".format (lib.logger.process_clock () - t0, + lib.logger.perf_counter () - w0)) From 1be4f63fb12c3f8acbf18344e6d8149ea058b793 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 13 Oct 2025 17:37:32 -0500 Subject: [PATCH 6/7] linting --- pyscf/csf_fci/spin_op.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyscf/csf_fci/spin_op.py b/pyscf/csf_fci/spin_op.py index 90577c79f..fe8bc7f9b 100644 --- a/pyscf/csf_fci/spin_op.py +++ b/pyscf/csf_fci/spin_op.py @@ -187,6 +187,6 @@ def print_line (c, ne): 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, + print ("Time elapsed {} clock ; {} wall".format (lib.logger.process_clock () - t0, lib.logger.perf_counter () - w0)) From b0c04958438925dac45162d17e5e2f2acc63afe9 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Tue, 14 Oct 2025 11:58:46 -0500 Subject: [PATCH 7/7] remove unnecessary seed Left over from debugging, oops --- pyscf/csf_fci/test/test_spin_op.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyscf/csf_fci/test/test_spin_op.py b/pyscf/csf_fci/test/test_spin_op.py index 06f7f50df..183b2cc57 100644 --- a/pyscf/csf_fci/test/test_spin_op.py +++ b/pyscf/csf_fci/test/test_spin_op.py @@ -9,7 +9,6 @@ from pyscf.fci import cistring from itertools import product -np.random.seed(1) 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 ())