diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index 69d06b61..c5df9ee3 100644 --- a/adcc/OperatorIntegrals.py +++ b/adcc/OperatorIntegrals.py @@ -255,7 +255,9 @@ 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 " @@ -263,6 +265,11 @@ def ssq_2p(self) -> TwoParticleOperator: "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( diff --git a/adcc/ReferenceState.py b/adcc/ReferenceState.py index d6d718c1..c0ad0ad5 100644 --- a/adcc/ReferenceState.py +++ b/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_ijkl = delta_ik delta_jl - delta_il delta_jk 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 + i, j, k, l = splitted # noqa: E741 + if i == k and j == l: + density[block].set_mask("ijij", 1) + if i == l and j == k: + 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 i == j and k == l and block in density.blocks_nonzero: + density[block].set_mask("iijj", 0) density.reference_state = self return density diff --git a/adcc/tests/LazyMp_test.py b/adcc/tests/LazyMp_test.py index b9d10281..447f5955 100644 --- a/adcc/tests/LazyMp_test.py +++ b/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 diff --git a/adcc/tests/ReferenceState_test.py b/adcc/tests/ReferenceState_test.py index e0520537..7ccac3d8 100644 --- a/adcc/tests/ReferenceState_test.py +++ b/adcc/tests/ReferenceState_test.py @@ -185,15 +185,25 @@ def test_list_fc_fv_cvs(self): systems = ["cn_sto3g", "cn_ccpvdz"] +# 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 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) diff --git a/adcc/tests/data/SHA256SUMS b/adcc/tests/data/SHA256SUMS index 3e9cfe6b..c7699f52 100644 --- a/adcc/tests/data/SHA256SUMS +++ b/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 diff --git a/adcc/tests/data/update_testdata.sh b/adcc/tests/data/update_testdata.sh index dc688dcc..46fdabbc 100755 --- a/adcc/tests/data/update_testdata.sh +++ b/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"