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