Fix <S^2> for frozen_core References

adcc/OperatorIntegrals.py

@@ -255,14 +255,21 @@ def ssq_1p(self) -> OneParticleOperator:
     @timed_member_call("_import_timer")
     def ssq_2p(self) -> TwoParticleOperator:
         """Returns the two-particle part of the S^2 operator"""
-        # currently not implemented for CVS (but fc and fv should be fine)
+        # NOTE: the implementation might also work for CVS. But double check
+        # once CVS 2p densities are implemented before removing the
+        # exception.
         if "o2" in self.mospaces.subspaces:
             raise NotImplementedError("The 2-particle part of the SSq operator is "
                                       "only implemented for the occupied and "
                                       "virtual spaces, i.e., "
                                       "CVS is not supported yet.")
         # Intermediates
         # S^aa, S^ab and S^bb (spin projected overlap matrices)
+        # NOTE: For UHF, the diagonal spin blocks of the overlap matrix
+        # (in the MO basis) are diagonal.
+        # Only the off-diagonal spin blocks contain off-diagonal elements!
+        # (For RHF, the off-diagonal spin blocks are also diagonal!)
+        # -> apply the UHF assumptions below
         ovlp_bb: Tensor = self.overlap_ao
         coeff_map = {}
         for sp in self.mospaces.subspaces:
@@ -284,39 +291,51 @@ def ssq_2p(self) -> TwoParticleOperator:
         transform_operator_ao2mo_spin_projected(ovlp_bb, S_bb, coeff_map, "bb",
                                                 self._conv_tol)
 
-        # additional intermediate
+        # additional intermediate (is diagonal)
         S_aa_minus_bb = S_aa - S_bb
 
         op = TwoParticleOperator(self.mospaces, symmetry=OperatorSymmetry.HERMITIAN)
