#!/usr/bin/env python
# Copyright 2014-2020 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>
#
'''
Restricted CCSD lambda equation solver which supports both real and complex
integrals. This code is slower than the pyscf.cc.ccsd_lambda implementation.
Note MO integrals are treated in chemist's notation
'''
import numpy as np
from pyscf import lib
from pyscf.lib import logger
from pyscf.cc import _ccsd
from pyscf.cc import ccsd_lambda
einsum = lib.einsum
[docs]
def kernel(mycc, eris=None, t1=None, t2=None, l1=None, l2=None,
max_cycle=50, tol=1e-8, verbose=logger.INFO):
if eris is None: eris = mycc.ao2mo()
return ccsd_lambda.kernel(mycc, eris, t1, t2, l1, l2, max_cycle, tol,
verbose, make_intermediates, update_lambda)
# update L1, L2
[docs]
def update_lambda(mycc, t1, t2, l1, l2, eris, imds):
time0 = logger.process_clock(), logger.perf_counter()
log = logger.Logger(mycc.stdout, mycc.verbose)
nocc, nvir = t1.shape
l1new = np.zeros_like(l1)
foo = eris.fock[:nocc,:nocc]
fov = eris.fock[:nocc,nocc:]
fvv = eris.fock[nocc:,nocc:]
tau = _ccsd.make_tau(t2, t1, t1)
theta = t2*2 - t2.transpose(0,1,3,2)
mvv = lib.einsum('klca,klcb->ba', l2, theta)
moo = lib.einsum('kicd,kjcd->ij', l2, theta)
mvv1 = lib.einsum('jc,jb->bc', l1, t1) + mvv
moo1 = lib.einsum('ic,kc->ik', l1, t1) + moo
# m3 = einsum('ijab,acbd->ijcd', l2, vvvv)
# = einsum('ijab,cadb->ijcd', l2.conj(), vvvv).conj()
m3 = mycc._add_vvvv(None, l2.conj(), eris, with_ovvv=False).conj()
m3 += lib.einsum('klab,ikjl->ijab', l2, imds.woooo)
m3 *= .5
ovov = np.asarray(eris.ovov)
l2tau = np.einsum('ijcd,klcd->ijkl', l2, tau)
m3 += np.einsum('kalb,ijkl->ijab', ovov, l2tau) * .5
l2tau = None
l2new = ovov.transpose(0,2,1,3) * .5
l2new += lib.einsum('ijac,cb->ijab', l2, imds.v1)
l2new -= lib.einsum('ikab,jk->ijab', l2, imds.v2)
l2new -= lib.einsum('ca,icjb->ijab', mvv1, ovov)
l2new -= lib.einsum('ik,kajb->ijab', moo1, ovov)
ovov = ovov * 2 - ovov.transpose(0,3,2,1)
l1new -= np.einsum('ik,ka->ia', moo, imds.v4)
l1new -= np.einsum('ca,ic->ia', mvv, imds.v4)
l2new += np.einsum('ia,jb->ijab', l1, imds.v4)
tmp = t1 + np.einsum('kc,kjcb->jb', l1, theta)
tmp -= lib.einsum('bd,jd->jb', mvv1, t1)
tmp -= lib.einsum('lj,lb->jb', moo, t1)
l1new += np.einsum('jbia,jb->ia', ovov, tmp)
ovov = tmp = None
ovvv = np.asarray(eris.get_ovvv())
l1new += np.einsum('iacb,bc->ia', ovvv, mvv1) * 2
l1new -= np.einsum('ibca,bc->ia', ovvv, mvv1)
l2new += lib.einsum('ic,jbca->jiba', l1, ovvv)
l2t1 = np.einsum('ijcd,kd->ijck', l2, t1)
m3 -= np.einsum('kbca,ijck->ijab', ovvv, l2t1)
l2t1 = ovvv = None
l2new += m3
l1new += np.einsum('ijab,jb->ia', m3, t1) * 2
l1new += np.einsum('jiba,jb->ia', m3, t1) * 2
l1new -= np.einsum('ijba,jb->ia', m3, t1)
l1new -= np.einsum('jiab,jb->ia', m3, t1)
ovoo = np.asarray(eris.ovoo)
l1new -= np.einsum('iajk,kj->ia', ovoo, moo1) * 2
l1new += np.einsum('jaik,kj->ia', ovoo, moo1)
l2new -= lib.einsum('ka,jbik->ijab', l1, ovoo)
ovoo = None
l2theta = l2*2 - l2.transpose(0,1,3,2)
l2new += lib.einsum('ikac,jbck->ijab', l2theta, imds.wovvo) * .5
tmp = lib.einsum('ikca,jbck->ijab', l2, imds.woVVo)
l2new += tmp * .5
l2new += tmp.transpose(1,0,2,3)
l2theta = None
l1new += fov
l1new += lib.einsum('ib,ba->ia', l1, imds.v1)
l1new -= lib.einsum('ja,ij->ia', l1, imds.v2)
l1new += np.einsum('jb,iabj->ia', l1, eris.ovvo) * 2
l1new -= np.einsum('jb,ijba->ia', l1, eris.oovv)
l1new -= lib.einsum('ijbc,bacj->ia', l2, imds.wvvvo)
l1new -= lib.einsum('kjca,ijck->ia', l2, imds.woovo)
l1new += np.einsum('ijab,bj->ia', l2, imds.w3) * 2
l1new -= np.einsum('ijba,bj->ia', l2, imds.w3)
eia = lib.direct_sum('i-j->ij', foo.diagonal(), fvv.diagonal() + mycc.level_shift)
l1new /= eia
l1new += l1
l2new = l2new + l2new.transpose(1,0,3,2)
l2new /= lib.direct_sum('ia+jb->ijab', eia, eia)
l2new += l2
time0 = log.timer_debug1('update l1 l2', *time0)
return l1new, l2new
if __name__ == '__main__':
from functools import reduce
from pyscf import gto
from pyscf import scf
from pyscf import ao2mo
from pyscf.cc import rccsd
mol = gto.Mole()
mol.verbose = 0
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.RHF(mol)
mf.conv_tol = 1e-16
mf.scf()
mcc = rccsd.RCCSD(mf)
mcc.conv_tol = 1e-12
ecc, t1, t2 = mcc.kernel()
eris = mcc.ao2mo()
l1, l2 = mcc.solve_lambda(t1, t2, eris=eris)
print(np.linalg.norm(l1)-0.0132626841292)
print(np.linalg.norm(l2)-0.212575609057)
from pyscf.cc import ccsd_rdm
dm1 = ccsd_rdm.make_rdm1(mcc, t1, t2, l1, l2)
dm2 = ccsd_rdm.make_rdm2(mcc, t1, t2, l1, l2)
h1 = reduce(np.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff))
nmo = h1.shape[0]
eri = ao2mo.full(mf._eri, mf.mo_coeff)
eri = ao2mo.restore(1, eri, nmo).reshape((nmo,)*4)
e1 = np.einsum('pq,pq', h1, dm1)
e2 = np.einsum('pqrs,pqrs', eri, dm2) * .5
print(e1+e2+mol.energy_nuc() - mf.e_tot - ecc)