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: Timothy Berkelbach <tim.berkelbach@gmail.com>
17#
18
19'''
20Restricted QCISD implementation
21The 4-index integrals are saved on disk entirely (without using any symmetry).
22
23Note MO integrals are treated in chemist's notation
24
25Ref:
26'''
27
28
29import numpy
30import numpy as np
31
32from pyscf import lib
33from pyscf.lib import logger
34from pyscf.cc import rccsd_slow as rccsd
35from pyscf.cc import rintermediates as imd
36from pyscf import __config__
37
38BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4)
39MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000)
40
41
42def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
43           tolnormt=1e-6, verbose=None):
44    '''Same as ccsd.kernel with strings modified to correct the method name'''
45    log = logger.new_logger(mycc, verbose)
46    if eris is None:
47        eris = mycc.ao2mo(mycc.mo_coeff)
48    if t1 is None and t2 is None:
49        t1, t2 = mycc.get_init_guess(eris)
50    elif t2 is None:
51        t2 = mycc.get_init_guess(eris)[1]
52
53    cput1 = cput0 = (logger.process_clock(), logger.perf_counter())
54    eold = 0
55    eccsd = mycc.energy(t1, t2, eris)
56    log.info('Init E_corr(QCISD) = %.15g', eccsd)
57
58    if isinstance(mycc.diis, lib.diis.DIIS):
59        adiis = mycc.diis
60    elif mycc.diis:
61        adiis = lib.diis.DIIS(mycc, mycc.diis_file, incore=mycc.incore_complete)
62        adiis.space = mycc.diis_space
63    else:
64        adiis = None
65
66    conv = False
67    for istep in range(max_cycle):
68        t1new, t2new = mycc.update_amps(t1, t2, eris)
69        tmpvec = mycc.amplitudes_to_vector(t1new, t2new)
70        tmpvec -= mycc.amplitudes_to_vector(t1, t2)
71        normt = numpy.linalg.norm(tmpvec)
72        tmpvec = None
73        if mycc.iterative_damping < 1.0:
74            alpha = mycc.iterative_damping
75            t1new = (1-alpha) * t1 + alpha * t1new
76            t2new *= alpha
77            t2new += (1-alpha) * t2
78        t1, t2 = t1new, t2new
79        t1new = t2new = None
80        t1, t2 = mycc.run_diis(t1, t2, istep, normt, eccsd-eold, adiis)
81        eold, eccsd = eccsd, mycc.energy(t1, t2, eris)
82        log.info('cycle = %d  E_corr(QCISD) = %.15g  dE = %.9g  norm(t1,t2) = %.6g',
83                 istep+1, eccsd, eccsd - eold, normt)
84        cput1 = log.timer('QCISD iter', *cput1)
85        if abs(eccsd-eold) < tol and normt < tolnormt:
86            conv = True
87            break
88    log.timer('QCISD', *cput0)
89    return conv, eccsd, t1, t2
90
91
92def update_amps(cc, t1, t2, eris):
93    # Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004) Eqs.(35)-(36)
94    nocc, nvir = t1.shape
95    fock = eris.fock
96
97    fov = fock[:nocc,nocc:].copy()
98    foo = fock[:nocc,:nocc].copy()
99    fvv = fock[nocc:,nocc:].copy()
100
101    Foo = imd.cc_Foo(0*t1,t2,eris)
102    Fvv = imd.cc_Fvv(0*t1,t2,eris)
103    Fov = imd.cc_Fov(t1,t2,eris)
104
105    Foo -= np.diag(np.diag(foo))
106    Fvv -= np.diag(np.diag(fvv))
107
108    # T1 equation
109    t1new = np.asarray(fov).conj().copy()
110    t1new +=   lib.einsum('ac,ic->ia', Fvv, t1)
111    t1new +=  -lib.einsum('ki,ka->ia', Foo, t1)
112    t1new += 2*lib.einsum('kc,kica->ia', Fov, t2)
113    t1new +=  -lib.einsum('kc,ikca->ia', Fov, t2)
114    t1new += 2*lib.einsum('kcai,kc->ia', eris.ovvo, t1)
115    t1new +=  -lib.einsum('kiac,kc->ia', eris.oovv, t1)
116    eris_ovvv = np.asarray(eris.ovvv)
117    t1new += 2*lib.einsum('kdac,ikcd->ia', eris_ovvv, t2)
118    t1new +=  -lib.einsum('kcad,ikcd->ia', eris_ovvv, t2)
119    t1new +=-2*lib.einsum('kilc,klac->ia', eris.ooov, t2)
120    t1new +=   lib.einsum('likc,klac->ia', eris.ooov, t2)
121
122    # T2 equation
123    t2new = np.asarray(eris.ovov).conj().transpose(0,2,1,3).copy()
124    Loo = imd.Loo(0*t1, t2, eris)
125    Lvv = imd.Lvv(0*t1, t2, eris)
126    Loo -= np.diag(np.diag(foo))
127    Lvv -= np.diag(np.diag(fvv))
128    Woooo = imd.cc_Woooo(0*t1, t2, eris)
129    Wvoov = imd.cc_Wvoov(0*t1, t2, eris)
130    Wvovo = imd.cc_Wvovo(0*t1, t2, eris)
131    Wvvvv = imd.cc_Wvvvv(0*t1, t2, eris)
132    t2new += lib.einsum('klij,klab->ijab', Woooo, t2)
133    t2new += lib.einsum('abcd,ijcd->ijab', Wvvvv, t2)
134    tmp = lib.einsum('ac,ijcb->ijab', Lvv, t2)
135    t2new += (tmp + tmp.transpose(1,0,3,2))
136    tmp = lib.einsum('ki,kjab->ijab', Loo, t2)
137    t2new -= (tmp + tmp.transpose(1,0,3,2))
138    tmp  = 2*lib.einsum('akic,kjcb->ijab', Wvoov, t2)
139    tmp -=   lib.einsum('akci,kjcb->ijab', Wvovo, t2)
140    t2new += (tmp + tmp.transpose(1,0,3,2))
141    tmp = lib.einsum('akic,kjbc->ijab', Wvoov, t2)
142    t2new -= (tmp + tmp.transpose(1,0,3,2))
143    tmp = lib.einsum('bkci,kjac->ijab', Wvovo, t2)
144    t2new -= (tmp + tmp.transpose(1,0,3,2))
145
146    tmp2 = np.asarray(eris.ovvv).conj().transpose(1,3,0,2)
147    tmp = lib.einsum('abic,jc->ijab', tmp2, t1)
148    t2new += (tmp + tmp.transpose(1,0,3,2))
149    tmp2 = np.asarray(eris.ooov).transpose(3,1,2,0).conj()
150    tmp = lib.einsum('akij,kb->ijab', tmp2, t1)
151    t2new -= (tmp + tmp.transpose(1,0,3,2))
152
153    mo_e = eris.fock.diagonal().real
154    eia = mo_e[:nocc,None] - mo_e[None,nocc:]
155    eijab = lib.direct_sum('ia,jb->ijab',eia,eia)
156    t1new /= eia
157    t2new /= eijab
158
159    return t1new, t2new
160
161
162class QCISD(rccsd.RCCSD):
163    '''restricted QCISD
164    '''
165
166    def kernel(self, t1=None, t2=None, eris=None):
167        return self.qcisd(t1, t2, eris)
168    def qcisd(self, t1=None, t2=None, eris=None):
169        assert(self.mo_coeff is not None)
170        assert(self.mo_occ is not None)
171
172        if self.verbose >= logger.WARN:
173            self.check_sanity()
174        self.dump_flags()
175
176        if eris is None:
177            eris = self.ao2mo(self.mo_coeff)
178
179        self.e_hf = getattr(eris, 'e_hf', None)
180        if self.e_hf is None:
181            self.e_hf = self._scf.e_tot
182
183        self.converged, self.e_corr, self.t1, self.t2 = \
184                kernel(self, eris, t1, t2, max_cycle=self.max_cycle,
185                       tol=self.conv_tol, tolnormt=self.conv_tol_normt,
186                       verbose=self.verbose)
187        self._finalize()
188        return self.e_corr, self.t1, self.t2
189
190    def energy(self, t1=None, t2=None, eris=None):
191        return rccsd.energy(self, t1*0, t2, eris)
192
193    update_amps = update_amps
194
195    def qcisd_t(self, t1=None, t2=None, eris=None):
196        from pyscf.cc import qcisd_t_slow as qcisd_t
197        if t1 is None: t1 = self.t1
198        if t2 is None: t2 = self.t2
199        if eris is None: eris = self.ao2mo(self.mo_coeff)
200        return qcisd_t.kernel(self, eris, t1, t2, self.verbose)
201
202    def density_fit(self, auxbasis=None, with_df=None):
203        raise NotImplementedError
204
205
206if __name__ == '__main__':
207    from pyscf import gto, scf
208
209    mol = gto.Mole()
210    mol.atom = """C  0.000  0.000  0.000
211                  H  0.637  0.637  0.637
212                  H -0.637 -0.637  0.637
213                  H -0.637  0.637 -0.637
214                  H  0.637 -0.637 -0.637"""
215    mol.basis = 'cc-pvdz'
216    mol.verbose = 7
217    mol.spin = 0
218    mol.build()
219    mf = scf.RHF(mol).run(conv_tol=1e-14)
220
221    mycc = QCISD(mf, frozen=1)
222    ecc, t1, t2 = mycc.kernel()
223    print(mycc.e_tot - -40.383989)
224    et = mycc.qcisd_t()
225    print(mycc.e_tot+et - -40.387679)
226