Source code for pyscf.cc.gccsd_t_rdm

#!/usr/bin/env python
# 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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

import numpy
from pyscf import lib
from pyscf.cc import gccsd_rdm

def _gamma1_intermediates(mycc, t1, t2, l1, l2, eris=None):
    doo, dov, dvo, dvv = gccsd_rdm._gamma1_intermediates(mycc, t1, t2, l1, l2)

    if eris is None: eris = mycc.ao2mo()

    nocc, nvir = t1.shape
    bcei = numpy.asarray(eris.ovvv).conj().transpose(3,2,1,0)
    majk = numpy.asarray(eris.ooov).conj().transpose(2,3,0,1)
    bcjk = numpy.asarray(eris.oovv).conj().transpose(2,3,0,1)

    mo_e = eris.mo_energy
    eia = mo_e[:nocc,None] - mo_e[nocc:]
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eia, eia, eia)

    t3c =(numpy.einsum('jkae,bcei->ijkabc', t2, bcei) -
          numpy.einsum('imbc,majk->ijkabc', t2, majk))
    t3c = t3c - t3c.transpose(0,1,2,4,3,5) - t3c.transpose(0,1,2,5,4,3)
    t3c = t3c - t3c.transpose(1,0,2,3,4,5) - t3c.transpose(2,1,0,3,4,5)
    t3c /= d3

    t3d = numpy.einsum('ia,bcjk->ijkabc', t1, bcjk)
    t3d += numpy.einsum('ai,jkbc->ijkabc', eris.fock[nocc:,:nocc], t2)
    t3d = t3d - t3d.transpose(0,1,2,4,3,5) - t3d.transpose(0,1,2,5,4,3)
    t3d = t3d - t3d.transpose(1,0,2,3,4,5) - t3d.transpose(2,1,0,3,4,5)
    t3d /= d3

    goo = numpy.einsum('iklabc,jklabc->ij', (t3c+t3d).conj(), t3c) * (1./12)
    gvv = numpy.einsum('ijkacd,ijkbcd->ab', t3c+t3d, t3c.conj()) * (1./12)
    doo[numpy.diag_indices(nocc)] -= goo.diagonal()
    dvv[numpy.diag_indices(nvir)] += gvv.diagonal()
    dvo += numpy.einsum('ijab,ijkabc->ck', t2.conj(), t3c) * (1./4)

    return doo, dov, dvo, dvv

# gamma2 intermediates in Chemist's notation
def _gamma2_intermediates(mycc, t1, t2, l1, l2, eris=None):
    dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = \
            gccsd_rdm._gamma2_intermediates(mycc, t1, t2, l1, l2)
    if eris is None: eris = mycc.ao2mo()

    nocc, nvir = t1.shape
    bcei = numpy.asarray(eris.ovvv).conj().transpose(3,2,1,0)
    majk = numpy.asarray(eris.ooov).conj().transpose(2,3,0,1)
    bcjk = numpy.asarray(eris.oovv).conj().transpose(2,3,0,1)

    mo_e = eris.mo_energy
    eia = mo_e[:nocc,None] - mo_e[nocc:]
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eia, eia, eia)

    t3c =(numpy.einsum('jkae,bcei->ijkabc', t2, bcei) -
          numpy.einsum('imbc,majk->ijkabc', t2, majk))
    t3c = t3c - t3c.transpose(0,1,2,4,3,5) - t3c.transpose(0,1,2,5,4,3)
    t3c = t3c - t3c.transpose(1,0,2,3,4,5) - t3c.transpose(2,1,0,3,4,5)
    t3c /= d3

    t3d = numpy.einsum('ia,bcjk->ijkabc', t1, bcjk)
    t3d += numpy.einsum('ai,jkbc->ijkabc', eris.fock[nocc:,:nocc], t2)
    t3d = t3d - t3d.transpose(0,1,2,4,3,5) - t3d.transpose(0,1,2,5,4,3)
    t3d = t3d - t3d.transpose(1,0,2,3,4,5) - t3d.transpose(2,1,0,3,4,5)
    t3d /= d3

    goovv = numpy.einsum('kc,ijkabc->ijab', t1.conj(), t3c).conj() * (1./4)
    dovov += goovv.transpose(0,2,1,3) - goovv.transpose(0,3,1,2)

    m3 = t3c * 2 + t3d
# *(1/8) instead of (1/4) because ooov appears 4 times in the 2pdm tensor due
# to symmetrization, and its contribution is scaled by 1/2 in Tr(H,2pdm)
    gooov  = numpy.einsum('imbc,ijkabc->jkma', t2, m3.conj()) * (1./8)
    dooov -= gooov.transpose(0,2,1,3) - gooov.transpose(1,2,0,3)

    govvv  = numpy.einsum('jkae,ijkabc->iecb', t2, m3.conj()) * (1./8)
    dovvv += govvv.transpose(0,2,1,3) - govvv.transpose(0,3,1,2)
    return dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov

[docs] def make_rdm1(mycc, t1, t2, l1, l2, eris=None, ao_repr=False): d1 = _gamma1_intermediates(mycc, t1, t2, l1, l2, eris) return gccsd_rdm._make_rdm1(mycc, d1, True, ao_repr=ao_repr)
# rdm2 in Chemist's notation
[docs] def make_rdm2(mycc, t1, t2, l1, l2, eris=None): d1 = _gamma1_intermediates(mycc, t1, t2, l1, l2, eris) d2 = _gamma2_intermediates(mycc, t1, t2, l1, l2, eris) return gccsd_rdm._make_rdm2(mycc, d1, d2, True, True)
if __name__ == '__main__': from functools import reduce from pyscf import gto from pyscf import scf from pyscf import ao2mo from pyscf import cc mol = gto.Mole() mol.atom = [ [8 , (0. , 0. , 0.)], [1 , (0. , -.957 , .587)], [1 , (0.2, .757 , .487)]] mol.basis = '631g' mol.build() mf0 = mf = scf.RHF(mol).run(conv_tol=1.) mf = mf.to_ghf() from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda from pyscf.cc import ccsd_t_rdm_slow as ccsd_t_rdm mycc0 = cc.CCSD(mf0) eris0 = mycc0.ao2mo() mycc0.kernel(eris=eris0) t1 = mycc0.t1 t2 = mycc0.t2 imds = ccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) l1, l2 = ccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) dm1ref = ccsd_t_rdm.make_rdm1(mycc0, t1, t2, l1, l2, eris0) dm2ref = ccsd_t_rdm.make_rdm2(mycc0, t1, t2, l1, l2, eris0) mycc = cc.GCCSD(mf) eris = mycc.ao2mo() t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) gdm1 = make_rdm1(mycc, t1, t2, l1, l2, eris) gdm2 = make_rdm2(mycc, t1, t2, l1, l2, eris) idxa = numpy.where(mycc.mo_coeff.orbspin == 0)[0] idxb = numpy.where(mycc.mo_coeff.orbspin == 1)[0] trdm1 = gdm1[idxa[:,None],idxa] trdm1+= gdm1[idxb[:,None],idxb] trdm2 = gdm2[idxa[:,None,None,None],idxa[:,None,None],idxa[:,None],idxa] trdm2+= gdm2[idxb[:,None,None,None],idxb[:,None,None],idxb[:,None],idxb] dm2ab = gdm2[idxa[:,None,None,None],idxa[:,None,None],idxb[:,None],idxb] trdm2+= dm2ab trdm2+= dm2ab.transpose(2,3,0,1) print(abs(trdm1 - dm1ref).max()) print(abs(trdm2 - dm2ref).max()) mol = gto.Mole() mol.atom = [ [8 , (0. , 0. , 0.)], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]] mol.basis = '631g' mol.spin = 2 mol.charge = 2 mol.build() mf0 = mf = scf.UHF(mol).run(conv_tol=1) mf = scf.addons.convert_to_ghf(mf) from pyscf.cc import uccsd_t_lambda from pyscf.cc import uccsd_t_rdm mycc0 = cc.UCCSD(mf0) eris0 = mycc0.ao2mo() mycc0.kernel(eris=eris0) t1 = mycc0.t1 t2 = mycc0.t2 imds = uccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) l1, l2 = uccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) dm1ref = uccsd_t_rdm.make_rdm1(mycc0, t1, t2, l1, l2, eris0) dm2ref = uccsd_t_rdm.make_rdm2(mycc0, t1, t2, l1, l2, eris0) mycc = cc.GCCSD(mf) eris = mycc.ao2mo() t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) gdm1 = make_rdm1(mycc, t1, t2, l1, l2, eris) gdm2 = make_rdm2(mycc, t1, t2, l1, l2, eris) idxa = numpy.where(mycc.mo_coeff.orbspin == 0)[0] idxb = numpy.where(mycc.mo_coeff.orbspin == 1)[0] print(abs(dm1ref[0] - gdm1[idxa[:,None],idxa]).max()) print(abs(dm1ref[1] - gdm1[idxb[:,None],idxb]).max()) print(abs(dm2ref[0] - gdm2[idxa[:,None,None,None],idxa[:,None,None],idxa[:,None],idxa]).max()) print(abs(dm2ref[1] - gdm2[idxa[:,None,None,None],idxa[:,None,None],idxb[:,None],idxb]).max()) print(abs(dm2ref[2] - gdm2[idxb[:,None,None,None],idxb[:,None,None],idxb[:,None],idxb]).max()) ecc, t1, t2 = mycc.kernel(eris=eris) e3ref = mycc.e_tot + mycc.ccsd_t() l1, l2 = mycc.solve_lambda(eris=eris) dm1 = make_rdm1(mycc, t1, t2, l1, l2, eris=eris) dm2 = make_rdm2(mycc, t1, t2, l1, l2, eris=eris) nao = mol.nao_nr() mo_a = mf.mo_coeff[:nao] mo_b = mf.mo_coeff[nao:] nmo = mo_a.shape[1] eri = ao2mo.kernel(mf._eri, mo_a+mo_b, compact=False).reshape([nmo]*4) orbspin = mf.mo_coeff.orbspin sym_forbid = (orbspin[:,None] != orbspin) eri[sym_forbid,:,:] = 0 eri[:,:,sym_forbid] = 0 hcore = scf.RHF(mol).get_hcore() h1 = reduce(numpy.dot, (mo_a.T.conj(), hcore, mo_a)) h1+= reduce(numpy.dot, (mo_b.T.conj(), hcore, mo_b)) e3 = numpy.einsum('ij,ji', h1, dm1) e3+= numpy.einsum('ijkl,ijkl', eri, dm2) * .5 e3+= mol.energy_nuc() print(e3 - e3ref)