Source code for pyscf.pbc.cc.kuccsd_rdm

#!/usr/bin/env python
# Copyright 2017-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.

import numpy
from pyscf import lib
import itertools

#einsum = numpy.einsum
einsum = lib.einsum

#TODO: optimize memory use

def _gamma1_intermediates(cc, t1, t2, l1=None, l2=None):
    if l1 is None:
        l1 = [amp.conj() for amp in t1]
    if l2 is None:
        l2 = [amp.conj() for amp in t2]
    t1a, t1b = t1
    t2aa, t2ab, t2bb = t2

    l1a, l1b = l1
    l2aa, l2ab, l2bb = l2
    nkpts, nocca, nvira = t1a.shape
    _, noccb, nvirb = t1b.shape

    kconserv = cc.khelper.kconserv

    dooa  = -einsum('xie,xje->xij', l1a, t1a)
    dooa -=  einsum('xyzimef,xyzjmef->xij', l2ab, t2ab)
    dooa -=  einsum('xyzimef,xyzjmef->xij', l2aa, t2aa) * .5
    doob  = -einsum('xie,xje->xij', l1b, t1b)
    doob -=  einsum('yxzmief,yxzmjef->xij', l2ab, t2ab)
    doob -=  einsum('xyzimef,xyzjmef->xij', l2bb, t2bb) * .5

    dvva  = einsum('xma,xmb->xab', t1a, l1a)
    dvva += einsum('xyzmnae,xyzmnbe->zab', t2ab, l2ab)
    dvva += einsum('xyzmnae,xyzmnbe->zab', t2aa, l2aa) * .5
    dvvb  = einsum('xma,xmb->xab', t1b, l1b)

    for km, kn, ke in itertools.product(range(nkpts),repeat=3):
        ka = kconserv[km,ke,kn]
        dvvb[ka] += einsum('mnea,mneb->ab', t2ab[km,kn,ke], l2ab[km,kn,ke])
    dvvb += einsum('xyzmnae,xyzmnbe->zab', t2bb, l2bb) * .5

    xt1a  = einsum('xyzmnef,xyzinef->xmi', l2aa, t2aa) * .5
    xt1a += einsum('xyzmnef,xyzinef->xmi', l2ab, t2ab)
    xt2a  = einsum('xyzmnaf,xyzmnef->zae', t2aa, l2aa) * .5
    xt2a += einsum('xyzmnaf,xyzmnef->zae', t2ab, l2ab)
    xt2a += einsum('xma,xme->xae', t1a, l1a)

    dvoa = t1a.copy().transpose(0,2,1)
    for ka, km in itertools.product(range(nkpts),repeat=2):
        dvoa[ka] += einsum('imae, me->ai', t2aa[ka,km,ka], l1a[km])
        dvoa[ka] += einsum('imae, me->ai', t2ab[ka,km,ka], l1b[km])
    dvoa -= einsum('xmi,xma->xai', xt1a, t1a)
    dvoa -= einsum('xie,xae->xai', t1a, xt2a)

    xt1b  = einsum('xyzmnef,xyzinef->xmi', l2bb, t2bb) * .5
    xt1b += einsum('xyznmef,xyznief->ymi', l2ab, t2ab)
    xt2b  = einsum('xyzmnaf,xyzmnef->zae', t2bb, l2bb) * .5
    for km, kn, kf in itertools.product(range(nkpts),repeat=3):
        ka = kconserv[km,kf,kn]
        xt2b[ka] += einsum('mnfa,mnfe->ae', t2ab[km,kn,kf], l2ab[km,kn,kf])
    xt2b += einsum('xma,xme->xae', t1b, l1b)

    dvob = t1b.copy().transpose(0,2,1)
    for ka, km in itertools.product(range(nkpts),repeat=2):
        dvob[ka] += einsum('imae,me->ai', t2bb[ka,km,ka], l1b[km])
        dvob[ka] += einsum('miea,me->ai', t2ab[km,ka,km], l1a[km])
    dvob -= einsum('xmi, xma->xai',xt1b, t1b)
    dvob -= einsum('xie,xae->xai', t1b, xt2b)

    dova = l1a
    dovb = l1b

    return ((dooa, doob), (dova, dovb), (dvoa, dvob), (dvva, dvvb))