-        op.oooo = (
-            + 1.0 * einsum("ik,jl->ijkl", S_aa_minus_bb.oo,  S_aa_minus_bb.oo)
-            # exploit symmetry S_ab.oo = S_ba.oo.T
-            + 4.0 * einsum("ik,jl->ijkl", S_ab.oo.T, S_ab.oo)
-        ).antisymmetrise(2, 3)
-        op.ooov = (
-            + 2.0 * einsum("ik,ja->ijka", S_ab.oo.T, S_ab.ov)
-            + 2.0 * einsum("ik,ja->ijka", S_ab.oo, S_ab.vo.T)
-        ).antisymmetrise(0, 1)
-        op.oovv = (
-            + 2.0 * einsum("ia,jb->ijab", S_ab.vo.T, S_ab.ov)
-            + 2.0 * einsum("ia,jb->ijab", S_ab.ov, S_ab.vo.T)
-        ).antisymmetrise(2, 3)
-        op.ovov = (
-            + 0.5 * einsum("ij,ab->iajb", S_aa_minus_bb.oo, S_aa_minus_bb.vv)
-            + 1.0 * einsum("ij,ab->iajb", S_ab.oo.T, S_ab.vv)
-            + 1.0 * einsum("ij,ab->iajb", S_ab.oo, S_ab.vv.T)
-            - 1.0 * einsum("ib,aj->iajb", S_ab.vo.T, S_ab.vo)
-            - 1.0 * einsum("ib,aj->iajb", S_ab.ov, S_ab.ov.T)
-        )
-        op.ovvv = (
-            + 2.0 * einsum("ib,ac->iabc", S_ab.vo.T, S_ab.vv)
-            + 2.0 * einsum("ib,ac->iabc", S_ab.ov, S_ab.vv.T)
-        ).antisymmetrise(2, 3)
-        op.vvvv = (
-            + 1.0 * einsum("ac,bd->abcd", S_aa_minus_bb.vv, S_aa_minus_bb.vv)
-            # exploit symmetry S_ab.vv = S_ba.vv.T
-            + 4.0 * einsum("ac,bd->abcd", S_ab.vv, S_ab.vv.T)
-        ).antisymmetrise(2, 3)
+        for block in op.canonical_blocks:
+            p, q, r, s = split_spaces(block)
+            # + S^ba_pr S^ab_qs + S^ba_qs S^ab_pr
+            # - S^ba_ps S^ab_qr - S^ba_qr S^ab_ps
+            # first generate the 2 terms with a positive sign:
+            # since we only have S^ab: S^ba = (S^ab).transpose
+            res = einsum("rp,qs->pqrs", S_ab[r + p], S_ab[q + s])
+            if p == q and r == s:
+                res = 2.0 * res.symmetrise((0, 1), (2, 3))
+            else:
+                res += 1.0 * einsum("sq,pr->pqrs", S_ab[s + q], S_ab[p + r])
+            # + 0.5 * D_pr D_qs - 0.5 * D_ps D_qr
+            # D = S^aa - S^bb is diagonal in the MO basis
+            # -> only consider when the spaces match
+            if p == r and q == s:
+                res += 0.5 * einsum(
+                    "pr,qs->pqrs", S_aa_minus_bb[p + r], S_aa_minus_bb[q + s]
+                )
+            # then deal with the terms with negative sign:
+            # 2 S^ab terms and the single D term
+            if p == q and r == s:
+                # NOTE: We only need 1 of the antisymmetrisations to get the
+                # correct result. The second one is needed for the correct
+                # permutational symmetry in the symmetry object.
+                res = 2.0 * res.antisymmetrise(0, 1).antisymmetrise(2, 3)
+            elif p == q:
+                res = 2.0 * res.antisymmetrise(0, 1)
+            elif r == s:
+                res = 2.0 * res.antisymmetrise(2, 3)
+            else:
+                # always subtract the S^ab contributions
+                res += (
+                    - einsum("sp,qr->pqrs", S_ab[s + p], S_ab[q + r])
+                    - einsum("rq,ps->pqrs", S_ab[r + q], S_ab[p + s])
+                )
+                # conditionally subtract the D contribution
+                if p == s and q == r:
+                    res -= 0.5 * einsum(
+                        "ps,qr->pqrs", S_aa_minus_bb[p + s], S_aa_minus_bb[q + r]
+                    )
+            op[block] = res
         return op
 
     def _import_dipole_like_operator(

adcc/ReferenceState.py

@@ -24,7 +24,7 @@
 
 from .misc import cached_property
 from .Tensor import Tensor
-from .MoSpaces import MoSpaces
+from .MoSpaces import split_spaces, MoSpaces
 from .backends import import_scf_results
 from .OperatorIntegrals import OperatorIntegrals
 from .OneParticleDensity import OneParticleDensity
@@ -224,16 +224,27 @@ def density_2p(self):
         """
         Return the two-particle Hartree-Fock density in the MO basis
         """
+        # NOTE: running over canonical blocks here, to avoid incorporating
+        # implicitly knowledge about the definition of a canonical block in
+        # the implementation
         density = TwoParticleDensity(self.mospaces,
                                      symmetry=OperatorSymmetry.HERMITIAN)
+        # gamma_pqrs = delta_pr delta_qs - delta_ps delta_qr for occupied p,q,r,s
         for block in density.canonical_blocks:
-            sym = libadcc.make_symmetry_operator(self.mospaces, block,
-                                                 density.symmetry.to_str(), "1")
-            density[block] = Tensor(sym)
-        for ss in self.mospaces.subspaces_occupied:
-            density[ss + ss + ss + ss].set_mask("ijij", 1)
-            density[ss + ss + ss + ss].set_mask("ijji", -1)
-            density[ss + ss + ss + ss].set_mask("iijj", 0)
+            splitted = split_spaces(block)
+            # skip any blocks containing virtual subspaces
+            if any(sp in self.mospaces.subspaces_virtual for sp in splitted):
+                continue
+            p, q, r, s = splitted
+            if p == r and q == s:
+                density[block].set_mask("ijij", 1)
+            if p == s and q == r:
+                density[block].set_mask("ijji", -1)
+            # only set elements to zero if one of the checks were true and
+            # __getitem__ initialized a tensor. There should be no need
+            # to initialize all canonical blocks here.
+            if p == q and r == s and block in density.blocks_nonzero:
+                density[block].set_mask("iijj", 0)
         density.reference_state = self
         return density
 

adcc/tests/LazyMp_test.py

@@ -37,9 +37,11 @@
 test_cases = testcases.get_by_filename(
     "h2o_sto3g", "cn_sto3g", "h2o_def2tzvp", "cn_ccpvdz"
 )
+# Include all test cases (not only CVS) to verify the SSQ operator implementation
+# for the other (fc, fv, ...) cases
 unrestricted_cases = [
     (case.file_name, c) for case in test_cases if not case.restricted
-    for c in ["gen", "cvs"]
+    for c in case.cases
 ]
 small_cases = [
     (case.file_name, c) for case in test_cases if not case.only_full_mode

adcc/tests/ReferenceState_test.py

@@ -185,15 +185,25 @@ def test_list_fc_fv_cvs(self):
 systems = ["cn_sto3g", "cn_ccpvdz"]
 
 
+# <S^2> of the HF determinant does not depend on the orbital space partitioning
 @pytest.mark.parametrize("system", systems)
 def test_ssq_reference_state(system):
-    system: testcases.TestCase = testcases.get_by_filename(system).pop()
+    testcase: testcases.TestCase = testcases.get_by_filename(system).pop()
     # we need to run the scf calculation since we don't store the <S^2> values
-    scfres = run_hf("pyscf", system.xyz, system.basis, multiplicity=2)
-    hf = ReferenceState(scfres)
-
+    scfres = run_hf(
+        "pyscf", testcase.xyz, testcase.basis, multiplicity=testcase.multiplicity
+    )
     ref_ssq, _ = scfres.spin_square()
-
-    hf_ssq = hf.ssq
-
-    assert hf_ssq == pytest.approx(ref_ssq)
+    for case in testcase.cases:
+        core_orbitals = testcase.core_orbitals if "cvs" in case else None
+        frozen_core = testcase.frozen_core if "fc" in case else None
+        frozen_virtual = testcase.frozen_virtual if "fv" in case else None
+        hf = ReferenceState(
+            scfres, core_orbitals=core_orbitals, frozen_core=frozen_core,
+            frozen_virtual=frozen_virtual
+        )
+        if "cvs" in case:
+            with pytest.raises(NotImplementedError):
+                hf.ssq
+        else:
+            assert hf.ssq == pytest.approx(ref_ssq)

adcc/tests/data/SHA256SUMS

@@ -1,11 +1,11 @@
 4c7edc6d80f2ba3ec31629ad15b0be7c43cc03b08e413dd7f243e794f8b344cc  ch2nh2_sto3g_hfdata.hdf5
 e8078f935a662077b256189ddb50030001da98039c53911936384bfa662f7d7f  ch2nh2_sto3g_hfimport.hdf5
-4712a85a147f4ab5234747e37fb224d9ebce4aabb5f5940b20a72c3e0422326f  cn_ccpvdz_adcc_adc0.hdf5
-fd87a47a784d075685f8784d6b1abc5c2978c8ee689d07b2297797168bb76279  cn_ccpvdz_adcc_adc1.hdf5
-2a4321a040980548abef981250025f37170827d1feca9c8ef03792ebe42eff55  cn_ccpvdz_adcc_adc2.hdf5
-7de3622e32c4afa3aeca7eb1df634ffad9b7f9d81903dfc7e30dccddfb07e20a  cn_ccpvdz_adcc_adc2x.hdf5
-97df4a82ef459d5e314626699a5a652daa9d66374b94dfa1ba1e33924d34488b  cn_ccpvdz_adcc_adc3.hdf5
-aa904f679e5a611c5991fcf6a78f801a6e5c5dfdcd6abf34869849c61e3969ed  cn_ccpvdz_adcc_mpdata.hdf5
+ef02a3d689870970472c70e6fbb40959fbd31c3bc3079201d29998f9c7cb4a17  cn_ccpvdz_adcc_adc0.hdf5
+3c591f63e4b9b8cb2c761c7ab6ee1dfb3f9c711d4bae4ef6464d9a4b7cc8395c  cn_ccpvdz_adcc_adc1.hdf5
+f8d58ee3e10a6f0852d9175491bbc5ac60dd96dc083f539c7e75825cfb391e61  cn_ccpvdz_adcc_adc2.hdf5
+0fdd78e8d31850979504d69e69d38c0c6b961e00161d5a44f2ef9a0beeac0a4c  cn_ccpvdz_adcc_adc2x.hdf5
+4849ded0374ccb596a01f597f46c04f512a5c07f3e242360ac352d93ceceeea8  cn_ccpvdz_adcc_adc3.hdf5
+d3b874d004b4486ce4d80c348d920b7690795426b05cb0b48b51a24ecdda68be  cn_ccpvdz_adcc_mpdata.hdf5
 5a2b2b104a503127358b499ffb33caec50b0ffd8b6c1acec3c2ea2cfd20fd7f4  cn_ccpvdz_adcman_adc0.hdf5
 e4ee38107952c24884b26065117fd931f69fa8b375f7d1abee7a4a1bfbc49926  cn_ccpvdz_adcman_adc1.hdf5
 f55705f2bb847bdedb89cd32ecc08f91496878b68d0f02de80f853c1dc1d58d1  cn_ccpvdz_adcman_adc2.hdf5
@@ -14,13 +14,13 @@ fac03e2b1d8f6af5628bf437338b7fdcd1b12e3a1e3db99c28fc183186a0c4e8  cn_ccpvdz_adcm
 c666e9eee29e337dc41416e0d0fea0b0bc5f5c159aa79f04f302fd1e4b8a6667  cn_ccpvdz_adcman_mpdata.hdf5
 e7e66c864f1aa75b5954c8cbddd43a4fe724ee190a83acadeb6a2fbb1c8a1b37  cn_ccpvdz_hfdata.hdf5
 81d8067f9a015caa2ca0c7b0905d135e576686966fd5b4f49c599d92928bdf36  cn_ccpvdz_hfimport.hdf5
-75bf2ad89119e309f63ebc7168fecdd020445b1f66b955be951212f857179a00  cn_sto3g_adcc_adc0.hdf5
-3c9d0b4767347d33e57eedf296d9719df3f6ce7b57b6086205a6c82ce39a3979  cn_sto3g_adcc_adc1.hdf5
-c2c176afcf391ffba6c4e75667ddf40157d6d1b947f58c6fe2a718b2d67b69d7  cn_sto3g_adcc_adc2.hdf5
-00dce1133c101ac2df3d33f2d8b0d15817a3f52e766f87193c3c867c635b6f89  cn_sto3g_adcc_adc2x.hdf5
-ee6d0d31a537db01520c8a080bfab8db28714ccce94d2498dc5e66a202087dec  cn_sto3g_adcc_adc3.hdf5
-d021f840bb114afa7c72193f30712de39c4379b350642935f3dfb696b678d75d  cn_sto3g_adcc_adc4.hdf5
-729df817c0fda1a7c8f72c21dfc9eebbaf2bf67e3fd48ae9d6f3b239749ce7ba  cn_sto3g_adcc_mpdata.hdf5
+14f04545518246cecbb3d1cf9ee11de38b56dfab85c511504729e6c5662b920a  cn_sto3g_adcc_adc0.hdf5
+f5bbe4f047c46ed16efffc4fbd40632c6e842d51a1f05932b8f9c339903f0b1b  cn_sto3g_adcc_adc1.hdf5
+d91bbcbfd2eba802ccc972e2eb161e55217f01bd4791d84334db1caee1e0dd6d  cn_sto3g_adcc_adc2.hdf5
+4ba3fb8e3bf4cb567ddbe9a7dfa261116968b7d6a6afd8b8a0ddeac3d8ba381a  cn_sto3g_adcc_adc2x.hdf5
+ae81e599cede490704fb99165f7b8282a8666a1dd9905a52ef98ce2e7e17e89f  cn_sto3g_adcc_adc3.hdf5
+7f171d8789200d690a33715ed6ca244b7afd6c88d507f11a39acdd3205008826  cn_sto3g_adcc_adc4.hdf5
+047e64f35eaf65e3822d8fc8c7d56015cd147707d4b38052d653f724be41419f  cn_sto3g_adcc_mpdata.hdf5
 30e5801cc6e352b8c9c40bbb5d0f28e03b93c16db459b7d8288e2f716ae83593  cn_sto3g_adcman_adc0.hdf5
 63ff2f2c18ee2592cbb37636d2ce9436dbe5b41388c62f305e1f932c67a9e621  cn_sto3g_adcman_adc1.hdf5
 8d5a0d387dca1f6fee65d4204f9aadf4832bf482508325b477c6b498e5ace5fb  cn_sto3g_adcman_adc2.hdf5
@@ -71,13 +71,13 @@ cc1c329590e97a0ecdbaa78a3ca260fc14306f996ba71cc5b3832271a8616d7d  h2o_sto3g_adcm
 6a19c32cc28124a38dbf9ece8a6e5f7498c74056f69d3a251cbb405eca70abac  h2o_sto3g_adcman_mpdata.hdf5
 985a3a6c01f4af2f8b848e27569a9c5926599a50af88694b8cdfbf49e17d9899  h2o_sto3g_hfdata.hdf5
 63bd479812e03cc6411f9b6553abe7567dfd1a36fd6809e4c8d2d01101a80254  h2o_sto3g_hfimport.hdf5
-7f9e54d19a28b00ce8026b01e099a3e6177039c656d970e9ea21f4731afe3ce2  hf_631g_adcc_adc0.hdf5
-76d08d15d17b996e8dda2f0a9c3a24b44de825a2c2c3c8488d9bb3016171f69d  hf_631g_adcc_adc1.hdf5
-e2c544817690412751adb8866e5ab2dcc5708f2c6e601dc8a763ed6e8c9c9e03  hf_631g_adcc_adc2.hdf5
-9aa4416cfd7fbff97099abbffe7bcfee8ca71ee1a205128599003b1ec352e3b5  hf_631g_adcc_adc2x.hdf5
-00ddb83f979c0641fad719fc10dc5079924d821b5dc26fd01c2de1b47965db30  hf_631g_adcc_adc3.hdf5
-8a534d5e30f5b300d5a4e29f612482027162b60f7fb3c0adcd1c154a5b7b5a4a  hf_631g_adcc_adc4.hdf5
-39b86f1b89a424b9e7b980bedd3697e71930407bf17f34f6f64b695a8b6ea513  hf_631g_adcc_mpdata.hdf5
+d8b88f3a4f7d9b29080565357a3fb521041b075c08744f4970bf3a69fbd478b1  hf_631g_adcc_adc0.hdf5
+de7a5929d2e08416c1fb306c67360faa882439687f86434ce9ae59d1080e94de  hf_631g_adcc_adc1.hdf5
+2dcf58102fa08054becf9fcf23d0ff34fa4fbad4c4d3684a45b9a1001bd553ac  hf_631g_adcc_adc2.hdf5
+3754748a27ddaba70324f8f030c6a20956718ef82cb2218a0332529e381652ea  hf_631g_adcc_adc2x.hdf5
+5ed7cccd392d81e8402b55fed960c3a72b64bd7c4f6d4ccd2241541d832c0d1d  hf_631g_adcc_adc3.hdf5
+004a334465a5a1eec65e7d2d11c3e36856def7c52c8f5b87fbc5a23ec52b06c6  hf_631g_adcc_adc4.hdf5
+f6d28bf919cda4d53dd7321e9b317ce2573c8f50e6ee13a0165473cc391b9cf4  hf_631g_adcc_mpdata.hdf5
 841ffa376340fbc082123fc025c6ecee5a1733467856f27518fc6b76e5e5cc8c  hf_631g_adcman_adc0.hdf5
 be59746e12d4bc91a65d9767fa67eb854caa0ba0493ed49e7cefe34461590d3a  hf_631g_adcman_adc1.hdf5
 c7af0bd44374b2d87a84ad6a19bd7a85491b1c2cd29a6f1f7ef483cb65e9d5ee  hf_631g_adcman_adc2.hdf5

adcc/tests/data/update_testdata.sh

@@ -4,7 +4,7 @@
 trap "exit 1" TERM
 export SCRIPT_PID=$$
 
-SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/adcc_test_data/0.8.3/"
+SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/adcc_test_data/0.8.4/"
 
 SHAFILE="SHA256SUMS"