diff --git a/examples/cc/21-dfccsd.py b/examples/cc/21-dfccsd.py index ea3f4439f7..dc194b402b 100644 --- a/examples/cc/21-dfccsd.py +++ b/examples/cc/21-dfccsd.py @@ -1,7 +1,7 @@ #!/usr/bin/env python ''' -DF-CCSD and DF-IP/EA-EOM-CCSD for RHF is available +DF-CCSD and DF-IP/EA-EOM-CCSD for RHF/UHF is available To use DF-RCCSD object, density fitting should be enabled at SCF level. DFCCSD uses the same auxiliary basis as the DF-SCF method. It does not diff --git a/pyscf/cc/__init__.py b/pyscf/cc/__init__.py index 34d043e1cf..a8a8c15d0c 100644 --- a/pyscf/cc/__init__.py +++ b/pyscf/cc/__init__.py @@ -120,15 +120,14 @@ def RCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): def UCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): from pyscf.df.df_jk import _DFHF + from pyscf.cc import dfuccsd mf = mf.remove_soscf() if not mf.istype('UHF'): mf = mf.to_uhf() if isinstance(mf, _DFHF) and mf.with_df: - # TODO: DF-UCCSD with memory-efficient particle-particle ladder, - # similar to dfccsd.RCCSD - return uccsd.UCCSD(mf, frozen, mo_coeff, mo_occ) + return dfuccsd.UCCSD(mf, frozen, mo_coeff, mo_occ) else: return uccsd.UCCSD(mf, frozen, mo_coeff, mo_occ) UCCSD.__doc__ = uccsd.UCCSD.__doc__ diff --git a/pyscf/cc/dfccsd.py b/pyscf/cc/dfccsd.py index 17a5eb5812..9da8c0f7ea 100644 --- a/pyscf/cc/dfccsd.py +++ b/pyscf/cc/dfccsd.py @@ -51,7 +51,7 @@ def _add_vvvv(self, t1, t2, eris, out=None, with_ovvv=False, t2sym=None): to_gpu = lib.to_gpu -def _contract_vvvv_t2(mycc, mol, vvL, t2, out=None, verbose=None): +def _contract_vvvv_t2(mycc, mol, vvL, VVL, t2, out=None, verbose=None): '''Ht2 = numpy.einsum('ijcd,acdb->ijab', t2, vvvv) Args: @@ -103,7 +103,7 @@ def contract_blk_(eri, i0, i1, j0, j1): ijL = vvL0[tril2sq[i0:i1,j0:j1] - off0].reshape(-1,naux) eri = numpy.ndarray(((i1-i0)*(j1-j0),nvir_pair), buffer=eribuf) for p0, p1 in lib.prange(0, nvir_pair, vvblk): - vvL1 = _cp(vvL[p0:p1]) + vvL1 = _cp(VVL[p0:p1]) eri[:,p0:p1] = lib.ddot(ijL, vvL1.T) vvL1 = None @@ -120,7 +120,7 @@ def contract_blk_(eri, i0, i1, j0, j1): class _ChemistsERIs(ccsd._ChemistsERIs): def _contract_vvvv_t2(self, mycc, t2, direct=False, out=None, verbose=None): assert (not direct) - return _contract_vvvv_t2(mycc, self.mol, self.vvL, t2, out, verbose) + return _contract_vvvv_t2(mycc, self.mol, self.vvL, self.vvL, t2, out, verbose) def _make_df_eris(cc, mo_coeff=None): assert cc._scf.istype('RHF') diff --git a/pyscf/cc/dfuccsd.py b/pyscf/cc/dfuccsd.py new file mode 100644 index 0000000000..ad7601dcdb --- /dev/null +++ b/pyscf/cc/dfuccsd.py @@ -0,0 +1,310 @@ +#!/usr/bin/env python +# Copyright 2014-2025 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. + + +import numpy as np +from pyscf import lib +from pyscf.lib import logger +from pyscf import df +from pyscf.cc import uccsd +from pyscf.cc import ccsd +from pyscf.cc import dfccsd +from pyscf import __config__ + +MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000) + +class UCCSD(uccsd.UCCSD): + _keys = {'with_df'} + + def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): + uccsd.UCCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) + if getattr(mf, 'with_df', None): + self.with_df = mf.with_df + else: + self.with_df = df.DF(mf.mol) + self.with_df.auxbasis = df.make_auxbasis(mf.mol, mp2fit=True) + + def reset(self, mol=None): + self.with_df.reset(mol) + return uccsd.UCCSD.reset(self, mol) + + def ao2mo(self, mo_coeff=None): + return _make_df_eris(self, mo_coeff) + + def _add_vvvv(self, t1, t2, eris, out=None, with_ovvv=False, t2sym=None): + assert (not self.direct) + return uccsd.UCCSD._add_vvvv(self, t1, t2, eris, out, with_ovvv, t2sym) + + def _add_vvVV(self, t1, t2, eris, out=None): + assert (not self.direct) + return uccsd.UCCSD._add_vvVV(self, t1, t2, eris, out) + + to_gpu = lib.to_gpu + +class _ChemistsERIs(uccsd._ChemistsERIs): + def _contract_vvvv_t2(self, mycc, t2, direct=False, out=None, verbose=None): + assert (not direct) + return dfccsd._contract_vvvv_t2(mycc, self.mol, self.vvL, self.vvL, t2, out, verbose) + def _contract_VVVV_t2(self, mycc, t2, direct=False, out=None, verbose=None): + assert (not direct) + return dfccsd._contract_vvvv_t2(mycc, self.mol, self.VVL, self.VVL, t2, out, verbose) + def _contract_vvVV_t2(self, mycc, t2, direct=False, out=None, verbose=None): + assert (not direct) + return dfccsd._contract_vvvv_t2(mycc, self.mol, self.vvL, self.VVL, t2, out, verbose) + +def _cp(a): + return np.asarray(a, order='C') + +def _make_df_eris(mycc, mo_coeff=None): + assert mycc._scf.istype('UHF') + cput0 = (logger.process_clock(), logger.perf_counter()) + log = logger.Logger(mycc.stdout, mycc.verbose) + eris = _ChemistsERIs() + eris._common_init_(mycc, mo_coeff) + + moa, mob = eris.mo_coeff + nocca, noccb = eris.nocc + nao = moa.shape[0] + nmoa = moa.shape[1] + nmob = mob.shape[1] + nvira = nmoa - nocca + nvirb = nmob - noccb + nvira_pair = nvira * (nvira + 1) // 2 + nvirb_pair = nvirb * (nvirb + 1) // 2 + with_df = mycc.with_df + naux = eris.naux = with_df.get_naoaux() + + # --- Three-center integrals + # (L|aa) + Loo = np.empty((naux, nocca, nocca)) + Lov = np.empty((naux, nocca, nvira)) + Lvo = np.empty((naux, nvira, nocca)) + # (L|bb) + LOO = np.empty((naux, noccb, noccb)) + LOV = np.empty((naux, noccb, nvirb)) + LVO = np.empty((naux, nvirb, noccb)) + + # --- Four-center integrals + eris.feri = lib.H5TmpFile() + # (aa|aa) + eris.oooo = eris.feri.create_dataset('oooo', (nocca,nocca,nocca,nocca), 'f8') + eris.oovv = eris.feri.create_dataset('oovv', (nocca,nocca,nvira,nvira), 'f8', chunks=(nocca,nocca,1,nvira)) + eris.ovoo = eris.feri.create_dataset('ovoo', (nocca,nvira,nocca,nocca), 'f8', chunks=(nocca,1,nocca,nocca)) + eris.ovvo = eris.feri.create_dataset('ovvo', (nocca,nvira,nvira,nocca), 'f8', chunks=(nocca,1,nvira,nocca)) + eris.ovov = eris.feri.create_dataset('ovov', (nocca,nvira,nocca,nvira), 'f8', chunks=(nocca,1,nocca,nvira)) + # (bb|bb) + eris.OOOO = eris.feri.create_dataset('OOOO', (noccb,noccb,noccb,noccb), 'f8') + eris.OOVV = eris.feri.create_dataset('OOVV', (noccb,noccb,nvirb,nvirb), 'f8', chunks=(noccb,noccb,1,nvirb)) + eris.OVOO = eris.feri.create_dataset('OVOO', (noccb,nvirb,noccb,noccb), 'f8', chunks=(noccb,1,noccb,noccb)) + eris.OVVO = eris.feri.create_dataset('OVVO', (noccb,nvirb,nvirb,noccb), 'f8', chunks=(noccb,1,nvirb,noccb)) + eris.OVOV = eris.feri.create_dataset('OVOV', (noccb,nvirb,noccb,nvirb), 'f8', chunks=(noccb,1,noccb,nvirb)) + # (aa|bb) + eris.ooOO = eris.feri.create_dataset('ooOO', (nocca,nocca,noccb,noccb), 'f8') + eris.ooVV = eris.feri.create_dataset('ooVV', (nocca,nocca,nvirb,nvirb), 'f8', chunks=(nocca,nocca,1,nvirb)) + eris.ovOO = eris.feri.create_dataset('ovOO', (nocca,nvira,noccb,noccb), 'f8', chunks=(nocca,1,noccb,noccb)) + eris.ovVO = eris.feri.create_dataset('ovVO', (nocca,nvira,nvirb,noccb), 'f8', chunks=(nocca,1,nvirb,noccb)) + eris.ovOV = eris.feri.create_dataset('ovOV', (nocca,nvira,noccb,nvirb), 'f8', chunks=(nocca,1,noccb,nvirb)) + # (bb|aa) + eris.OOvv = eris.feri.create_dataset('OOvv', (noccb,noccb,nvira,nvira), 'f8', chunks=(noccb,noccb,1,nvira)) + eris.OVoo = eris.feri.create_dataset('OVoo', (noccb,nvirb,nocca,nocca), 'f8', chunks=(noccb,1,nocca,nocca)) + eris.OVvo = eris.feri.create_dataset('OVvo', (noccb,nvirb,nvira,nocca), 'f8', chunks=(noccb,1,nvira,nocca)) + + # nrow ~ 4e9/8/blockdim to ensure hdf5 chunk < 4GB + chunks = (min(nvira_pair, int(4e8 / with_df.blockdim)), min(naux, with_df.blockdim)) + eris.vvL = eris.feri.create_dataset('vvL', (nvira_pair,naux), 'f8', chunks=chunks) + chunks = (min(nvirb_pair, int(4e8 / with_df.blockdim)), min(naux, with_df.blockdim)) + eris.VVL = eris.feri.create_dataset('VVL', (nvirb_pair,naux), 'f8', chunks=chunks) + + # Transform three-center integrals to MO basis + p1 = 0 + for eri1 in with_df.loop(): + eri1 = lib.unpack_tril(eri1).reshape(-1, nao, nao) + # (L|aa) + Lpq = lib.einsum('Lab,ap,bq->Lpq', eri1, moa, moa) + p0, p1 = p1, p1 + Lpq.shape[0] + blk = np.s_[p0:p1] + Loo[blk] = Lpq[:, :nocca, :nocca] + Lov[blk] = Lpq[:, :nocca, nocca:] + Lvo[blk] = Lpq[:, nocca:, :nocca] + eris.vvL[:, p0:p1] = lib.pack_tril(Lpq[:, nocca:, nocca:]).T + # (L|bb) + Lpq = None + Lpq = lib.einsum('Lab,ap,bq->Lpq', eri1, mob, mob) + LOO[blk] = Lpq[:, :noccb, :noccb] + LOV[blk] = Lpq[:, :noccb, noccb:] + LVO[blk] = Lpq[:, noccb:, :noccb] + eris.VVL[:, p0:p1] = lib.pack_tril(Lpq[:, noccb:, noccb:]).T + Lpq = None + + Loo = Loo.reshape(naux, nocca * nocca) + Lov = Lov.reshape(naux, nocca * nvira) + Lvo = Lvo.reshape(naux, nocca * nvira) + LOO = LOO.reshape(naux, noccb * noccb) + LOV = LOV.reshape(naux, noccb * nvirb) + LVO = LVO.reshape(naux, noccb * nvirb) + + eris.oooo[:] = lib.ddot(Loo.T, Loo).reshape(nocca,nocca,nocca,nocca) + eris.ovoo[:] = lib.ddot(Lov.T, Loo).reshape(nocca,nvira,nocca,nocca) + eris.ovvo[:] = lib.ddot(Lov.T, Lvo).reshape(nocca,nvira,nvira,nocca) + eris.ovov[:] = lib.ddot(Lov.T, Lov).reshape(nocca,nvira,nocca,nvira) + + eris.OVoo[:] = lib.ddot(LOV.T, Loo).reshape(noccb,nvirb,nocca,nocca) + eris.OVvo[:] = lib.ddot(LOV.T, Lvo).reshape(noccb,nvirb,nvira,nocca) + + Lvo = None + + eris.OOOO[:] = lib.ddot(LOO.T, LOO).reshape(noccb,noccb,noccb,noccb) + eris.OVOO[:] = lib.ddot(LOV.T, LOO).reshape(noccb,nvirb,noccb,noccb) + eris.OVVO[:] = lib.ddot(LOV.T, LVO).reshape(noccb,nvirb,nvirb,noccb) + eris.OVOV[:] = lib.ddot(LOV.T, LOV).reshape(noccb,nvirb,noccb,nvirb) + + eris.ooOO[:] = lib.ddot(Loo.T, LOO).reshape(nocca,nocca,noccb,noccb) + eris.ovOO[:] = lib.ddot(Lov.T, LOO).reshape(nocca,nvira,noccb,noccb) + eris.ovVO[:] = lib.ddot(Lov.T, LVO).reshape(nocca,nvira,nvirb,noccb) + eris.ovOV[:] = lib.ddot(Lov.T, LOV).reshape(nocca,nvira,noccb,nvirb) + + LVO = None + + mem_now = lib.current_memory()[0] + max_memory = max(0, mycc.max_memory - mem_now) + + # eris.oovv + blksize = max(ccsd.BLKMIN, int((max_memory*.9e6/8-nocca**2*nvira_pair)/(nocca**2+naux))) + oovv_tril = np.empty((nocca * nocca, nvira_pair)) + for p0, p1 in lib.prange(0, nvira_pair, blksize): + oovv_tril[:, p0:p1] = lib.ddot(Loo.T, _cp(eris.vvL[p0:p1]).T) + eris.oovv[:] = lib.unpack_tril(oovv_tril).reshape(nocca, nocca, nvira, nvira) + oovv_tril = None + + # eris.ooVV + blksize = max(ccsd.BLKMIN, int((max_memory*.9e6/8-nocca**2*nvirb_pair)/(nocca**2+naux))) + oovv_tril = np.empty((nocca * nocca, nvirb_pair)) + for p0, p1 in lib.prange(0, nvirb_pair, blksize): + oovv_tril[:, p0:p1] = lib.ddot(Loo.T, _cp(eris.VVL[p0:p1]).T) + eris.ooVV[:] = lib.unpack_tril(oovv_tril).reshape(nocca, nocca, nvirb, nvirb) + oovv_tril = Loo = None + + mem_now = lib.current_memory()[0] + max_memory = max(0, mycc.max_memory - mem_now) + + # eris.OOvv + blksize = max(ccsd.BLKMIN, int((max_memory*.9e6/8-noccb**2*nvira_pair)/(noccb**2+naux))) + oovv_tril = np.empty((noccb * noccb, nvira_pair)) + for p0, p1 in lib.prange(0, nvira_pair, blksize): + oovv_tril[:, p0:p1] = lib.ddot(LOO.T, _cp(eris.vvL[p0:p1]).T) + eris.OOvv[:] = lib.unpack_tril(oovv_tril).reshape(noccb, noccb, nvira, nvira) + oovv_tril = None + + # eris.OOVV + blksize = max(ccsd.BLKMIN, int((max_memory*.9e6/8-noccb**2*nvirb_pair)/(noccb**2+naux))) + oovv_tril = np.empty((noccb * noccb, nvirb_pair)) + for p0, p1 in lib.prange(0, nvirb_pair, blksize): + oovv_tril[:, p0:p1] = lib.ddot(LOO.T, _cp(eris.VVL[p0:p1]).T) + eris.OOVV[:] = lib.unpack_tril(oovv_tril).reshape(noccb, noccb, nvirb, nvirb) + oovv_tril = LOO = None + + mem_now = lib.current_memory()[0] + max_memory = max(0, mycc.max_memory - mem_now) + + Lov = Lov.reshape(naux, nocca, nvira) + LOV = LOV.reshape(naux, noccb, nvirb) + + # eris.ovvv + vblk = max(nocca, int((max_memory*.15e6/8)/(nocca*nvira_pair))) + vvblk = int(min(nvira_pair, 4e8/nocca, max(4, (max_memory*.8e6/8)/(vblk*nocca+naux)))) + eris.ovvv = eris.feri.create_dataset('ovvv', (nocca, nvira, nvira_pair), 'f8', chunks=(nocca, 1, vvblk)) + for q0, q1 in lib.prange(0, nvira_pair, vvblk): + vvL = _cp(eris.vvL[q0:q1]) + for p0, p1 in lib.prange(0, nvira, vblk): + tmpLov = Lov[:, :, p0:p1].reshape(naux, -1) + eris.ovvv[:, p0:p1, q0:q1] = lib.ddot(tmpLov.T, vvL.T).reshape(nocca, p1 - p0, q1 - q0) + vvL = None + + # eris.ovVV + vblk = max(nocca, int((max_memory*.15e6/8)/(nocca*nvirb_pair))) + vvblk = int(min(nvirb_pair, 4e8/nocca, max(4, (max_memory*.8e6/8)/(vblk*nocca+naux)))) + eris.ovVV = eris.feri.create_dataset('ovVV', (nocca, nvira, nvirb_pair), 'f8', chunks=(nocca, 1, vvblk)) + for q0, q1 in lib.prange(0, nvirb_pair, vvblk): + vvL = _cp(eris.VVL[q0:q1]) + for p0, p1 in lib.prange(0, nvira, vblk): + tmpLov = Lov[:, :, p0:p1].reshape(naux, -1) + eris.ovVV[:, p0:p1, q0:q1] = lib.ddot(tmpLov.T, vvL.T).reshape(nocca, p1 - p0, q1 - q0) + vvL = None + Lov = None + + mem_now = lib.current_memory()[0] + max_memory = max(0, mycc.max_memory - mem_now) + + # eris.OVvv + vblk = max(noccb, int((max_memory*.15e6/8)/(noccb*nvira_pair))) + vvblk = int(min(nvira_pair, 4e8/noccb, max(4, (max_memory*.8e6/8)/(vblk*noccb+naux)))) + eris.OVvv = eris.feri.create_dataset('OVvv', (noccb, nvirb, nvira_pair), 'f8', chunks=(noccb, 1, vvblk)) + for q0, q1 in lib.prange(0, nvira_pair, vvblk): + vvL = _cp(eris.vvL[q0:q1]) + for p0, p1 in lib.prange(0, nvirb, vblk): + tmpLov = LOV[:, :, p0:p1].reshape(naux, -1) + eris.OVvv[:, p0:p1, q0:q1] = lib.ddot(tmpLov.T, vvL.T).reshape(noccb, p1 - p0, q1 - q0) + vvL = None + + # eris.OVVV + vblk = max(noccb, int((max_memory*.15e6/8)/(noccb*nvirb_pair))) + vvblk = int(min(nvirb_pair, 4e8/noccb, max(4, (max_memory*.8e6/8)/(vblk*noccb+naux)))) + eris.OVVV = eris.feri.create_dataset('OVVV', (noccb, nvirb, nvirb_pair), 'f8', chunks=(noccb, 1, vvblk)) + for q0, q1 in lib.prange(0, nvirb_pair, vvblk): + vvL = _cp(eris.VVL[q0:q1]) + for p0, p1 in lib.prange(0, nvirb, vblk): + tmpLov = LOV[:, :, p0:p1].reshape(naux, -1) + eris.OVVV[:, p0:p1, q0:q1] = lib.ddot(tmpLov.T, vvL.T).reshape(noccb, p1 - p0, q1 - q0) + vvL = None + LOV = None + + log.timer('DF-UCCSD integral transformation', *cput0) + return eris + +if __name__ == '__main__': + from pyscf import gto + from pyscf import scf + + mol = gto.Mole() + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + + mol.basis = 'cc-pvdz' + mol.build() + mf = scf.UHF(mol).density_fit('weigend').run() + mycc = UCCSD(mf).run() + print(mycc.e_corr - -0.2133709727796199) + + print("IP energies... (right eigenvector)") + e,v = mycc.ipccsd(nroots=8) + print(e) + print(e[0] - 0.4336428577342009) + print(e[2] - 0.5188000951518845) + print(e[4] - 0.6785158684375829) + + print("EA energies... (right eigenvector)") + e,v = mycc.eaccsd(nroots=8) + print(e) + print(e[0] - 0.1673013569134136) + print(e[2] - 0.2399984284491973) + print(e[4] - 0.5096018470162480) + + e, v = mycc.eeccsd(nroots=4) + print(e[0] - 0.2757563806054133) + print(e[1] - 0.2757563806171079) + print(e[2] - 0.2757563806183815) + print(e[3] - 0.3006896721085447) diff --git a/pyscf/cc/test/test_dfuccsd.py b/pyscf/cc/test/test_dfuccsd.py new file mode 100644 index 0000000000..296b131ddd --- /dev/null +++ b/pyscf/cc/test/test_dfuccsd.py @@ -0,0 +1,254 @@ +#!/usr/bin/env python +# Copyright 2014-2025 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. + +import unittest +import numpy + +from pyscf import gto, lib +from pyscf import scf +from pyscf import cc +from pyscf.cc import dfuccsd, eom_uccsd + +def setUpModule(): + global mol, rhf, mf, ucc, ucc1, nocca, nvira, noccb, nvirb + mol = gto.Mole() + mol.verbose = 7 + mol.output = '/dev/null' + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + + mol.basis = '631g' + mol.build() + rhf = scf.RHF(mol).density_fit(auxbasis='weigend') + rhf.conv_tol_grad = 1e-8 + rhf.kernel() + mf = scf.addons.convert_to_uhf(rhf) + + ucc = dfuccsd.UCCSD(mf).run(conv_tol=1e-10) + + mol1 = mol.copy() + mol1.spin = 2 + mol1.build() + mf0 = scf.UHF(mol1).density_fit(auxbasis='weigend') + mf0 = mf0.run(conv_tol=1e-12) + mf1 = mf0.copy() + + nocca, noccb = mol1.nelec + nmo = mol1.nao_nr() + nvira, nvirb = nmo - nocca, nmo-noccb + numpy.random.seed(12) + mf1.mo_coeff = numpy.random.random((2, nmo, nmo)) - 0.5 + gmf = scf.addons.convert_to_ghf(mf1) + orbspin = gmf.mo_coeff.orbspin + + ucc1 = dfuccsd.UCCSD(mf1) + + numpy.random.seed(11) + no = nocca + noccb + nv = nvira + nvirb + r1 = numpy.random.random((no, nv)) - .9 + r2 = numpy.random.random((no, no, nv, nv)) - .9 + r2 = r2 - r2.transpose(1, 0, 2, 3) + r2 = r2 - r2.transpose(0, 1, 3, 2) + r1 = cc.addons.spin2spatial(r1, orbspin) + r2 = cc.addons.spin2spatial(r2, orbspin) + r1, r2 = eom_uccsd.vector_to_amplitudes_ee( + eom_uccsd.amplitudes_to_vector_ee(r1, r2), ucc1.nmo, ucc1.nocc) + ucc1.t1 = r1 + ucc1.t2 = r2 + +def tearDownModule(): + global mol, rhf, mf, ucc, ucc1 + mol.stdout.close() + del mol, rhf, mf, ucc, ucc1 + +class KnownValues(unittest.TestCase): + def test_with_df(self): + self.assertAlmostEqual(ucc.e_tot, -76.118403942938741, 6) + numpy.random.seed(1) + mo_coeff = numpy.random.random(mf.mo_coeff.shape) + eris = cc.uccsd.UCCSD(mf).ao2mo(mo_coeff) + self.assertAlmostEqual(lib.fp(eris.oooo), 4.96203346086189, 11) + self.assertAlmostEqual(lib.fp(eris.ovoo), -1.36660785172450, 11) + self.assertAlmostEqual(lib.fp(eris.ovov), 125.80972789115610, 11) + self.assertAlmostEqual(lib.fp(eris.oovv), 55.12252557132108, 11) + self.assertAlmostEqual(lib.fp(eris.ovvo), 133.48517302161093, 11) + self.assertAlmostEqual(lib.fp(eris.ovvv), 59.41874702857587, 11) + self.assertAlmostEqual(lib.fp(eris.vvvv), 43.56245722797580, 11) + self.assertAlmostEqual(lib.fp(eris.OOOO),-407.06688974686460, 11) + self.assertAlmostEqual(lib.fp(eris.OVOO), 56.27143890752203, 11) + self.assertAlmostEqual(lib.fp(eris.OVOV),-287.70639282707754, 11) + self.assertAlmostEqual(lib.fp(eris.OOVV), -85.47974577569636, 11) + self.assertAlmostEqual(lib.fp(eris.OVVO),-228.18507149174869, 11) + self.assertAlmostEqual(lib.fp(eris.OVVV), -10.71277146732805, 11) + self.assertAlmostEqual(lib.fp(eris.VVVV), -89.90585731601215, 11) + self.assertAlmostEqual(lib.fp(eris.ooOO),-336.66684771631481, 11) + self.assertAlmostEqual(lib.fp(eris.ovOO), -16.41425875389384, 11) + self.assertAlmostEqual(lib.fp(eris.ovOV), 231.59582076434012, 11) + self.assertAlmostEqual(lib.fp(eris.ooVV), 20.33912803565480, 11) + self.assertAlmostEqual(lib.fp(eris.ovVO), 206.47963466976572, 11) + self.assertAlmostEqual(lib.fp(eris.ovVV), -71.27033893003525, 11) + self.assertAlmostEqual(lib.fp(eris.vvVV), 172.46980782354174, 11) + self.assertAlmostEqual(lib.fp(eris.OVoo), -19.93667104163490, 11) + self.assertAlmostEqual(lib.fp(eris.OOvv), -27.76467892422170, 11) + self.assertAlmostEqual(lib.fp(eris.OVvo),-140.08585937729470, 11) + self.assertAlmostEqual(lib.fp(eris.OVvv), 40.69754434429569, 11) + + def test_df_ipccsd(self): + eom = ucc.eomip_method() + e,v = eom.kernel(nroots=1, koopmans=False) + self.assertAlmostEqual(e, 0.42788191076505, 5) + e,v = ucc.ipccsd(nroots=8) + self.assertAlmostEqual(e[0], 0.42788191078779, 5) + self.assertAlmostEqual(e[2], 0.50229583182010, 5) + self.assertAlmostEqual(e[4], 0.68557653039451, 5) + + e,v = ucc.ipccsd(nroots=4, guess=v[:4]) + self.assertAlmostEqual(e[2], 0.50229583182010, 5) + + def test_df_ipccsd_koopmans(self): + e,v = ucc.ipccsd(nroots=8, koopmans=True) + self.assertAlmostEqual(e[0], 0.42788191078779, 5) + self.assertAlmostEqual(e[2], 0.50229583182010, 5) + self.assertAlmostEqual(e[4], 0.68557653039451, 5) + + def test_df_eaccsd(self): + eom = ucc.eomea_method() + e,v = eom.kernel(nroots=1, koopmans=False) + self.assertAlmostEqual(e, 0.19038860308234, 5) + e,v = ucc.eaccsd(nroots=8) + self.assertAlmostEqual(e[0], 0.19038860335165, 5) + self.assertAlmostEqual(e[2], 0.28339727179550, 5) + self.assertAlmostEqual(e[4], 0.52224978073482, 5) + + e,v = ucc.eaccsd(nroots=4, guess=v[:4]) + self.assertAlmostEqual(e[2], 0.28339727179551, 5) + + def test_df_eaccsd_koopmans(self): + e,v = ucc.eaccsd(nroots=6, koopmans=True) + self.assertAlmostEqual(e[0], 0.19038860335165, 5) + self.assertAlmostEqual(e[2], 0.28339727179550, 5) + + def test_eomee(self): + self.assertAlmostEqual(ucc.e_corr, -0.13519304930252, 5) + eom = ucc.eomee_method() + e,v = eom.kernel(nroots=1, koopmans=False) + self.assertAlmostEqual(e, 0.28107576231548, 5) + + e,v = ucc.eeccsd(nroots=4) + self.assertAlmostEqual(e[0], 0.28107576231548, 5) + self.assertAlmostEqual(e[1], 0.28107576231548, 5) + self.assertAlmostEqual(e[2], 0.28107576231548, 5) + self.assertAlmostEqual(e[3], 0.30810935830310, 5) + + e,v = ucc.eeccsd(nroots=4, guess=v[:4]) + self.assertAlmostEqual(e[3], 0.30810935830310, 5) + + def test_eomee_ccsd_spin_keep(self): + e, v = ucc.eomee_ccsd(nroots=2, koopmans=False) + self.assertAlmostEqual(e[0], 0.28107576231548, 5) + self.assertAlmostEqual(e[1], 0.30810935830310, 5) + + e, v = ucc.eomee_ccsd(nroots=2, koopmans=True) + self.assertAlmostEqual(e[0], 0.28107576231548, 5) + self.assertAlmostEqual(e[1], 0.30810935830310, 5) + + def test_df_eomee_ccsd_matvec(self): + numpy.random.seed(10) + r1 = [numpy.random.random((nocca,nvira))-.9, + numpy.random.random((noccb,nvirb))-.9] + r2 = [numpy.random.random((nocca,nocca,nvira,nvira))-.9, + numpy.random.random((nocca,noccb,nvira,nvirb))-.9, + numpy.random.random((noccb,noccb,nvirb,nvirb))-.9] + r2[0] = r2[0] - r2[0].transpose(1,0,2,3) + r2[0] = r2[0] - r2[0].transpose(0,1,3,2) + r2[2] = r2[2] - r2[2].transpose(1,0,2,3) + r2[2] = r2[2] - r2[2].transpose(0,1,3,2) + + uee1 = eom_uccsd.EOMEESpinKeep(ucc1) + vec = uee1.amplitudes_to_vector(r1,r2) + vec1 = uee1.matvec(vec) + + r1, r2 = uee1.vector_to_amplitudes(vec1) + self.assertAlmostEqual(lib.fp(vec1), 49.6181987143722836, 9) + + def test_df_eomee_ccsd_diag(self): + vec1, vec2 = eom_uccsd.EOMEE(ucc1).get_diag() + self.assertAlmostEqual(lib.fp(vec1), 67.5604889872366812, 9) + self.assertAlmostEqual(lib.fp(vec2), 161.2361154461494834, 9) + + def test_df_eomee_init_guess(self): + uee = eom_uccsd.EOMEESpinKeep(ucc1) + diag = uee.get_diag()[0] + guess = uee.get_init_guess(nroots=1, koopmans=False, diag=diag) + self.assertAlmostEqual(lib.fp(guess[0]), -0.4558886175539251, 9) + + guess = uee.get_init_guess(nroots=1, koopmans=True, diag=diag) + self.assertAlmostEqual(lib.fp(guess[0]), -0.8438701329927379, 9) + + guess = uee.get_init_guess(nroots=4, koopmans=False, diag=diag) + self.assertAlmostEqual(lib.fp(guess), 0.3749885550995292, 9) + + guess = uee.get_init_guess(nroots=4, koopmans=True, diag=diag) + self.assertAlmostEqual(lib.fp(guess), -0.3812403236695565, 9) + + def test_df_eomsf_ccsd_matvec(self): + numpy.random.seed(10) + myeom = eom_uccsd.EOMEESpinFlip(ucc1) + vec = numpy.random.random(myeom.vector_size()) - .9 + vec1 = myeom.matvec(vec) + self.assertAlmostEqual(lib.fp(vec1), -1655.6005996668061471, 8) + + def test_ao2mo(self): + numpy.random.seed(2) + mo = numpy.random.random(mf.mo_coeff.shape) + mycc = cc.CCSD(mf).density_fit(auxbasis='ccpvdz-ri') + mycc.max_memory = 0 + eri_df = mycc.ao2mo(mo) + + self.assertAlmostEqual(lib.fp(eri_df.oooo),-493.98003157749906, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovoo),-203.89515661847452, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovov), -57.62195194777571, 9) + self.assertAlmostEqual(lib.fp(eri_df.oovv), -91.84858398271636, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovvo), -14.88387735916913, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovvv), -24.35941895353339, 9) + self.assertTrue(eri_df.vvvv is None) + self.assertAlmostEqual(lib.fp(eri_df.OOOO), 144.14457267205376, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVOO),-182.57213907795114, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVOV), 462.82041314370520, 9) + self.assertAlmostEqual(lib.fp(eri_df.OOVV), 165.48054495591805, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVVO), 499.35368006930844, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVVV), 117.15437910286980, 9) + self.assertTrue(eri_df.VVVV is None) + self.assertAlmostEqual(lib.fp(eri_df.ooOO), -13.89000382034967, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovOO),-256.60884544897004, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovOV), -93.37973542764470, 9) + self.assertAlmostEqual(lib.fp(eri_df.ooVV), -35.14359736260144, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovVO), -80.53112424767198, 9) + self.assertAlmostEqual(lib.fp(eri_df.ovVV), -1.10606382100687, 9) + self.assertTrue(eri_df.vvVV is None) + self.assertAlmostEqual(lib.fp(eri_df.OVoo),-364.17910423297951, 9) + self.assertAlmostEqual(lib.fp(eri_df.OOvv), -60.94153644542936, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVvo), 355.78017614135643, 9) + self.assertAlmostEqual(lib.fp(eri_df.OVvv), 57.10096926407320, 9) + self.assertAlmostEqual(lib.fp(eri_df.vvL), -0.51651775168057, 9) + self.assertAlmostEqual(lib.fp(eri_df.VVL), -5.60414552806429, 9) + + +if __name__ == "__main__": + print("Full Tests for DFUCCSD") + unittest.main() diff --git a/pyscf/cc/uccsd.py b/pyscf/cc/uccsd.py index fd078a3c96..bf90ee48a5 100644 --- a/pyscf/cc/uccsd.py +++ b/pyscf/cc/uccsd.py @@ -87,7 +87,7 @@ def update_amps(cc, t1, t2, eris): for p0,p1 in lib.prange(0, nocca, blksize): ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1] ovvv = ovvv - ovvv.transpose(0,3,2,1) - Fvva += np.einsum('mf,mfae->ae', t1a[p0:p1], ovvv) + Fvva += np.einsum('mf,mfae->ae', t1a[p0:p1], ovvv, optimize=True) wovvo[p0:p1] += lib.einsum('jf,mebf->mbej', t1a, ovvv) u1a += 0.5*lib.einsum('mief,meaf->ia', t2aa[p0:p1], ovvv) u2aa[:,p0:p1] += lib.einsum('ie,mbea->imab', t1a, ovvv.conj()) @@ -100,7 +100,7 @@ def update_amps(cc, t1, t2, eris): for p0,p1 in lib.prange(0, noccb, blksize): OVVV = eris.get_OVVV(slice(p0,p1)) # OVVV = eris.OVVV[p0:p1] OVVV = OVVV - OVVV.transpose(0,3,2,1) - Fvvb += np.einsum('mf,mfae->ae', t1b[p0:p1], OVVV) + Fvvb += np.einsum('mf,mfae->ae', t1b[p0:p1], OVVV, optimize=True) wOVVO[p0:p1] = lib.einsum('jf,mebf->mbej', t1b, OVVV) u1b += 0.5*lib.einsum('MIEF,MEAF->IA', t2bb[p0:p1], OVVV) u2bb[:,p0:p1] += lib.einsum('ie,mbea->imab', t1b, OVVV.conj()) @@ -112,7 +112,7 @@ def update_amps(cc, t1, t2, eris): blksize = max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvira*nvirb**2*3+1))) for p0,p1 in lib.prange(0, nocca, blksize): ovVV = eris.get_ovVV(slice(p0,p1)) # ovVV = eris.ovVV[p0:p1] - Fvvb += np.einsum('mf,mfAE->AE', t1a[p0:p1], ovVV) + Fvvb += np.einsum('mf,mfAE->AE', t1a[p0:p1], ovVV, optimize=True) woVvO[p0:p1] = lib.einsum('JF,meBF->mBeJ', t1b, ovVV) woVVo[p0:p1] = lib.einsum('jf,mfBE->mBEj',-t1a, ovVV) u1b += lib.einsum('mIeF,meAF->IA', t2ab[p0:p1], ovVV) @@ -125,7 +125,7 @@ def update_amps(cc, t1, t2, eris): blksize = max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvirb*nvira**2*3+1))) for p0,p1 in lib.prange(0, noccb, blksize): OVvv = eris.get_OVvv(slice(p0,p1)) # OVvv = eris.OVvv[p0:p1] - Fvva += np.einsum('MF,MFae->ae', t1b[p0:p1], OVvv) + Fvva += np.einsum('MF,MFae->ae', t1b[p0:p1], OVvv, optimize=True) wOvVo[p0:p1] = lib.einsum('jf,MEbf->MbEj', t1a, OVvv) wOvvO[p0:p1] = lib.einsum('JF,MFbe->MbeJ',-t1b, OVvv) u1a += lib.einsum('iMfE,MEaf->ia', t2ab[:,p0:p1], OVvv) @@ -143,7 +143,7 @@ def update_amps(cc, t1, t2, eris): u2aa += lib.einsum('mnab,mnij->ijab', tauaa, Woooo*.5) Woooo = tauaa = None ovoo = eris_ovoo - eris_ovoo.transpose(2,1,0,3) - Fooa += np.einsum('ne,nemi->mi', t1a, ovoo) + Fooa += np.einsum('ne,nemi->mi', t1a, ovoo, optimize=True) u1a += 0.5*lib.einsum('mnae,meni->ia', t2aa, ovoo) wovvo += lib.einsum('nb,nemj->mbej', t1a, ovoo) ovoo = eris_ovoo = None @@ -152,7 +152,7 @@ def update_amps(cc, t1, t2, eris): ovov = eris_ovov - eris_ovov.transpose(0,3,2,1) Fvva -= .5 * lib.einsum('mnaf,menf->ae', tilaa, ovov) Fooa += .5 * lib.einsum('inef,menf->mi', tilaa, ovov) - Fova = np.einsum('nf,menf->me',t1a, ovov) + Fova = np.einsum('nf,menf->me',t1a, ovov, optimize=True) u2aa += ovov.conj().transpose(0,2,1,3) * .5 wovvo -= 0.5*lib.einsum('jnfb,menf->mbej', t2aa, ovov) woVvO += 0.5*lib.einsum('nJfB,menf->mBeJ', t2ab, ovov) @@ -169,7 +169,7 @@ def update_amps(cc, t1, t2, eris): u2bb += lib.einsum('mnab,mnij->ijab', taubb, WOOOO*.5) WOOOO = taubb = None OVOO = eris_OVOO - eris_OVOO.transpose(2,1,0,3) - Foob += np.einsum('ne,nemi->mi', t1b, OVOO) + Foob += np.einsum('ne,nemi->mi', t1b, OVOO, optimize=True) u1b += 0.5*lib.einsum('mnae,meni->ia', t2bb, OVOO) wOVVO += lib.einsum('nb,nemj->mbej', t1b, OVOO) OVOO = eris_OVOO = None @@ -178,7 +178,7 @@ def update_amps(cc, t1, t2, eris): OVOV = eris_OVOV - eris_OVOV.transpose(0,3,2,1) Fvvb -= .5 * lib.einsum('MNAF,MENF->AE', tilbb, OVOV) Foob += .5 * lib.einsum('inef,menf->mi', tilbb, OVOV) - Fovb = np.einsum('nf,menf->me',t1b, OVOV) + Fovb = np.einsum('nf,menf->me',t1b, OVOV, optimize=True) u2bb += OVOV.conj().transpose(0,2,1,3) * .5 wOVVO -= 0.5*lib.einsum('jnfb,menf->mbej', t2bb, OVOV) wOvVo += 0.5*lib.einsum('jNbF,MENF->MbEj', t2ab, OVOV) @@ -188,11 +188,11 @@ def update_amps(cc, t1, t2, eris): eris_OVoo = np.asarray(eris.OVoo) eris_ovOO = np.asarray(eris.ovOO) - Fooa += np.einsum('NE,NEmi->mi', t1b, eris_OVoo) + Fooa += np.einsum('NE,NEmi->mi', t1b, eris_OVoo, optimize=True) u1a -= lib.einsum('nMaE,MEni->ia', t2ab, eris_OVoo) wOvVo -= lib.einsum('nb,MEnj->MbEj', t1a, eris_OVoo) woVVo += lib.einsum('NB,NEmj->mBEj', t1b, eris_OVoo) - Foob += np.einsum('ne,neMI->MI', t1a, eris_ovOO) + Foob += np.einsum('ne,neMI->MI', t1a, eris_ovOO, optimize=True) u1b -= lib.einsum('mNeA,meNI->IA', t2ab, eris_ovOO) woVvO -= lib.einsum('NB,meNJ->mBeJ', t1b, eris_ovOO) wOvvO += lib.einsum('nb,neMJ->MbeJ', t1a, eris_ovOO) @@ -211,8 +211,8 @@ def update_amps(cc, t1, t2, eris): Fvvb -= lib.einsum('nMfA,nfME->AE', tilab, eris_ovOV) Fooa += lib.einsum('iNeF,meNF->mi', tilab, eris_ovOV) Foob += lib.einsum('nIfE,nfME->MI', tilab, eris_ovOV) - Fova += np.einsum('NF,meNF->me',t1b, eris_ovOV) - Fovb += np.einsum('nf,nfME->ME',t1a, eris_ovOV) + Fova += np.einsum('NF,meNF->me',t1b, eris_ovOV, optimize=True) + Fovb += np.einsum('nf,nfME->ME',t1a, eris_ovOV, optimize=True) u2ab += eris_ovOV.conj().transpose(0,2,1,3) wovvo += 0.5*lib.einsum('jNbF,meNF->mbej', t2ab, eris_ovOV) wOVVO += 0.5*lib.einsum('nJfB,nfME->MBEJ', t2ab, eris_ovOV) @@ -231,23 +231,23 @@ def update_amps(cc, t1, t2, eris): Fova += fova Fovb += fovb u1a += fova.conj() - u1a += np.einsum('ie,ae->ia', t1a, Fvva) - u1a -= np.einsum('ma,mi->ia', t1a, Fooa) - u1a -= np.einsum('imea,me->ia', t2aa, Fova) - u1a += np.einsum('iMaE,ME->ia', t2ab, Fovb) + u1a += np.einsum('ie,ae->ia', t1a, Fvva, optimize=True) + u1a -= np.einsum('ma,mi->ia', t1a, Fooa, optimize=True) + u1a -= np.einsum('imea,me->ia', t2aa, Fova, optimize=True) + u1a += np.einsum('iMaE,ME->ia', t2ab, Fovb, optimize=True) u1b += fovb.conj() - u1b += np.einsum('ie,ae->ia',t1b,Fvvb) - u1b -= np.einsum('ma,mi->ia',t1b,Foob) - u1b -= np.einsum('imea,me->ia', t2bb, Fovb) - u1b += np.einsum('mIeA,me->IA', t2ab, Fova) + u1b += np.einsum('ie,ae->ia', t1b, Fvvb, optimize=True) + u1b -= np.einsum('ma,mi->ia', t1b, Foob, optimize=True) + u1b -= np.einsum('imea,me->ia', t2bb, Fovb, optimize=True) + u1b += np.einsum('mIeA,me->IA', t2ab, Fova, optimize=True) eris_oovv = np.asarray(eris.oovv) eris_ovvo = np.asarray(eris.ovvo) wovvo -= eris_oovv.transpose(0,2,3,1) wovvo += eris_ovvo.transpose(0,2,1,3) oovv = eris_oovv - eris_ovvo.transpose(0,3,2,1) - u1a-= np.einsum('nf,niaf->ia', t1a, oovv) - tmp1aa = lib.einsum('ie,mjbe->mbij', t1a, oovv) + u1a-= np.einsum('nf,niaf->ia', t1a, oovv, optimize=True) + tmp1aa = lib.einsum('ie,mjbe->mbij', t1a, oovv) u2aa += 2*lib.einsum('ma,mbij->ijab', t1a, tmp1aa) eris_ovvo = eris_oovv = oovv = tmp1aa = None @@ -256,8 +256,8 @@ def update_amps(cc, t1, t2, eris): wOVVO -= eris_OOVV.transpose(0,2,3,1) wOVVO += eris_OVVO.transpose(0,2,1,3) OOVV = eris_OOVV - eris_OVVO.transpose(0,3,2,1) - u1b-= np.einsum('nf,niaf->ia', t1b, OOVV) - tmp1bb = lib.einsum('ie,mjbe->mbij', t1b, OOVV) + u1b-= np.einsum('nf,niaf->ia', t1b, OOVV, optimize=True) + tmp1bb = lib.einsum('ie,mjbe->mbij', t1b, OOVV) u2bb += 2*lib.einsum('ma,mbij->ijab', t1b, tmp1bb) eris_OVVO = eris_OOVV = OOVV = None @@ -265,7 +265,7 @@ def update_amps(cc, t1, t2, eris): eris_ovVO = np.asarray(eris.ovVO) woVVo -= eris_ooVV.transpose(0,2,3,1) woVvO += eris_ovVO.transpose(0,2,1,3) - u1b+= np.einsum('nf,nfAI->IA', t1a, eris_ovVO) + u1b+= np.einsum('nf,nfAI->IA', t1a, eris_ovVO, optimize=True) tmp1ab = lib.einsum('ie,meBJ->mBiJ', t1a, eris_ovVO) tmp1ab+= lib.einsum('IE,mjBE->mBjI', t1b, eris_ooVV) u2ab -= lib.einsum('ma,mBiJ->iJaB', t1a, tmp1ab) @@ -275,7 +275,7 @@ def update_amps(cc, t1, t2, eris): eris_OVvo = np.asarray(eris.OVvo) wOvvO -= eris_OOvv.transpose(0,2,3,1) wOvVo += eris_OVvo.transpose(0,2,1,3) - u1a+= np.einsum('NF,NFai->ia', t1b, eris_OVvo) + u1a+= np.einsum('NF,NFai->ia', t1b, eris_OVvo, optimize=True) tmp1ba = lib.einsum('IE,MEbj->MbIj', t1b, eris_OVvo) tmp1ba+= lib.einsum('ie,MJbe->MbJi', t1a, eris_OOvv) u2ab -= lib.einsum('MA,MbIj->jIbA', t1b, tmp1ba) @@ -354,18 +354,18 @@ def energy(cc, t1=None, t2=None, eris=None): eris_ovOV = np.asarray(eris.ovOV) fova = eris.focka[:nocca,nocca:] fovb = eris.fockb[:noccb,noccb:] - e = np.einsum('ia,ia', fova, t1a) - e += np.einsum('ia,ia', fovb, t1b) - e += 0.25*np.einsum('ijab,iajb',t2aa,eris_ovov) - e -= 0.25*np.einsum('ijab,ibja',t2aa,eris_ovov) - e += 0.25*np.einsum('ijab,iajb',t2bb,eris_OVOV) - e -= 0.25*np.einsum('ijab,ibja',t2bb,eris_OVOV) - e += np.einsum('iJaB,iaJB',t2ab,eris_ovOV) - e += 0.5*np.einsum('ia,jb,iajb',t1a,t1a,eris_ovov) - e -= 0.5*np.einsum('ia,jb,ibja',t1a,t1a,eris_ovov) - e += 0.5*np.einsum('ia,jb,iajb',t1b,t1b,eris_OVOV) - e -= 0.5*np.einsum('ia,jb,ibja',t1b,t1b,eris_OVOV) - e += np.einsum('ia,jb,iajb',t1a,t1b,eris_ovOV) + e = np.einsum('ia,ia', fova, t1a, optimize=True) + e += np.einsum('ia,ia', fovb, t1b, optimize=True) + e += 0.25*np.einsum('ijab,iajb',t2aa, eris_ovov, optimize=True) + e -= 0.25*np.einsum('ijab,ibja',t2aa, eris_ovov, optimize=True) + e += 0.25*np.einsum('ijab,iajb',t2bb, eris_OVOV, optimize=True) + e -= 0.25*np.einsum('ijab,ibja',t2bb, eris_OVOV, optimize=True) + e += np.einsum('iJaB,iaJB',t2ab, eris_ovOV, optimize=True) + e += 0.5*np.einsum('ia,jb,iajb',t1a, t1a, eris_ovov, optimize=True) + e -= 0.5*np.einsum('ia,jb,ibja',t1a, t1a, eris_ovov, optimize=True) + e += 0.5*np.einsum('ia,jb,iajb',t1b, t1b, eris_OVOV, optimize=True) + e -= 0.5*np.einsum('ia,jb,ibja',t1b, t1b, eris_OVOV, optimize=True) + e += np.einsum('ia,jb,iajb',t1a, t1b, eris_ovOV, optimize=True) if abs(e.imag) > 1e-4: logger.warn(cc, 'Non-zero imaginary part found in UCCSD energy %s', e) return e.real @@ -579,11 +579,11 @@ def init_amps(self, eris=None): t2bb = eris_OVOV.transpose(0,2,1,3).conj() / lib.direct_sum('ia+jb->ijab', eia_b, eia_b) t2aa = t2aa - t2aa.transpose(0,1,3,2) t2bb = t2bb - t2bb.transpose(0,1,3,2) - e = np.einsum('iJaB,iaJB', t2ab, eris_ovOV) - e += 0.25*np.einsum('ijab,iajb', t2aa, eris_ovov) - e -= 0.25*np.einsum('ijab,ibja', t2aa, eris_ovov) - e += 0.25*np.einsum('ijab,iajb', t2bb, eris_OVOV) - e -= 0.25*np.einsum('ijab,ibja', t2bb, eris_OVOV) + e = np.einsum('iJaB,iaJB', t2ab, eris_ovOV, optimize=True) + e += 0.25*np.einsum('ijab,iajb', t2aa, eris_ovov, optimize=True) + e -= 0.25*np.einsum('ijab,ibja', t2aa, eris_ovov, optimize=True) + e += 0.25*np.einsum('ijab,iajb', t2bb, eris_OVOV, optimize=True) + e -= 0.25*np.einsum('ijab,ibja', t2bb, eris_OVOV, optimize=True) self.emp2 = e.real logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2) logger.timer(self, 'init mp2', *time0) @@ -687,11 +687,10 @@ def ao2mo(self, mo_coeff=None): return _make_eris_incore(self, mo_coeff) elif getattr(self._scf, 'with_df', None): - # TODO: Uncomment once there is an unrestricted DF-CCSD implementation - #logger.warn(self, 'UCCSD detected DF being used in the HF object. ' - # 'MO integrals are computed based on the DF 3-index tensors.\n' - # 'It\'s recommended to use dfccsd.CCSD for the ' - # 'DF-CCSD calculations') + logger.warn(self, 'UCCSD detected DF being used in the HF object. ' + 'MO integrals are computed based on the DF 3-index tensors.\n' + 'It\'s recommended to use dfuccsd.UCCSD for the ' + 'DF-UCCSD calculations') return _make_df_eris_outcore(self, mo_coeff) else: return _make_eris_outcore(self, mo_coeff) @@ -732,6 +731,16 @@ def eomee_method(self): from pyscf.cc import eom_uccsd return eom_uccsd.EOMEE(self) + def density_fit(self, auxbasis=None, with_df=None): + from pyscf.cc import dfuccsd + mycc = dfuccsd.UCCSD(self._scf, self.frozen, self.mo_coeff, self.mo_occ) + if with_df is not None: + mycc.with_df = with_df + if mycc.with_df.auxbasis != auxbasis: + mycc.with_df = mycc.with_df.copy() + mycc.with_df.auxbasis = auxbasis + return mycc + def nuc_grad_method(self): from pyscf.grad import uccsd return uccsd.Gradients(self) @@ -1183,8 +1192,8 @@ def make_tau(t2, t1, r1, fac=1, out=None): return tau1aa, tau1ab, tau1bb def make_tau_aa(t2aa, t1a, r1a, fac=1, out=None): - tau1aa = np.einsum('ia,jb->ijab', t1a, r1a) - tau1aa-= np.einsum('ia,jb->jiab', t1a, r1a) + tau1aa = np.einsum('ia,jb->ijab', t1a, r1a, optimize=True) + tau1aa-= np.einsum('ia,jb->jiab', t1a, r1a, optimize=True) tau1aa = tau1aa - tau1aa.transpose(0,1,3,2) tau1aa *= fac * .5 tau1aa += t2aa @@ -1193,8 +1202,8 @@ def make_tau_aa(t2aa, t1a, r1a, fac=1, out=None): def make_tau_ab(t2ab, t1, r1, fac=1, out=None): t1a, t1b = t1 r1a, r1b = r1 - tau1ab = np.einsum('ia,jb->ijab', t1a, r1b) - tau1ab+= np.einsum('ia,jb->ijab', r1a, t1b) + tau1ab = np.einsum('ia,jb->ijab', t1a, r1b, optimize=True) + tau1ab+= np.einsum('ia,jb->ijab', r1a, t1b, optimize=True) tau1ab *= fac * .5 tau1ab += t2ab return tau1ab