#!/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: Timothy Berkelbach <tim.berkelbach@gmail.com>
#
'''
Restricted CCSD implementation which supports both real and complex integrals.
The 4-index integrals are saved on disk entirely (without using any symmetry).
This code is slower than the pyscf.cc.ccsd implementation.
Note MO integrals are treated in chemist's notation
Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004)
'''
import numpy as np
from pyscf import lib
from pyscf import ao2mo
from pyscf.lib import logger
from pyscf.cc import ccsd
from pyscf.cc import rintermediates as imd
from pyscf.mp import mp2
from pyscf import __config__
BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4)
MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000)
[docs]
def update_amps(cc, t1, t2, eris):
# Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004) Eqs.(35)-(36)
assert (isinstance(eris, ccsd._ChemistsERIs))
nocc, nvir = t1.shape
fock = eris.fock
mo_e_o = eris.mo_energy[:nocc]
mo_e_v = eris.mo_energy[nocc:] + cc.level_shift
fov = fock[:nocc,nocc:].copy()
foo = fock[:nocc,:nocc].copy()
fvv = fock[nocc:,nocc:].copy()
Foo = imd.cc_Foo(t1,t2,eris)
Fvv = imd.cc_Fvv(t1,t2,eris)
Fov = imd.cc_Fov(t1,t2,eris)
# Move energy terms to the other side
Foo[np.diag_indices(nocc)] -= mo_e_o
Fvv[np.diag_indices(nvir)] -= mo_e_v
# T1 equation
t1new =-2*np.einsum('kc,ka,ic->ia', fov, t1, t1)
t1new += np.einsum('ac,ic->ia', Fvv, t1)
t1new += -np.einsum('ki,ka->ia', Foo, t1)
t1new += 2*np.einsum('kc,kica->ia', Fov, t2)
t1new += -np.einsum('kc,ikca->ia', Fov, t2)
t1new += np.einsum('kc,ic,ka->ia', Fov, t1, t1)
t1new += fov.conj()
t1new += 2*np.einsum('kcai,kc->ia', eris.ovvo, t1)
t1new += -np.einsum('kiac,kc->ia', eris.oovv, t1)
eris_ovvv = np.asarray(eris.get_ovvv())
t1new += 2*lib.einsum('kdac,ikcd->ia', eris_ovvv, t2)
t1new += -lib.einsum('kcad,ikcd->ia', eris_ovvv, t2)
t1new += 2*lib.einsum('kdac,kd,ic->ia', eris_ovvv, t1, t1)
t1new += -lib.einsum('kcad,kd,ic->ia', eris_ovvv, t1, t1)
eris_ovoo = np.asarray(eris.ovoo, order='C')
t1new +=-2*lib.einsum('lcki,klac->ia', eris_ovoo, t2)
t1new += lib.einsum('kcli,klac->ia', eris_ovoo, t2)
t1new +=-2*lib.einsum('lcki,lc,ka->ia', eris_ovoo, t1, t1)
t1new += lib.einsum('kcli,lc,ka->ia', eris_ovoo, t1, t1)
# T2 equation
tmp2 = lib.einsum('kibc,ka->abic', eris.oovv, -t1)
tmp2 += np.asarray(eris_ovvv).conj().transpose(1,3,0,2)
tmp = lib.einsum('abic,jc->ijab', tmp2, t1)
t2new = tmp + tmp.transpose(1,0,3,2)
tmp2 = lib.einsum('kcai,jc->akij', eris.ovvo, t1)
tmp2 += eris_ovoo.transpose(1,3,0,2).conj()
tmp = lib.einsum('akij,kb->ijab', tmp2, t1)
t2new -= tmp + tmp.transpose(1,0,3,2)
t2new += np.asarray(eris.ovov).conj().transpose(0,2,1,3)
if cc.cc2:
Woooo2 = np.asarray(eris.oooo).transpose(0,2,1,3).copy()
Woooo2 += lib.einsum('lcki,jc->klij', eris_ovoo, t1)
Woooo2 += lib.einsum('kclj,ic->klij', eris_ovoo, t1)
Woooo2 += lib.einsum('kcld,ic,jd->klij', eris.ovov, t1, t1)
t2new += lib.einsum('klij,ka,lb->ijab', Woooo2, t1, t1)
Wvvvv = lib.einsum('kcbd,ka->abcd', eris_ovvv, -t1)
Wvvvv = Wvvvv + Wvvvv.transpose(1,0,3,2)
Wvvvv += np.asarray(eris.vvvv).transpose(0,2,1,3)
t2new += lib.einsum('abcd,ic,jd->ijab', Wvvvv, t1, t1)
Lvv2 = fvv - np.einsum('kc,ka->ac', fov, t1)
Lvv2 -= np.diag(np.diag(fvv))
tmp = lib.einsum('ac,ijcb->ijab', Lvv2, t2)
t2new += (tmp + tmp.transpose(1,0,3,2))
Loo2 = foo + np.einsum('kc,ic->ki', fov, t1)
Loo2 -= np.diag(np.diag(foo))
tmp = lib.einsum('ki,kjab->ijab', Loo2, t2)
t2new -= (tmp + tmp.transpose(1,0,3,2))
else:
Loo = imd.Loo(t1, t2, eris)
Lvv = imd.Lvv(t1, t2, eris)
Loo[np.diag_indices(nocc)] -= mo_e_o
Lvv[np.diag_indices(nvir)] -= mo_e_v
Woooo = imd.cc_Woooo(t1, t2, eris)
Wvoov = imd.cc_Wvoov(t1, t2, eris)
Wvovo = imd.cc_Wvovo(t1, t2, eris)
Wvvvv = imd.cc_Wvvvv(t1, t2, eris)
tau = t2 + np.einsum('ia,jb->ijab', t1, t1)
t2new += lib.einsum('klij,klab->ijab', Woooo, tau)
t2new += lib.einsum('abcd,ijcd->ijab', Wvvvv, tau)
tmp = lib.einsum('ac,ijcb->ijab', Lvv, t2)
t2new += (tmp + tmp.transpose(1,0,3,2))
tmp = lib.einsum('ki,kjab->ijab', Loo, t2)
t2new -= (tmp + tmp.transpose(1,0,3,2))
tmp = 2*lib.einsum('akic,kjcb->ijab', Wvoov, t2)
tmp -= lib.einsum('akci,kjcb->ijab', Wvovo, t2)
t2new += (tmp + tmp.transpose(1,0,3,2))
tmp = lib.einsum('akic,kjbc->ijab', Wvoov, t2)
t2new -= (tmp + tmp.transpose(1,0,3,2))
tmp = lib.einsum('bkci,kjac->ijab', Wvovo, t2)
t2new -= (tmp + tmp.transpose(1,0,3,2))
eia = mo_e_o[:,None] - mo_e_v
eijab = lib.direct_sum('ia,jb->ijab',eia,eia)
t1new /= eia
t2new /= eijab
return t1new, t2new
[docs]
def energy(cc, t1=None, t2=None, eris=None):
'''RCCSD correlation energy'''
if t1 is None: t1 = cc.t1
if t2 is None: t2 = cc.t2
if eris is None: eris = cc.ao2mo()
nocc, nvir = t1.shape
fock = eris.fock
e = 2*np.einsum('ia,ia', fock[:nocc,nocc:], t1)
tau = np.einsum('ia,jb->ijab',t1,t1)
tau += t2
eris_ovov = np.asarray(eris.ovov)
e += 2*np.einsum('ijab,iajb', tau, eris_ovov)
e += -np.einsum('ijab,ibja', tau, eris_ovov)
if abs(e.imag) > 1e-4:
logger.warn(cc, 'Non-zero imaginary part found in RCCSD energy %s', e)
return e.real
[docs]
class RCCSD(ccsd.CCSD):
'''restricted CCSD with IP-EOM, EA-EOM, EE-EOM, and SF-EOM capabilities
Ground-state CCSD is performed in optimized ccsd.CCSD and EOM is performed here.
'''
[docs]
def kernel(self, t1=None, t2=None, eris=None, mbpt2=False):
return self.ccsd(t1, t2, eris, mbpt2)
[docs]
def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False):
'''Ground-state CCSD.
Kwargs:
mbpt2 : bool
Use one-shot MBPT2 approximation to CCSD.
'''
if mbpt2:
pt = mp2.MP2(self._scf, self.frozen, self.mo_coeff, self.mo_occ)
self.e_corr, self.t2 = pt.kernel(eris=eris)
nocc, nvir = self.t2.shape[1:3]
self.t1 = np.zeros((nocc,nvir))
return self.e_corr, self.t1, self.t2
if eris is None:
eris = self.ao2mo(self.mo_coeff)
return ccsd.CCSDBase.ccsd(self, t1, t2, eris)
[docs]
def ao2mo(self, mo_coeff=None):
nmo = self.nmo
nao = self.mo_coeff.shape[0]
nmo_pair = nmo * (nmo+1) // 2
nao_pair = nao * (nao+1) // 2
mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6
mem_now = lib.current_memory()[0]
if (self._scf._eri is not None and
(mem_incore+mem_now < self.max_memory) or self.mol.incore_anyway):
return _make_eris_incore(self, mo_coeff)
elif getattr(self._scf, 'with_df', None):
logger.warn(self, 'CCSD detected DF being used in the HF object. '
'MO integrals are computed based on the DF 3-index tensors.\n'
'It\'s recommended to use dfccsd.CCSD for the '
'DF-CCSD calculations')
raise NotImplementedError
#return _make_df_eris_outcore(self, mo_coeff)
else:
return _make_eris_outcore(self, mo_coeff)
energy = energy
update_amps = update_amps
[docs]
def solve_lambda(self, t1=None, t2=None, l1=None, l2=None,
eris=None):
from pyscf.cc import rccsd_lambda
if t1 is None: t1 = self.t1
if t2 is None: t2 = self.t2
if eris is None: eris = self.ao2mo(self.mo_coeff)
self.converged_lambda, self.l1, self.l2 = \
rccsd_lambda.kernel(self, eris, t1, t2, l1, l2,
max_cycle=self.max_cycle,
tol=self.conv_tol_normt,
verbose=self.verbose)
return self.l1, self.l2
class _ChemistsERIs(ccsd._ChemistsERIs):
def get_ovvv(self, *slices):
'''To access a subblock of ovvv tensor'''
if slices:
return self.ovvv[slices]
else:
return self.ovvv
def _make_eris_incore(mycc, mo_coeff=None, ao2mofn=None):
cput0 = (logger.process_clock(), logger.perf_counter())
eris = _ChemistsERIs()
eris._common_init_(mycc, mo_coeff)
nocc = eris.nocc
nmo = eris.fock.shape[0]
if callable(ao2mofn):
eri1 = ao2mofn(eris.mo_coeff).reshape([nmo]*4)
else:
eri1 = ao2mo.incore.full(mycc._scf._eri, eris.mo_coeff)
eri1 = ao2mo.restore(1, eri1, nmo)
eris.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy()
eris.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy()
eris.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy()
eris.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy()
eris.ovvo = eri1[:nocc,nocc:,nocc:,:nocc].copy()
eris.ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy()
eris.vvvv = eri1[nocc:,nocc:,nocc:,nocc:].copy()
logger.timer(mycc, 'CCSD integral transformation', *cput0)
return eris
def _make_eris_outcore(mycc, mo_coeff=None):
from pyscf.scf.hf import RHF
assert isinstance(mycc._scf, RHF)
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(mycc.stdout, mycc.verbose)
eris = _ChemistsERIs()
eris._common_init_(mycc, mo_coeff)
mol = mycc.mol
mo_coeff = eris.mo_coeff
nocc = eris.nocc
nao, nmo = mo_coeff.shape
nvir = nmo - nocc
eris.feri1 = lib.H5TmpFile()
eris.oooo = eris.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8')
eris.ovoo = eris.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), 'f8', chunks=(nocc,1,nocc,nocc))
eris.ovov = eris.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), 'f8', chunks=(nocc,1,nocc,nvir))
eris.ovvo = eris.feri1.create_dataset('ovvo', (nocc,nvir,nvir,nocc), 'f8', chunks=(nocc,1,nvir,nocc))
eris.ovvv = eris.feri1.create_dataset('ovvv', (nocc,nvir,nvir,nvir), 'f8')
eris.oovv = eris.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), 'f8', chunks=(nocc,nocc,1,nvir))
eris.vvvv = eris.feri1.create_dataset('vvvv', (nvir,nvir,nvir,nvir), 'f8')
max_memory = max(MEMORYMIN, mycc.max_memory-lib.current_memory()[0])
ftmp = lib.H5TmpFile()
ao2mo.full(mol, mo_coeff, ftmp, max_memory=max_memory, verbose=log)
eri = ftmp['eri_mo']
nocc_pair = nocc*(nocc+1)//2
tril2sq = lib.square_mat_in_trilu_indices(nmo)
oo = eri[:nocc_pair]
eris.oooo[:] = ao2mo.restore(1, oo[:,:nocc_pair], nocc)
oovv = lib.take_2d(oo, tril2sq[:nocc,:nocc].ravel(), tril2sq[nocc:,nocc:].ravel())
eris.oovv[:] = oovv.reshape(nocc,nocc,nvir,nvir)
oo = oovv = None
tril2sq = lib.square_mat_in_trilu_indices(nmo)
blksize = min(nvir, max(BLKMIN, int(max_memory*1e6/8/nmo**3/2)))
for p0, p1 in lib.prange(0, nvir, blksize):
q0, q1 = p0+nocc, p1+nocc
off0 = q0*(q0+1)//2
off1 = q1*(q1+1)//2
buf = lib.unpack_tril(eri[off0:off1])
tmp = buf[ tril2sq[q0:q1,:nocc] - off0 ]
eris.ovoo[:,p0:p1] = tmp[:,:,:nocc,:nocc].transpose(1,0,2,3)
eris.ovvo[:,p0:p1] = tmp[:,:,nocc:,:nocc].transpose(1,0,2,3)
eris.ovov[:,p0:p1] = tmp[:,:,:nocc,nocc:].transpose(1,0,2,3)
eris.ovvv[:,p0:p1] = tmp[:,:,nocc:,nocc:].transpose(1,0,2,3)
tmp = buf[ tril2sq[q0:q1,nocc:q1] - off0 ]
eris.vvvv[p0:p1,:p1] = tmp[:,:,nocc:,nocc:]
if p0 > 0:
eris.vvvv[:p0,p0:p1] = tmp[:,:p0,nocc:,nocc:].transpose(1,0,2,3)
buf = tmp = None
log.timer('CCSD integral transformation', *cput0)
return eris
if __name__ == '__main__':
from pyscf import scf
from pyscf import gto
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 = 'cc-pvdz'
#mol.basis = '3-21G'
mol.verbose = 0
mol.spin = 0
mol.build()
mf = scf.RHF(mol).run(conv_tol=1e-14)
mycc = RCCSD(mf)
mycc.max_memory = 0
eris = mycc.ao2mo()
emp2, t1, t2 = mycc.init_amps(eris)
print(lib.finger(t2) - 0.044540097905897198)
np.random.seed(1)
t1 = np.random.random(t1.shape)*.1
t2 = np.random.random(t2.shape)*.1
t2 = t2 + t2.transpose(1,0,3,2)
t1, t2 = update_amps(mycc, t1, t2, eris)
print(lib.finger(t1) - 0.25118555558133576)
print(lib.finger(t2) - 0.02352137419932243)
ecc, t1, t2 = mycc.kernel()
print(ecc - -0.21334326214236796)
e, v = mycc.ipccsd(nroots=3)
print(e[0] - 0.43356041409195489)
print(e[1] - 0.51876598058509493)
print(e[2] - 0.6782879569941862 )
e, v = mycc.eeccsd(nroots=4)
print(e[0] - 0.2757159395886167)
print(e[1] - 0.2757159395886167)
print(e[2] - 0.2757159395886167)
print(e[3] - 0.3005716731825082)
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 = '631g'
mol.build()
mf = scf.RHF(mol)
mf.conv_tol = 1e-16
mf.scf()
mo_coeff = mf.mo_coeff + np.sin(mf.mo_coeff) * .01j
nao = mo_coeff.shape[0]
eri = ao2mo.restore(1, mf._eri, nao)
eri0 = lib.einsum('pqrs,pi,qj,rk,sl->ijkl', eri, mo_coeff.conj(), mo_coeff,
mo_coeff.conj(), mo_coeff)
nocc, nvir = 5, nao-5
eris = _ChemistsERIs(mol)
eris.oooo = eri0[:nocc,:nocc,:nocc,:nocc].copy()
eris.ovoo = eri0[:nocc,nocc:,:nocc,:nocc].copy()
eris.oovv = eri0[:nocc,:nocc,nocc:,nocc:].copy()
eris.ovvo = eri0[:nocc,nocc:,nocc:,:nocc].copy()
eris.ovov = eri0[:nocc,nocc:,:nocc,nocc:].copy()
eris.ovvv = eri0[:nocc,nocc:,nocc:,nocc:].copy()
eris.vvvv = eri0[nocc:,nocc:,nocc:,nocc:].copy()
eris.fock = np.diag(mf.mo_energy)
eris.mo_energy = mf.mo_energy
np.random.seed(1)
t1 = np.random.random((nocc,nvir)) + np.random.random((nocc,nvir))*.1j - .5
t2 = np.random.random((nocc,nocc,nvir,nvir)) - .5
t2 = t2 + np.sin(t2) * .1j
t2 = t2 + t2.transpose(1,0,3,2)
mycc = RCCSD(mf)
t1new_ref, t2new_ref = update_amps(mycc, t1, t2, eris)
orbspin = np.zeros(nao*2, dtype=int)
orbspin[1::2] = 1
eri1 = np.zeros([nao*2]*4, dtype=np.complex128)
eri1[0::2,0::2,0::2,0::2] = eri1[0::2,0::2,1::2,1::2] = eri0
eri1[1::2,1::2,0::2,0::2] = eri1[1::2,1::2,1::2,1::2] = eri0
eri1 = eri1.transpose(0,2,1,3) - eri1.transpose(0,2,3,1)
erig = gccsd._PhysicistsERIs(mol)
nocc *= 2
nvir *= 2
erig.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy()
erig.ooov = eri1[:nocc,:nocc,:nocc,nocc:].copy()
erig.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy()
erig.ovvo = eri1[:nocc,nocc:,nocc:,:nocc].copy()
erig.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy()
erig.ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy()
erig.vvvv = eri1[nocc:,nocc:,nocc:,nocc:].copy()
mo_e = np.array([mf.mo_energy]*2)
erig.fock = np.diag(mo_e.T.ravel())
myccg = gccsd.GCCSD(scf.addons.convert_to_ghf(mf))
t1, t2 = myccg.amplitudes_from_ccsd(t1, t2)
t1new, t2new = gccsd.update_amps(myccg, t1, t2, erig)
print(abs(t1new[0::2,0::2]-t1new_ref).max())
t2aa = t2new[0::2,0::2,0::2,0::2]
t2ab = t2new[0::2,1::2,0::2,1::2]
print(abs(t2ab-t2new_ref).max())
print(abs(t2ab-t2ab.transpose(1,0,2,3) - t2aa).max())