[docs] def make_rdm1(mycc, t1, t2, l1=None, l2=None, ao_repr=False): r''' One-particle spin density matrices dm1a, dm1b in MO basis (the occupied-virtual blocks due to the orbital response contribution are not included). dm1a[p,q] = <q_alpha^\dagger p_alpha> dm1b[p,q] = <q_beta^\dagger p_beta> The convention of 1-pdm is based on McWeeney's book, Eq (5.4.20). ''' d1 = _gamma1_intermediates(mycc, t1, t2, l1, l2) return _make_rdm1(mycc, d1, with_frozen=True, ao_repr=ao_repr)
def _make_rdm1(mycc, d1, with_frozen=True, ao_repr=False): doo, dOO = d1[0] dov, dOV = d1[1] dvo, dVO = d1[2] dvv, dVV = d1[3] nkpts, nocca, nvira = dov.shape _, noccb, nvirb = dOV.shape nmoa = nocca + nvira nmob = noccb + nvirb dtype = numpy.result_type(doo, dOO, dov, dOV, dvo, dVO, dvv, dVV) dm1a = numpy.empty((nkpts,nmoa,nmoa), dtype=dtype) dm1a[:,:nocca,:nocca] = doo + doo.conj().transpose(0,2,1) dm1a[:,:nocca,nocca:] = dov + dvo.conj().transpose(0,2,1) dm1a[:,nocca:,:nocca] = dm1a[:,:nocca,nocca:].conj().transpose(0,2,1) dm1a[:,nocca:,nocca:] = dvv + dvv.conj().transpose(0,2,1) dm1a *= .5 for k in range(nkpts): dm1a[k][numpy.diag_indices(nocca)] +=1 dm1b = numpy.empty((nkpts,nmob,nmob), dtype=dtype) dm1b[:,:noccb,:noccb] = dOO + dOO.conj().transpose(0,2,1) dm1b[:,:noccb,noccb:] = dOV + dVO.conj().transpose(0,2,1) dm1b[:,noccb:,:noccb] = dm1b[:,:noccb,noccb:].conj().transpose(0,2,1) dm1b[:,noccb:,noccb:] = dVV + dVV.conj().transpose(0,2,1) dm1b *= .5 for k in range(nkpts): dm1b[k][numpy.diag_indices(noccb)] +=1 if with_frozen and mycc.frozen is not None: raise NotImplementedError _, nmoa = mycc.mo_occ[0].size _, nmob = mycc.mo_occ[1].size nocca = numpy.count_nonzero(mycc.mo_occ[0] > 0) noccb = numpy.count_nonzero(mycc.mo_occ[1] > 0) rdm1a = numpy.zeros((nkpts,nmoa,nmoa), dtype=dm1a.dtype) rdm1b = numpy.zeros((nkpts,nmob,nmob), dtype=dm1b.dtype) rdm1a[:,numpy.diag_indices(nocca)] = 1 rdm1b[:,numpy.diag_indices(noccb)] = 1 moidx = mycc.get_frozen_mask() moidxa = numpy.where(moidx[0])[0] moidxb = numpy.where(moidx[1])[0] rdm1a[moidxa[:,None],moidxa] = dm1a rdm1b[moidxb[:,None],moidxb] = dm1b dm1a = rdm1a dm1b = rdm1b if ao_repr: mo_a, mo_b = mycc.mo_coeff dm1a = lib.einsum('xpi,xij,xqj->xpq', mo_a, dm1a, mo_a.conj()) dm1b = lib.einsum('xpi,xij,xqj->xpq', mo_b, dm1b, mo_b.conj()) return dm1a, dm1b