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