diff --git a/cucc/README.md b/cucc/README.md index 3222f53..36033e8 100644 --- a/cucc/README.md +++ b/cucc/README.md @@ -3,3 +3,9 @@ This package follows the interfaces in PySCF packages, and some codes are adapted from it to keep the same interface. The examples are in the `byteqc/cucc/test` folder. The `buffer.py` and `culib.py` modules enable automatic backend determination. The arrays in `cucc` can be stored on the GPU, CPU, or disk, all with the same interface. Users and developers can work with these arrays without having to consider the underlying backend. + +`cucc.CCSD` accepts PySCF UHF references and returns unrestricted amplitudes +in the usual `(t1a, t1b)` and `(t2aa, t2ab, t2bb)` form. The UCCSD and +UCCSD(T) contractions run on the GPU in an antisymmetrized spin-orbital +representation. The initial implementation keeps transformed integrals +in-core and checks their footprint against `gpulim * mem_ratio` before upload. diff --git a/cucc/__init__.py b/cucc/__init__.py index 18b632f..0fe084e 100644 --- a/cucc/__init__.py +++ b/cucc/__init__.py @@ -13,23 +13,30 @@ # See the License for the specific language governing permissions and # limitations under the License. -from pyscf import scf -from byteqc.cucc import ccsd -from byteqc.cucc import dfccsd from pyscf.lib import param def CCSD(mf, frozen=None, mo_coeff=None, mo_occ=None, gpulim=None, cpulim=None, pool=None, path=param.TMPDIR, mem_ratio=0.65): - if isinstance(mf, scf.uhf.UHF) or isinstance(mf, scf.ghf.GHF): - AssertionError('Not implement') - else: - return RCCSD(mf, frozen, mo_coeff, mo_occ, gpulim, cpulim, + if mf.istype('UHF'): + return UCCSD(mf, frozen, mo_coeff, mo_occ, gpulim, cpulim, pool=pool, path=path, mem_ratio=mem_ratio) + if mf.istype('GHF'): + raise NotImplementedError('GHF-CCSD is not implemented') + return RCCSD(mf, frozen, mo_coeff, mo_occ, gpulim, cpulim, + pool=pool, path=path, mem_ratio=mem_ratio) + + +def UCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None, gpulim=None, + cpulim=None, pool=None, path=param.TMPDIR, mem_ratio=0.65): + from byteqc.cucc import uccsd + return uccsd.UCCSD(mf, frozen, mo_coeff, mo_occ, gpulim, cpulim, + pool=pool, path=path, mem_ratio=mem_ratio) def RCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None, gpulim=None, cpulim=None, pool=None, path=param.TMPDIR, mem_ratio=0.65): + from byteqc.cucc import ccsd, dfccsd if getattr(mf, 'with_df', None): return dfccsd.RCCSD(mf, frozen, mo_coeff, mo_occ, gpulim, cpulim, pool=pool, path=path, mem_ratio=mem_ratio) diff --git a/cucc/test/test_uccsd.py b/cucc/test/test_uccsd.py new file mode 100644 index 0000000..116b77c --- /dev/null +++ b/cucc/test/test_uccsd.py @@ -0,0 +1,60 @@ +# Copyright (c) 2024 Bytedance Ltd. and/or its affiliates +# This file is part of ByteQC. +# +# Licensed under the Apache License, Version 2.0 (the "License") +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https: // www.apache.org/licenses/LICENSE-2.0 + +import unittest + +import cupy +import numpy +from byteqc import cucc +from pyscf import cc, gto, scf + + +class UCCSDTest(unittest.TestCase): + + @classmethod + def setUpClass(cls): + mol = gto.M( + atom='Li 0 0 0; H 0 0 1.6', + basis='sto-3g', + charge=1, + spin=1, + verbose=0, + ) + cls.mf = scf.UHF(mol).run(conv_tol=1e-12) + + def test_uccsd_and_triples_match_pyscf(self): + ref = cc.UCCSD(self.mf).run(conv_tol=1e-10, + conv_tol_normt=1e-8) + ref_et = ref.ccsd_t() + + gpu = cucc.CCSD(self.mf) + gpu.conv_tol = 1e-10 + gpu.conv_tol_normt = 1e-8 + gpu.kernel() + gpu_et = gpu.ccsd_t() + + self.assertTrue(gpu.converged) + self.assertAlmostEqual(gpu.e_corr, ref.e_corr, 9) + self.assertAlmostEqual(gpu_et, ref_et, 9) + self.assertEqual(len(gpu.t1), 2) + self.assertEqual(len(gpu.t2), 3) + for actual, expected in zip(gpu.t1, ref.t1): + numpy.testing.assert_allclose(cupy.asnumpy(actual), expected, + atol=1e-8, rtol=1e-7) + for actual, expected in zip(gpu.t2, ref.t2): + numpy.testing.assert_allclose(cupy.asnumpy(actual), expected, + atol=1e-8, rtol=1e-7) + + def test_ghf_failure_is_explicit(self): + with self.assertRaisesRegex(NotImplementedError, 'GHF-CCSD'): + cucc.CCSD(self.mf.to_ghf()) + + +if __name__ == '__main__': + unittest.main() diff --git a/cucc/uccsd.py b/cucc/uccsd.py new file mode 100644 index 0000000..e628dc8 --- /dev/null +++ b/cucc/uccsd.py @@ -0,0 +1,458 @@ +# Copyright (c) 2024 Bytedance Ltd. and/or its affiliates +# This file is part of ByteQC. +# +# Licensed under the Apache License, Version 2.0 (the "License") +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https: // www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# +# ByteQC includes code adapted from PySCF (https://github.com/pyscf/pyscf), +# which is licensed under the Apache License 2.0. The original copyright: +# Copyright 2014-2021 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""GPU-accelerated unrestricted CCSD and CCSD(T). + +The public amplitudes use PySCF's spatial-orbital UCCSD convention: +``t1 == (t1a, t1b)`` and ``t2 == (t2aa, t2ab, t2bb)``. Internally the +equations are evaluated in an antisymmetrized spin-orbital representation. +This keeps one contraction path for all spin blocks and maps naturally to +CuPy while preserving the established unrestricted interface. +""" + +from functools import reduce + +import cupy +import numpy +from byteqc.cucc.buffer import BufferPool +from pyscf import ao2mo +from pyscf import lib as pyscf_lib +from pyscf.cc import addons +from pyscf.cc import gccsd as pyscf_gccsd +from pyscf.cc import uccsd as pyscf_uccsd +from pyscf.lib import logger, param + + +def _einsum(subscripts, *operands): + return cupy.einsum(subscripts, *operands, optimize=True) + + +def _make_tau(t2, t1a, t1b, fac=1): + t1t1 = _einsum('ia,jb->ijab', fac * .5 * t1a, t1b) + t1t1 = t1t1 - t1t1.transpose(1, 0, 2, 3) + tau = t1t1 - t1t1.transpose(0, 1, 3, 2) + tau += t2 + return tau + + +def _cc_fvv(t1, t2, eris): + nocc = t1.shape[0] + fov = eris.fock[:nocc, nocc:] + fvv = eris.fock[nocc:, nocc:] + tau = _make_tau(t2, t1, t1, fac=.5) + out = fvv - .5 * _einsum('me,ma->ae', fov, t1) + out += _einsum('mf,amef->ae', t1, eris.ovvv.transpose(1, 0, 3, 2)) + out -= .5 * _einsum('mnaf,mnef->ae', tau, eris.oovv) + return out + + +def _cc_foo(t1, t2, eris): + nocc = t1.shape[0] + fov = eris.fock[:nocc, nocc:] + foo = eris.fock[:nocc, :nocc] + tau = _make_tau(t2, t1, t1, fac=.5) + out = foo + .5 * _einsum('me,ie->mi', fov, t1) + out += _einsum('ne,mnie->mi', t1, eris.ooov) + out += .5 * _einsum('inef,mnef->mi', tau, eris.oovv) + return out + + +def _cc_fov(t1, t2, eris): + nocc = t1.shape[0] + return (eris.fock[:nocc, nocc:] + + _einsum('nf,mnef->me', t1, eris.oovv)) + + +def _cc_woooo(t1, t2, eris): + tau = _make_tau(t2, t1, t1) + tmp = _einsum('je,mnie->mnij', t1, eris.ooov) + out = eris.oooo + tmp - tmp.transpose(0, 1, 3, 2) + out += .25 * _einsum('ijef,mnef->mnij', tau, eris.oovv) + return out + + +def _cc_wvvvv(t1, t2, eris): + tau = _make_tau(t2, t1, t1) + tmp = _einsum('mb,mafe->bafe', t1, eris.ovvv) + out = eris.vvvv - tmp + tmp.transpose(1, 0, 2, 3) + out += _einsum('mnab,mnef->abef', tau, .25 * eris.oovv) + return out + + +def _cc_wovvo(t1, t2, eris): + eris_ovvo = -eris.ovov.transpose(0, 1, 3, 2) + eris_oovo = -eris.ooov.transpose(0, 1, 3, 2) + out = _einsum('jf,mbef->mbej', t1, eris.ovvv) + out -= _einsum('nb,mnej->mbej', t1, eris_oovo) + out -= .5 * _einsum('jnfb,mnef->mbej', t2, eris.oovv) + out -= _einsum('jf,nb,mnef->mbej', t1, t1, eris.oovv) + out += eris_ovvo + return out + + +def update_amps(cc, t1, t2, eris): + """Evaluate the spin-orbital CCSD amplitude equations on the GPU.""" + nocc, nvir = t1.shape + fov = eris.fock[:nocc, nocc:] + mo_e_o = eris.mo_energy[:nocc] + mo_e_v = eris.mo_energy[nocc:] + cc.level_shift + + tau = _make_tau(t2, t1, t1) + fvv = _cc_fvv(t1, t2, eris) + foo = _cc_foo(t1, t2, eris) + fov_i = _cc_fov(t1, t2, eris) + woooo = _cc_woooo(t1, t2, eris) + wvvvv = _cc_wvvvv(t1, t2, eris) + wovvo = _cc_wovvo(t1, t2, eris) + + fvv[cupy.diag_indices(nvir)] -= mo_e_v + foo[cupy.diag_indices(nocc)] -= mo_e_o + + t1new = _einsum('ie,ae->ia', t1, fvv) + t1new -= _einsum('ma,mi->ia', t1, foo) + t1new += _einsum('imae,me->ia', t2, fov_i) + t1new -= _einsum('nf,naif->ia', t1, eris.ovov) + t1new -= .5 * _einsum('imef,maef->ia', t2, eris.ovvv) + t1new -= .5 * _einsum('mnae,mnie->ia', t2, eris.ooov) + t1new += fov.conj() + + ftmp = fvv - .5 * _einsum('mb,me->be', t1, fov_i) + tmp = _einsum('ijae,be->ijab', t2, ftmp) + t2new = tmp - tmp.transpose(0, 1, 3, 2) + ftmp = foo + .5 * _einsum('je,me->mj', t1, fov_i) + tmp = _einsum('imab,mj->ijab', t2, ftmp) + t2new -= tmp - tmp.transpose(1, 0, 2, 3) + t2new += eris.oovv.conj() + t2new += .5 * _einsum('mnab,mnij->ijab', tau, woooo) + t2new += .5 * _einsum('ijef,abef->ijab', tau, wvvvv) + tmp = _einsum('imae,mbej->ijab', t2, wovvo) + tmp += _einsum('ie,ma,mbje->ijab', t1, t1, eris.ovov) + tmp = tmp - tmp.transpose(1, 0, 2, 3) + tmp = tmp - tmp.transpose(0, 1, 3, 2) + t2new += tmp + tmp = _einsum('ie,jeba->ijab', t1, eris.ovvv.conj()) + t2new += tmp - tmp.transpose(1, 0, 2, 3) + tmp = _einsum('ma,ijmb->ijab', t1, eris.ooov.conj()) + t2new -= tmp - tmp.transpose(0, 1, 3, 2) + + eia = mo_e_o[:, None] - mo_e_v + t1new /= eia + t2new /= eia[:, None, :, None] + eia[None, :, None, :] + return t1new, t2new + + +def energy(cc, t1, t2, eris): + nocc = t1.shape[0] + out = _einsum('ia,ia', eris.fock[:nocc, nocc:], t1) + out += .25 * _einsum('ijab,ijab', t2, eris.oovv) + out += .5 * _einsum('ia,jb,ijab', t1, t1, eris.oovv) + out = out.item() + if abs(out.imag) > 1e-4: + logger.warn(cc, 'Non-zero imaginary part found in UCCSD energy %s', out) + return out.real + + +def ccsd_t(cc, eris, t1, t2): + """Evaluate the spin-orbital (T) correction on the GPU.""" + nocc, nvir = t1.shape + bcei = eris.ovvv.conj().transpose(3, 2, 1, 0) + majk = eris.ooov.conj().transpose(2, 3, 0, 1) + bcjk = eris.oovv.conj().transpose(2, 3, 0, 1) + fvo = eris.fock[nocc:, :nocc] + eijk = (eris.mo_energy[:nocc, None, None] + + eris.mo_energy[None, :nocc, None] + + eris.mo_energy[None, None, :nocc]) + ev = eris.mo_energy[nocc:] + t2t = t2.transpose(2, 3, 0, 1) + t1t = t1.T + + def get_wv(a, b, c): + w = _einsum('ejk,ei->ijk', t2t[a], bcei[b, c]) + w -= _einsum('im,mjk->ijk', t2t[b, c], majk[:, a]) + v = _einsum('i,jk->ijk', t1t[a], bcjk[b, c]) + v += _einsum('i,jk->ijk', fvo[a], t2t[b, c]) + v += w + w = w + w.transpose(2, 0, 1) + w.transpose(1, 2, 0) + return w, v + + et = cupy.asarray(0, dtype=cupy.result_type(t1, t2, eris.oovv)) + for a in range(nvir): + for b in range(a): + for c in range(b): + wabc, vabc = get_wv(a, b, c) + wcab, vcab = get_wv(c, a, b) + wbac, vbac = get_wv(b, a, c) + w = wabc + wcab - wbac + v = vabc + vcab - vbac + w /= eijk - ev[a] - ev[b] - ev[c] + et += _einsum('ijk,ijk', w, v.conj()) + return (et * .5).real.item() + + +class _PhysicistsERIs: + """Antisymmetrized spin-orbital integrals resident on the GPU.""" + + def __init__(self): + self.mo_coeff = None + self.mo_energy = None + self.orbspin = None + self.nocc = None + self.fock = None + self.oooo = None + self.ooov = None + self.oovv = None + self.ovov = None + self.ovvv = None + self.vvvv = None + + +class UCCSD(pyscf_uccsd.UCCSD): + """Unrestricted CCSD with GPU spin-orbital contractions. + + The first implementation is in-core on the GPU. ``gpulim`` and + ``mem_ratio`` are enforced before integral upload so an oversized problem + fails with a useful message instead of silently spilling or overcommitting + device memory. + """ + + update_amps = update_amps + + def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None, + gpulim=None, cpulim=None, pool=None, path=param.TMPDIR, + mem_ratio=.65): + super().__init__(mf, frozen=frozen, mo_coeff=mo_coeff, mo_occ=mo_occ) + if not 0 < mem_ratio <= 1: + raise ValueError('mem_ratio must be in the interval (0, 1]') + self.pool = (pool if pool is not None else + BufferPool(gpulim, cpulim, path, verbose=self.verbose)) + self.mem_ratio = mem_ratio + self._t1_spin = None + self._t2_spin = None + self._keys.update(['pool', 'mem_ratio']) + + @property + def _spin_slices(self): + nocca, noccb = self.nocc + nmoa, nmob = self.nmo + nvira, nvirb = nmoa - nocca, nmob - noccb + return (slice(0, nocca), slice(nocca, nocca + noccb), + slice(0, nvira), slice(nvira, nvira + nvirb)) + + def _spin_to_spatial(self, t1, t2): + oa, ob, va, vb = self._spin_slices + return ((t1[oa, va], t1[ob, vb]), + (t2[oa, oa, va, va], t2[oa, ob, va, vb], + t2[ob, ob, vb, vb])) + + def _spatial_to_spin(self, t1, t2): + if getattr(t1, 'ndim', None) == 2: + return cupy.asarray(t1), cupy.asarray(t2) + orbspin = numpy.hstack(( + numpy.zeros(self.nocc[0], dtype=int), + numpy.ones(self.nocc[1], dtype=int), + numpy.zeros(self.nmo[0] - self.nocc[0], dtype=int), + numpy.ones(self.nmo[1] - self.nocc[1], dtype=int))) + t1 = tuple(cupy.asnumpy(x) if isinstance(x, cupy.ndarray) else x + for x in t1) + t2 = tuple(cupy.asnumpy(x) if isinstance(x, cupy.ndarray) else x + for x in t2) + return (cupy.asarray(addons.spatial2spin(t1, orbspin)), + cupy.asarray(addons.spatial2spin(t2, orbspin))) + + def ao2mo(self, mo_coeff=None): + if mo_coeff is None: + mo_coeff = self.mo_coeff + maska, maskb = self.get_frozen_mask() + occa = numpy.asarray(self.mo_occ[0]) > 0 + occb = numpy.asarray(self.mo_occ[1]) > 0 + groups = ((0, numpy.where(maska & occa)[0]), + (1, numpy.where(maskb & occb)[0]), + (0, numpy.where(maska & ~occa)[0]), + (1, numpy.where(maskb & ~occb)[0])) + coeff = [] + orbspin = [] + for spin, idx in groups: + coeff.append(numpy.asarray(mo_coeff[spin])[:, idx]) + orbspin.extend([spin] * len(idx)) + coeff = numpy.hstack(coeff) + orbspin = numpy.asarray(orbspin, dtype=int) + nmo = coeff.shape[1] + nocc = sum(self.nocc) + + dm = self._scf.make_rdm1(self.mo_coeff, self.mo_occ) + vhf = self._scf.get_veff(self.mol, dm) + fockao = self._scf.get_fock(vhf=vhf, dm=dm) + fock = numpy.zeros((nmo, nmo), dtype=numpy.result_type(*fockao)) + starts = numpy.cumsum([0] + [len(idx) for _, idx in groups]) + for p, (spinp, idxp) in enumerate(groups): + cp = numpy.asarray(mo_coeff[spinp])[:, idxp] + for q, (spinq, idxq) in enumerate(groups): + if spinp == spinq: + cq = numpy.asarray(mo_coeff[spinq])[:, idxq] + fock[starts[p]:starts[p + 1], starts[q]:starts[q + 1]] = \ + reduce(numpy.dot, (cp.conj().T, fockao[spinp], cq)) + + eri = ao2mo.kernel(self.mol, coeff, compact=False).reshape([nmo] * 4) + eri = eri.reshape(nmo * nmo, nmo * nmo) + forbidden = (orbspin[:, None] != orbspin).ravel() + eri[forbidden, :] = 0 + eri[:, forbidden] = 0 + eri = eri.reshape([nmo] * 4) + eri = eri.transpose(0, 2, 1, 3) - eri.transpose(0, 2, 3, 1) + + cpu_blocks = { + 'fock': fock, + 'oooo': eri[:nocc, :nocc, :nocc, :nocc].copy(), + 'ooov': eri[:nocc, :nocc, :nocc, nocc:].copy(), + 'oovv': eri[:nocc, :nocc, nocc:, nocc:].copy(), + 'ovov': eri[:nocc, nocc:, :nocc, nocc:].copy(), + 'ovvv': eri[:nocc, nocc:, nocc:, nocc:].copy(), + 'vvvv': eri[nocc:, nocc:, nocc:, nocc:].copy(), + } + required = sum(x.nbytes for x in cpu_blocks.values()) + available = int(self.pool.free_memory * self.mem_ratio) + if required > available: + raise MemoryError( + f'UCCSD requires {required / 1e9:.2f} GB of in-core GPU ' + f'integral storage; the configured limit is ' + f'{available / 1e9:.2f} GB') + + eris = _PhysicistsERIs() + eris.mo_coeff = mo_coeff + eris.orbspin = orbspin + eris.nocc = nocc + for name, block in cpu_blocks.items(): + setattr(eris, name, self.pool.asarray(block)) + eris.mo_energy = eris.fock.diagonal().real.copy() + return eris + + def init_amps(self, eris=None): + if eris is None: + eris = self.ao2mo() + nocc = eris.nocc + eia = (eris.mo_energy[:nocc, None] + - eris.mo_energy[None, nocc:]) + t1 = eris.fock[:nocc, nocc:].conj() / eia + t2 = eris.oovv.conj() / (eia[:, None, :, None] + + eia[None, :, None, :]) + self.emp2 = (.25 * _einsum('ijab,ijab', t2, + eris.oovv.conj())).real.item() + logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2) + return self.emp2, t1, t2 + + def _iterate(self, eris, t1, t2): + log = logger.new_logger(self) + eold = energy(self, t1, t2, eris) + adiis = pyscf_lib.diis.DIIS(self) if self.diis else None + if adiis is not None: + adiis.space = self.diis_space + converged = False + for cycle in range(self.max_cycle): + t1new, t2new = update_amps(self, t1, t2, eris) + norm1 = cupy.linalg.norm(t1new - t1) + norm2 = cupy.linalg.norm(t2new - t2) + normt = cupy.sqrt(norm1 * norm1 + norm2 * norm2).item() + if self.callback is not None: + self.callback(locals()) + if self.iterative_damping < 1: + alpha = self.iterative_damping + t1new = alpha * t1new + (1 - alpha) * t1 + t2new = alpha * t2new + (1 - alpha) * t2 + if adiis is not None and cycle >= self.diis_start_cycle: + nocc, nvir = t1new.shape + io = cupy.tril_indices(nocc, -1) + iv = cupy.tril_indices(nvir, -1) + t2packed = t2new[ + io[0][:, None], io[1][:, None], + iv[0][None, :], iv[1][None, :]].ravel() + vec = cupy.concatenate((t1new.ravel(), t2packed)).get() + vec = adiis.update(vec) + t1new, t2new = pyscf_gccsd.vector_to_amplitudes( + vec, nocc + nvir, nocc) + t1new = cupy.asarray(t1new) + t2new = cupy.asarray(t2new) + enew = energy(self, t1new, t2new, eris) + log.info('cycle = %d E_corr(UCCSD) = %.15g dE = %.9g ' + 'norm(t1,t2) = %.6g', cycle + 1, enew, + enew - eold, normt) + if abs(enew - eold) < self.conv_tol and normt < self.conv_tol_normt: + converged = True + t1, t2, eold = t1new, t2new, enew + break + t1, t2, eold = t1new, t2new, enew + return converged, eold, t1, t2 + + def kernel(self, t1=None, t2=None, eris=None, mbpt2=False): + return self.ccsd(t1, t2, eris, mbpt2=mbpt2) + + def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False): + if eris is None: + eris = self.ao2mo() + self.e_hf = self._scf.e_tot + if mbpt2 or t1 is None or t2 is None: + _, t1spin, t2spin = self.init_amps(eris) + else: + t1spin, t2spin = self._spatial_to_spin(t1, t2) + if mbpt2: + self.e_corr = self.emp2 + self.converged = True + else: + self.converged, self.e_corr, t1spin, t2spin = \ + self._iterate(eris, t1spin, t2spin) + self._t1_spin = self.pool.add('t1', t1spin) + self._t2_spin = self.pool.add('t2', t2spin) + self.t1, self.t2 = self._spin_to_spatial( + self._t1_spin, self._t2_spin) + self._finalize() + return self.e_corr, self.t1, self.t2 + + def energy(self, t1=None, t2=None, eris=None): + if eris is None: + eris = self.ao2mo() + if t1 is None or t2 is None: + t1spin, t2spin = self._t1_spin, self._t2_spin + else: + t1spin, t2spin = self._spatial_to_spin(t1, t2) + return energy(self, t1spin, t2spin, eris) + + def ccsd_t(self, t1=None, t2=None, eris=None): + if eris is None: + eris = self.ao2mo() + if t1 is None or t2 is None: + t1spin, t2spin = self._t1_spin, self._t2_spin + else: + t1spin, t2spin = self._spatial_to_spin(t1, t2) + et = ccsd_t(self, eris, t1spin, t2spin) + logger.note(self, 'UCCSD(T) correction = %.15g', et) + return et + + uccsd_t = ccsd_t