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