UCCSD with density fitting

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

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__

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')

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)