1#!/usr/bin/env python 2# Copyright 2014-2021 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, max_memory=2000, 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 35 nocc, nvir = t1.shape 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 40 mo_e = eris.mo_energy 41 eia = mo_e[:nocc,None] - mo_e[nocc:] 42 d3 = lib.direct_sum('ia+jb+kc->ijkabc', eia, eia, eia) 43 44 t3c =(numpy.einsum('jkae,bcei->ijkabc', t2, bcei) - 45 numpy.einsum('imbc,majk->ijkabc', t2, majk)) 46 t3c = t3c - t3c.transpose(0,1,2,4,3,5) - t3c.transpose(0,1,2,5,4,3) 47 t3c = t3c - t3c.transpose(1,0,2,3,4,5) - t3c.transpose(2,1,0,3,4,5) 48 t3c /= d3 49 # e4 = numpy.einsum('ijkabc,ijkabc,ijkabc', t3c.conj(), d3, t3c) / 36 50 # sia = numpy.einsum('jkbc,ijkabc->ia', eris.oovv, t3c) * .25 51 # e5 = numpy.einsum('ia,ia', sia, t1.conj()) 52 # et = e4 + e5 53 # return et 54 t3d = numpy.einsum('ia,bcjk->ijkabc', t1, bcjk) 55 t3d += numpy.einsum('ai,jkbc->ijkabc', eris.fock[nocc:,:nocc], t2) 56 t3d = t3d - t3d.transpose(0,1,2,4,3,5) - t3d.transpose(0,1,2,5,4,3) 57 t3d = t3d - t3d.transpose(1,0,2,3,4,5) - t3d.transpose(2,1,0,3,4,5) 58 t3d /= d3 59 et = numpy.einsum('ijkabc,ijkabc,ijkabc', (t3c+t3d).conj(), d3, t3c) / 36 60 return et 61 62 63if __name__ == '__main__': 64 from pyscf import gto 65 from pyscf import scf 66 from pyscf import cc 67 68 mol = gto.Mole() 69 mol.atom = [ 70 [8 , (0. , 0. , 0.)], 71 [1 , (0. , -.957 , .587)], 72 [1 , (0.2, .757 , .487)]] 73 74 mol.basis = '631g' 75 mol.build() 76 mf = scf.RHF(mol).run(conv_tol=1e-1) 77 mycc = cc.CCSD(mf).set(conv_tol=1e-11).run() 78 et = mycc.ccsd_t() 79 80 mycc = cc.GCCSD(scf.addons.convert_to_ghf(mf)).set(conv_tol=1e-11).run() 81 eris = mycc.ao2mo() 82 print(kernel(mycc, eris) - et) 83