#!/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>
# Jun Yang
#
import numpy
from pyscf import lib
#einsum = numpy.einsum
einsum = lib.einsum
def _gamma1_intermediates(mycc, t1, t2, l1, l2):
doo =-einsum('ie,je->ij', l1, t1)
doo -= einsum('imef,jmef->ij', l2, t2) * .5
dvv = einsum('ma,mb->ab', t1, l1)
dvv += einsum('mnea,mneb->ab', t2, l2) * .5
xt1 = einsum('mnef,inef->mi', l2, t2) * .5
xt2 = einsum('mnfa,mnfe->ae', t2, l2) * .5
xt2 += einsum('ma,me->ae', t1, l1)
dvo = einsum('imae,me->ai', t2, l1)
dvo -= einsum('mi,ma->ai', xt1, t1)
dvo -= einsum('ie,ae->ai', t1, xt2)
dvo += t1.T
dov = l1
return doo, dov, dvo, dvv
# gamma2 intermediates in Chemist's notation
# When computing intermediates, the convention
# dm2[q,p,s,r] = <p^\dagger r^\dagger s q> is assumed in this function.
# It changes to dm2[p,q,r,s] = <p^\dagger r^\dagger s q> in _make_rdm2
def _gamma2_intermediates(mycc, t1, t2, l1, l2):
tau = t2 + einsum('ia,jb->ijab', t1, t1) * 2
miajb = einsum('ikac,kjcb->iajb', l2, t2)
goovv = 0.25 * (l2.conj() + tau)
tmp = einsum('kc,kica->ia', l1, t2)
goovv += einsum('ia,jb->ijab', tmp, t1)
tmp = einsum('kc,kb->cb', l1, t1)
goovv += einsum('cb,ijca->ijab', tmp, t2) * .5
tmp = einsum('kc,jc->kj', l1, t1)
goovv += einsum('kiab,kj->ijab', tau, tmp) * .5
tmp = numpy.einsum('ldjd->lj', miajb)
goovv -= einsum('lj,liba->ijab', tmp, tau) * .25
tmp = numpy.einsum('ldlb->db', miajb)
goovv -= einsum('db,jida->ijab', tmp, tau) * .25
goovv -= einsum('ldia,ljbd->ijab', miajb, tau) * .5
tmp = einsum('klcd,ijcd->ijkl', l2, tau) * .25**2
goovv += einsum('ijkl,klab->ijab', tmp, tau)
goovv = goovv.conj()
gvvvv = einsum('ijab,ijcd->abcd', tau, l2) * 0.125
goooo = einsum('klab,ijab->klij', l2, tau) * 0.125
gooov = einsum('jkba,ib->jkia', tau, l1) * -0.25
gooov += einsum('iljk,la->jkia', goooo, t1)
tmp = numpy.einsum('icjc->ij', miajb) * .25
gooov -= einsum('ij,ka->jkia', tmp, t1)
gooov += einsum('icja,kc->jkia', miajb, t1) * .5
gooov = gooov.conj()
gooov += einsum('jkab,ib->jkia', l2, t1) * .25
govvo = einsum('ia,jb->ibaj', l1, t1)
govvo += numpy.einsum('iajb->ibaj', miajb)
govvo -= einsum('ikac,jc,kb->ibaj', l2, t1, t1)
govvv = einsum('ja,ijcb->iacb', l1, tau) * .25
govvv += einsum('bcad,id->iabc', gvvvv, t1)
tmp = numpy.einsum('kakb->ab', miajb) * .25
govvv += einsum('ab,ic->iacb', tmp, t1)
govvv += einsum('kaib,kc->iabc', miajb, t1) * .5
govvv = govvv.conj()
govvv += einsum('ijbc,ja->iabc', l2, t1) * .25
dovov = goovv.transpose(0,2,1,3) - goovv.transpose(0,3,1,2)
dvvvv = gvvvv.transpose(0,2,1,3) - gvvvv.transpose(0,3,1,2)
doooo = goooo.transpose(0,2,1,3) - goooo.transpose(0,3,1,2)
dovvv = govvv.transpose(0,2,1,3) - govvv.transpose(0,3,1,2)
dooov = gooov.transpose(0,2,1,3) - gooov.transpose(1,2,0,3)
dovvo = govvo.transpose(0,2,1,3)
dovov =(dovov + dovov.transpose(2,3,0,1)) * .5
dvvvv = dvvvv + dvvvv.transpose(1,0,3,2).conj()
doooo = doooo + doooo.transpose(1,0,3,2).conj()
dovvo =(dovvo + dovvo.transpose(3,2,1,0).conj()) * .5
doovv = None # = -dovvo.transpose(0,3,2,1)
dvvov = None
return (dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov)
[docs]
def make_rdm1(mycc, t1, t2, l1, l2, ao_repr=False, with_frozen=True, with_mf=True):
r'''
One-particle density matrix in the molecular spin-orbital representation
(the occupied-virtual blocks from the orbital response contribution are
not included).
dm1[p,q] = <q^\dagger p> (p,q are spin-orbitals)
The convention of 1-pdm is based on McWeeney's book, Eq (5.4.20).
The contraction between 1-particle Hamiltonian and rdm1 is
E = einsum('pq,qp', h1, rdm1)
'''
d1 = _gamma1_intermediates(mycc, t1, t2, l1, l2)
return _make_rdm1(mycc, d1, with_frozen=with_frozen, ao_repr=ao_repr,
with_mf=with_mf)
[docs]
def make_rdm2(mycc, t1, t2, l1, l2, ao_repr=False, with_frozen=True, with_dm1=True):
r'''
Two-particle density matrix in the molecular spin-orbital representation
dm2[p,q,r,s] = <p^\dagger r^\dagger s q>
where p,q,r,s are spin-orbitals. p,q correspond to one particle and r,s
correspond to another particle. The contraction between ERIs (in
Chemist's notation) and rdm2 is
E = einsum('pqrs,pqrs', eri, rdm2)
'''
d1 = _gamma1_intermediates(mycc, t1, t2, l1, l2)
d2 = _gamma2_intermediates(mycc, t1, t2, l1, l2)
return _make_rdm2(mycc, d1, d2, with_dm1=with_dm1, with_frozen=with_frozen,
ao_repr=ao_repr)
def _make_rdm1(mycc, d1, with_frozen=True, ao_repr=False, with_mf=True):
r'''
One-particle density matrix in the molecular spin-orbital representation
(the occupied-virtual blocks from the orbital response contribution are
not included).
dm1[p,q] = <q^\dagger p> (p,q are spin-orbitals)
The convention of 1-pdm is based on McWeeney's book, Eq (5.4.20).
The contraction between 1-particle Hamiltonian and rdm1 is
E = einsum('pq,qp', h1, rdm1)
'''
doo, dov, dvo, dvv = d1
nocc, nvir = dov.shape
nmo = nocc + nvir
dm1 = numpy.empty((nmo,nmo), dtype=doo.dtype)
dm1[:nocc,:nocc] = doo + doo.conj().T
dm1[:nocc,nocc:] = dov + dvo.conj().T
dm1[nocc:,:nocc] = dm1[:nocc,nocc:].conj().T
dm1[nocc:,nocc:] = dvv + dvv.conj().T
dm1 *= .5
if with_mf:
dm1[numpy.diag_indices(nocc)] += 1
if with_frozen and mycc.frozen is not None:
nmo = mycc.mo_occ.size
nocc = numpy.count_nonzero(mycc.mo_occ > 0)
rdm1 = numpy.zeros((nmo,nmo), dtype=dm1.dtype)
if with_mf:
rdm1[numpy.diag_indices(nocc)] = 1
moidx = numpy.where(mycc.get_frozen_mask())[0]
rdm1[moidx[:,None],moidx] = dm1
dm1 = rdm1
if ao_repr:
mo = mycc.mo_coeff
dm1 = lib.einsum('pi,ij,qj->pq', mo, dm1, mo.conj())
return dm1
def _make_rdm2(mycc, d1, d2, with_dm1=True, with_frozen=True, ao_repr=False):
r'''
dm2[p,q,r,s] = <p^\dagger r^\dagger s q>
Note the contraction between ERIs (in Chemist's notation) and rdm2 is
E = einsum('pqrs,pqrs', eri, rdm2)
'''
dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2
nocc, nvir = dovov.shape[:2]
nmo = nocc + nvir
dm2 = numpy.empty((nmo,nmo,nmo,nmo), dtype=doooo.dtype)
dovov = numpy.asarray(dovov)
dm2[:nocc,nocc:,:nocc,nocc:] = dovov
dm2[nocc:,:nocc,nocc:,:nocc] = dm2[:nocc,nocc:,:nocc,nocc:].transpose(1,0,3,2).conj()
dovov = None
dovvo = numpy.asarray(dovvo)
dm2[:nocc,:nocc,nocc:,nocc:] =-dovvo.transpose(0,3,2,1)
dm2[nocc:,nocc:,:nocc,:nocc] =-dovvo.transpose(2,1,0,3)
dm2[:nocc,nocc:,nocc:,:nocc] = dovvo
dm2[nocc:,:nocc,:nocc,nocc:] = dovvo.transpose(1,0,3,2).conj()
dovvo = None
dm2[nocc:,nocc:,nocc:,nocc:] = dvvvv
dm2[:nocc,:nocc,:nocc,:nocc] = doooo
dovvv = numpy.asarray(dovvv)
dm2[:nocc,nocc:,nocc:,nocc:] = dovvv
dm2[nocc:,nocc:,:nocc,nocc:] = dovvv.transpose(2,3,0,1)
dm2[nocc:,nocc:,nocc:,:nocc] = dovvv.transpose(3,2,1,0).conj()
dm2[nocc:,:nocc,nocc:,nocc:] = dovvv.transpose(1,0,3,2).conj()
dovvv = None
dooov = numpy.asarray(dooov)
dm2[:nocc,:nocc,:nocc,nocc:] = dooov
dm2[:nocc,nocc:,:nocc,:nocc] = dooov.transpose(2,3,0,1)
dm2[:nocc,:nocc,nocc:,:nocc] = dooov.transpose(1,0,3,2).conj()
dm2[nocc:,:nocc,:nocc,:nocc] = dooov.transpose(3,2,1,0).conj()
if with_frozen and mycc.frozen is not None:
nmo, nmo0 = mycc.mo_occ.size, nmo
nocc = numpy.count_nonzero(mycc.mo_occ > 0)
rdm2 = numpy.zeros((nmo,nmo,nmo,nmo), dtype=dm2.dtype)
moidx = numpy.where(mycc.get_frozen_mask())[0]
idx = (moidx.reshape(-1,1) * nmo + moidx).ravel()
lib.takebak_2d(rdm2.reshape(nmo**2,nmo**2),
dm2.reshape(nmo0**2,nmo0**2), idx, idx)
dm2 = rdm2
if with_dm1:
dm1 = _make_rdm1(mycc, d1, with_frozen)
dm1[numpy.diag_indices(nocc)] -= 1
for i in range(nocc):
# Be careful with the convention of dm1 and dm2.transpose
# at the end
dm2[i,i,:,:] += dm1
dm2[:,:,i,i] += dm1
dm2[:,i,i,:] -= dm1
dm2[i,:,:,i] -= dm1.T
for i in range(nocc):
for j in range(nocc):
dm2[i,i,j,j] += 1
dm2[i,j,j,i] -= 1
# dm2 was computed as dm2[p,q,r,s] = < p^\dagger r^\dagger s q > in the
# above. Transposing it so that it be contracted with ERIs (in Chemist's
# notation):
# E = einsum('pqrs,pqrs', eri, rdm2)
dm2 = dm2.transpose(1,0,3,2)
if ao_repr:
from pyscf.cc import ccsd_rdm
dm2 = ccsd_rdm._rdm2_mo2ao(dm2, mycc.mo_coeff)
return dm2
if __name__ == '__main__':
from functools import reduce
from pyscf import gto
from pyscf import scf
from pyscf import ao2mo
from pyscf.cc import gccsd
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.build()
mf = scf.UHF(mol).run(conv_tol=1.)
mf = scf.addons.convert_to_ghf(mf)
mycc = gccsd.GCCSD(mf)
ecc, t1, t2 = mycc.kernel()
l1, l2 = mycc.solve_lambda()
dm1 = make_rdm1(mycc, t1, t2, l1, l2)
dm2 = make_rdm2(mycc, t1, t2, l1, l2)
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))
e1 = numpy.einsum('ij,ji', h1, dm1)
e1+= numpy.einsum('ijkl,ijkl', eri, dm2) * .5
e1+= mol.energy_nuc()
print(e1 - mycc.e_tot)
#TODO: test 1pdm, 2pdm against FCI
