1#!/usr/bin/env python 2# Copyright 2014-2020 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18 19''' 20GHF-CCSD(T) with spin-orbital integrals 21''' 22 23import numpy 24from pyscf import lib 25from pyscf.lib import logger 26from pyscf.cc import gccsd 27 28# spin-orbital formula 29# JCP 98, 8718 (1993); DOI:10.1063/1.464480 30def kernel(cc, eris, t1=None, t2=None, verbose=logger.INFO): 31 assert(isinstance(eris, gccsd._PhysicistsERIs)) 32 if t1 is None or t2 is None: 33 t1, t2 = cc.t1, cc.t2 34 nocc, nvir = t1.shape 35 36 bcei = numpy.asarray(eris.ovvv).conj().transpose(3,2,1,0) 37 majk = numpy.asarray(eris.ooov).conj().transpose(2,3,0,1) 38 bcjk = numpy.asarray(eris.oovv).conj().transpose(2,3,0,1) 39 fvo = eris.fock[nocc:,:nocc] 40 mo_e = eris.mo_energy 41 eijk = lib.direct_sum('i+j+k->ijk', mo_e[:nocc], mo_e[:nocc], mo_e[:nocc]) 42 eabc = lib.direct_sum('a+b+c->abc', mo_e[nocc:], mo_e[nocc:], mo_e[nocc:]) 43 44 t2T = t2.transpose(2,3,0,1) 45 t1T = t1.T 46 def get_wv(a, b, c): 47 w = numpy.einsum('ejk,ei->ijk', t2T[a,:], bcei[b,c]) 48 w -= numpy.einsum('im,mjk->ijk', t2T[b,c], majk[:,a]) 49 v = numpy.einsum('i,jk->ijk', t1T[a], bcjk[b,c]) 50 v += numpy.einsum('i,jk->ijk', fvo[a], t2T [b,c]) 51 v += w 52 w = w + w.transpose(2,0,1) + w.transpose(1,2,0) 53 return w, v 54 55 et = 0 56 for a in range(nvir): 57 for b in range(a): 58 for c in range(b): 59 wabc, vabc = get_wv(a, b, c) 60 wcab, vcab = get_wv(c, a, b) 61 wbac, vbac = get_wv(b, a, c) 62 63 w = wabc + wcab - wbac 64 v = vabc + vcab - vbac 65 w /= eijk - eabc[a,b,c] 66 et += numpy.einsum('ijk,ijk', w, v.conj()) 67 et /= 2 68 return et 69 70 71if __name__ == '__main__': 72 from pyscf import gto 73 from pyscf import scf 74 from pyscf import cc 75 76 mol = gto.Mole() 77 mol.atom = [ 78 [8 , (0. , 0. , 0.)], 79 [1 , (0. , -.957 , .587)], 80 [1 , (0.2, .757 , .487)]] 81 82 mol.basis = '631g' 83 mol.build() 84 mf = scf.RHF(mol).run(conv_tol=1e-1) 85 mycc = cc.CCSD(mf).set(conv_tol=1e-11).run() 86 et = mycc.ccsd_t() 87 88 mycc = cc.GCCSD(scf.addons.convert_to_ghf(mf)).set(conv_tol=1e-11).run() 89 eris = mycc.ao2mo() 90 print(kernel(mycc, eris) - et) 91 92 numpy.random.seed(1) 93 mf.mo_coeff = numpy.random.random(mf.mo_coeff.shape) - .9 94 mycc = cc.GCCSD(scf.addons.convert_to_ghf(mf)) 95 eris = mycc.ao2mo() 96 nocc = 10 97 nvir = mol.nao_nr() * 2 - nocc 98 t1 = numpy.random.random((nocc,nvir))*.1 - .2 99 t2 = numpy.random.random((nocc,nocc,nvir,nvir))*.1 - .2 100 t2 = t2 - t2.transpose(1,0,2,3) 101 t2 = t2 - t2.transpose(0,1,3,2) 102 print(kernel(mycc, eris, t1, t2) - 263713.3945021223) 103