1# Copyright 2014-2019 The PySCF Developers. All Rights Reserved. 2# 3# Licensed under the Apache License, Version 2.0 (the "License"); 4# you may not use this file except in compliance with the License. 5# You may obtain a copy of the License at 6# 7# http://www.apache.org/licenses/LICENSE-2.0 8# 9# Unless required by applicable law or agreed to in writing, software 10# distributed under the License is distributed on an "AS IS" BASIS, 11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12# See the License for the specific language governing permissions and 13# limitations under the License. 14# 15# Author: Samragni Banerjee <samragnibanerjee4@gmail.com> 16# Alexander Sokolov <alexander.y.sokolov@gmail.com> 17# 18 19''' 20Restricted algebraic diagrammatic construction 21''' 22import numpy as np 23import pyscf.ao2mo as ao2mo 24from pyscf import lib 25from pyscf.lib import logger 26from pyscf.adc import radc_ao2mo 27from pyscf.adc import dfadc 28from pyscf import __config__ 29from pyscf import df 30from pyscf import symm 31 32def kernel(adc, nroots=1, guess=None, eris=None, verbose=None): 33 34 adc.method = adc.method.lower() 35 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 36 raise NotImplementedError(adc.method) 37 38 cput0 = (logger.process_clock(), logger.perf_counter()) 39 log = logger.Logger(adc.stdout, adc.verbose) 40 if adc.verbose >= logger.WARN: 41 adc.check_sanity() 42 adc.dump_flags() 43 44 if eris is None: 45 eris = adc.transform_integrals() 46 47 imds = adc.get_imds(eris) 48 matvec, diag = adc.gen_matvec(imds, eris) 49 50 guess = adc.get_init_guess(nroots, diag, ascending = True) 51 52 conv, adc.E, U = lib.linalg_helper.davidson_nosym1( 53 lambda xs : [matvec(x) for x in xs], 54 guess, diag, nroots=nroots, verbose=log, tol=adc.conv_tol, 55 max_cycle=adc.max_cycle, max_space=adc.max_space, tol_residual=adc.tol_residual) 56 57 adc.U = np.array(U).T.copy() 58 59 if adc.compute_properties: 60 adc.P,adc.X = adc.get_properties(nroots) 61 62 nfalse = np.shape(conv)[0] - np.sum(conv) 63 64 header = ("\n*************************************************************" 65 "\n ADC calculation summary" 66 "\n*************************************************************") 67 logger.info(adc, header) 68 69 if nfalse >= 1: 70 logger.warn(adc, "Davidson iterations for " + str(nfalse) + " root(s) not converged\n") 71 72 for n in range(nroots): 73 print_string = ('%s root %d | Energy (Eh) = %14.10f | Energy (eV) = %12.8f ' % 74 (adc.method, n, adc.E[n], adc.E[n]*27.2114)) 75 if adc.compute_properties: 76 print_string += ("| Spec factors = %10.8f " % adc.P[n]) 77 print_string += ("| conv = %s" % conv[n]) 78 logger.info(adc, print_string) 79 80 log.timer('ADC', *cput0) 81 82 return adc.E, adc.U, adc.P, adc.X 83 84 85def compute_amplitudes_energy(myadc, eris, verbose=None): 86 87 t1, t2, myadc.imds.t2_1_vvvv = myadc.compute_amplitudes(eris) 88 e_corr = myadc.compute_energy(t2, eris) 89 90 return e_corr, t1, t2 91 92 93def compute_amplitudes(myadc, eris): 94 95 cput0 = (logger.process_clock(), logger.perf_counter()) 96 log = logger.Logger(myadc.stdout, myadc.verbose) 97 98 if myadc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 99 raise NotImplementedError(myadc.method) 100 101 nocc = myadc._nocc 102 nvir = myadc._nvir 103 104 eris_oooo = eris.oooo 105 eris_ovoo = eris.ovoo 106 eris_oovv = eris.oovv 107 eris_ovvo = eris.ovvo 108 109 # Compute first-order doubles t2 (tijab) 110 111 v2e_oovv = eris_ovvo[:].transpose(0,3,1,2).copy() 112 113 e = myadc.mo_energy 114 d_ij = e[:nocc][:,None] + e[:nocc] 115 d_ab = e[nocc:][:,None] + e[nocc:] 116 117 D2 = d_ij.reshape(-1,1) - d_ab.reshape(-1) 118 D2 = D2.reshape((nocc,nocc,nvir,nvir)) 119 120 D1 = e[:nocc][:None].reshape(-1,1) - e[nocc:].reshape(-1) 121 D1 = D1.reshape((nocc,nvir)) 122 123 t2_1 = v2e_oovv/D2 124 if not isinstance(eris.oooo, np.ndarray): 125 t2_1 = radc_ao2mo.write_dataset(t2_1) 126 127 del v2e_oovv 128 del D2 129 130 cput0 = log.timer_debug1("Completed t2_1 amplitude calculation", *cput0) 131 132 # Compute second-order singles t1 (tij) 133 134 if isinstance(eris.ovvv, type(None)): 135 chnk_size = radc_ao2mo.calculate_chunk_size(myadc) 136 else: 137 chnk_size = nocc 138 a = 0 139 t1_2 = np.zeros((nocc,nvir)) 140 141 for p in range(0,nocc,chnk_size): 142 if getattr(myadc, 'with_df', None): 143 eris_ovvv = dfadc.get_ovvv_df(myadc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 144 else: 145 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 146 k = eris_ovvv.shape[0] 147 148 t1_2 += 0.5*lib.einsum('kdac,ikcd->ia',eris_ovvv,t2_1[:,a:a+k],optimize=True) 149 t1_2 -= 0.5*lib.einsum('kdac,kicd->ia',eris_ovvv,t2_1[a:a+k,:],optimize=True) 150 t1_2 -= 0.5*lib.einsum('kcad,ikcd->ia',eris_ovvv,t2_1[:,a:a+k],optimize=True) 151 t1_2 += 0.5*lib.einsum('kcad,kicd->ia',eris_ovvv,t2_1[a:a+k,:],optimize=True) 152 153 t1_2 += lib.einsum('kdac,ikcd->ia',eris_ovvv,t2_1[:,a:a+k],optimize=True) 154 del eris_ovvv 155 a += k 156 157 t1_2 -= 0.5*lib.einsum('lcki,klac->ia',eris_ovoo,t2_1[:],optimize=True) 158 t1_2 += 0.5*lib.einsum('lcki,lkac->ia',eris_ovoo,t2_1[:],optimize=True) 159 t1_2 -= 0.5*lib.einsum('kcli,lkac->ia',eris_ovoo,t2_1[:],optimize=True) 160 t1_2 += 0.5*lib.einsum('kcli,klac->ia',eris_ovoo,t2_1[:],optimize=True) 161 t1_2 -= lib.einsum('lcki,klac->ia',eris_ovoo,t2_1[:],optimize=True) 162 163 t1_2 = t1_2/D1 164 165 cput0 = log.timer_debug1("Completed t1_2 amplitude calculation", *cput0) 166 167 t2_2 = None 168 t1_3 = None 169 t2_1_vvvv = None 170 171 if (myadc.method == "adc(2)-x" or myadc.method == "adc(3)"): 172 173 # Compute second-order doubles t2 (tijab) 174 175 eris_oooo = eris.oooo 176 eris_ovvo = eris.ovvo 177 178 if isinstance(eris.vvvv, np.ndarray): 179 eris_vvvv = eris.vvvv 180 temp = t2_1.reshape(nocc*nocc,nvir*nvir) 181 t2_1_vvvv = np.dot(temp,eris_vvvv.T).reshape(nocc,nocc,nvir,nvir) 182 elif isinstance(eris.vvvv, list): 183 t2_1_vvvv = contract_ladder(myadc,t2_1[:],eris.vvvv) 184 else: 185 t2_1_vvvv = contract_ladder(myadc,t2_1[:],eris.Lvv) 186 187 if not isinstance(eris.oooo, np.ndarray): 188 t2_1_vvvv = radc_ao2mo.write_dataset(t2_1_vvvv) 189 190 t2_2 = t2_1_vvvv[:].copy() 191 192 t2_2 += lib.einsum('kilj,klab->ijab',eris_oooo,t2_1[:],optimize=True) 193 t2_2 += lib.einsum('kcbj,kica->ijab',eris_ovvo,t2_1[:],optimize=True) 194 t2_2 -= lib.einsum('kcbj,ikca->ijab',eris_ovvo,t2_1[:],optimize=True) 195 t2_2 += lib.einsum('kcbj,ikac->ijab',eris_ovvo,t2_1[:],optimize=True) 196 t2_2 -= lib.einsum('kjbc,ikac->ijab',eris_oovv,t2_1[:],optimize=True) 197 t2_2 -= lib.einsum('kibc,kjac->ijab',eris_oovv,t2_1[:],optimize=True) 198 t2_2 -= lib.einsum('kjac,ikcb->ijab',eris_oovv,t2_1[:],optimize=True) 199 t2_2 += lib.einsum('kcai,kjcb->ijab',eris_ovvo,t2_1[:],optimize=True) 200 t2_2 -= lib.einsum('kcai,jkcb->ijab',eris_ovvo,t2_1[:],optimize=True) 201 t2_2 += lib.einsum('kcai,kjcb->ijab',eris_ovvo,t2_1[:],optimize=True) 202 t2_2 -= lib.einsum('kiac,kjcb->ijab',eris_oovv,t2_1[:],optimize=True) 203 204 D2 = d_ij.reshape(-1,1) - d_ab.reshape(-1) 205 D2 = D2.reshape((nocc,nocc,nvir,nvir)) 206 207 t2_2 = t2_2/D2 208 if not isinstance(eris.oooo, np.ndarray): 209 t2_2 = radc_ao2mo.write_dataset(t2_2) 210 del D2 211 212 cput0 = log.timer_debug1("Completed t2_2 amplitude calculation", *cput0) 213 214 if (myadc.method == "adc(3)"): 215 216 eris_ovoo = eris.ovoo 217 218 t1_3 = lib.einsum('d,ilad,ld->ia',e[nocc:],t2_1[:],t1_2,optimize=True) 219 t1_3 -= lib.einsum('d,liad,ld->ia',e[nocc:],t2_1[:],t1_2,optimize=True) 220 t1_3 += lib.einsum('d,ilad,ld->ia',e[nocc:],t2_1[:],t1_2,optimize=True) 221 222 t1_3 -= lib.einsum('l,ilad,ld->ia',e[:nocc],t2_1[:], t1_2,optimize=True) 223 t1_3 += lib.einsum('l,liad,ld->ia',e[:nocc],t2_1[:], t1_2,optimize=True) 224 t1_3 -= lib.einsum('l,ilad,ld->ia',e[:nocc],t2_1[:],t1_2,optimize=True) 225 226 t1_3 += 0.5*lib.einsum('a,ilad,ld->ia',e[nocc:],t2_1[:], t1_2,optimize=True) 227 t1_3 -= 0.5*lib.einsum('a,liad,ld->ia',e[nocc:],t2_1[:], t1_2,optimize=True) 228 t1_3 += 0.5*lib.einsum('a,ilad,ld->ia',e[nocc:],t2_1[:],t1_2,optimize=True) 229 230 t1_3 -= 0.5*lib.einsum('i,ilad,ld->ia',e[:nocc],t2_1[:], t1_2,optimize=True) 231 t1_3 += 0.5*lib.einsum('i,liad,ld->ia',e[:nocc],t2_1[:], t1_2,optimize=True) 232 t1_3 -= 0.5*lib.einsum('i,ilad,ld->ia',e[:nocc],t2_1[:],t1_2,optimize=True) 233 234 t1_3 += lib.einsum('ld,iadl->ia',t1_2,eris_ovvo,optimize=True) 235 t1_3 -= lib.einsum('ld,ladi->ia',t1_2,eris_ovvo,optimize=True) 236 t1_3 += lib.einsum('ld,iadl->ia',t1_2,eris_ovvo,optimize=True) 237 238 t1_3 += lib.einsum('ld,ldai->ia',t1_2,eris_ovvo ,optimize=True) 239 t1_3 -= lib.einsum('ld,liad->ia',t1_2,eris_oovv ,optimize=True) 240 t1_3 += lib.einsum('ld,ldai->ia',t1_2,eris_ovvo,optimize=True) 241 242 t1_3 -= 0.5*lib.einsum('lmad,mdli->ia',t2_2[:],eris_ovoo,optimize=True) 243 t1_3 += 0.5*lib.einsum('mlad,mdli->ia',t2_2[:],eris_ovoo,optimize=True) 244 t1_3 += 0.5*lib.einsum('lmad,ldmi->ia',t2_2[:],eris_ovoo,optimize=True) 245 t1_3 -= 0.5*lib.einsum('mlad,ldmi->ia',t2_2[:],eris_ovoo,optimize=True) 246 t1_3 -= lib.einsum('lmad,mdli->ia',t2_2[:],eris_ovoo,optimize=True) 247 248 if isinstance(eris.ovvv, type(None)): 249 chnk_size = radc_ao2mo.calculate_chunk_size(myadc) 250 else: 251 chnk_size = nocc 252 a = 0 253 254 for p in range(0,nocc,chnk_size): 255 if getattr(myadc, 'with_df', None): 256 eris_ovvv = dfadc.get_ovvv_df(myadc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 257 else: 258 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 259 k = eris_ovvv.shape[0] 260 261 t1_3 += 0.5*lib.einsum('ilde,lead->ia', t2_2[:,a:a+k],eris_ovvv,optimize=True) 262 t1_3 -= 0.5*lib.einsum('lide,lead->ia', t2_2[a:a+k],eris_ovvv,optimize=True) 263 264 t1_3 -= 0.5*lib.einsum('ilde,ldae->ia', t2_2[:,a:a+k],eris_ovvv,optimize=True) 265 t1_3 += 0.5*lib.einsum('lide,ldae->ia', t2_2[a:a+k],eris_ovvv,optimize=True) 266 267 t1_3 -= lib.einsum('ildf,mefa,lmde->ia',t2_1[:], eris_ovvv, t2_1[:,a:a+k] ,optimize=True) 268 t1_3 += lib.einsum('ildf,mefa,mlde->ia',t2_1[:], eris_ovvv, t2_1[a:a+k] ,optimize=True) 269 t1_3 += lib.einsum('lidf,mefa,lmde->ia',t2_1[:], eris_ovvv, t2_1[:,a:a+k] ,optimize=True) 270 t1_3 -= lib.einsum('lidf,mefa,mlde->ia',t2_1[:], eris_ovvv, t2_1[a:a+k] ,optimize=True) 271 272 t1_3 += lib.einsum('ildf,mafe,lmde->ia',t2_1[:], eris_ovvv, t2_1[:,a:a+k] ,optimize=True) 273 t1_3 -= lib.einsum('ildf,mafe,mlde->ia',t2_1[:], eris_ovvv, t2_1[a:a+k] ,optimize=True) 274 t1_3 -= lib.einsum('lidf,mafe,lmde->ia',t2_1[:], eris_ovvv, t2_1[:,a:a+k] ,optimize=True) 275 t1_3 += lib.einsum('lidf,mafe,mlde->ia',t2_1[:], eris_ovvv, t2_1[a:a+k] ,optimize=True) 276 277 t1_3 += lib.einsum('ilfd,mefa,mled->ia', t2_1[:],eris_ovvv, t2_1[a:a+k],optimize=True) 278 t1_3 -= lib.einsum('ilfd,mafe,mled->ia', t2_1[:],eris_ovvv, t2_1[a:a+k],optimize=True) 279 280 t1_3 += 0.5*lib.einsum('ilaf,mefd,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 281 t1_3 -= 0.5*lib.einsum('ilaf,mefd,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 282 t1_3 -= 0.5*lib.einsum('liaf,mefd,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 283 t1_3 += 0.5*lib.einsum('liaf,mefd,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 284 285 t1_3 -= 0.5*lib.einsum('ilaf,mdfe,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 286 t1_3 += 0.5*lib.einsum('ilaf,mdfe,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 287 t1_3 += 0.5*lib.einsum('liaf,mdfe,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 288 t1_3 -= 0.5*lib.einsum('liaf,mdfe,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 289 290 t1_3[a:a+k] += 0.5*lib.einsum('lmdf,iaef,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 291 t1_3[a:a+k] -= 0.5*lib.einsum('lmdf,iaef,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 292 t1_3[a:a+k] -= 0.5*lib.einsum('mldf,iaef,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 293 t1_3[a:a+k] += 0.5*lib.einsum('mldf,iaef,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 294 295 t1_3[a:a+k] -= 0.5*lib.einsum('lmdf,ifea,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 296 t1_3[a:a+k] += 0.5*lib.einsum('lmdf,ifea,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 297 t1_3[a:a+k] += 0.5*lib.einsum('mldf,ifea,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 298 t1_3[a:a+k] -= 0.5*lib.einsum('mldf,ifea,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 299 300 t1_3[a:a+k] += lib.einsum('mlfd,iaef,mled->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 301 t1_3[a:a+k] -= lib.einsum('mlfd,ifea,mled->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 302 303 t1_3[a:a+k] -= 0.25*lib.einsum('lmef,iedf,lmad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 304 t1_3[a:a+k] += 0.25*lib.einsum('lmef,iedf,mlad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 305 t1_3[a:a+k] += 0.25*lib.einsum('mlef,iedf,lmad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 306 t1_3[a:a+k] -= 0.25*lib.einsum('mlef,iedf,mlad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 307 308 t1_3[a:a+k] += 0.25*lib.einsum('lmef,ifde,lmad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 309 t1_3[a:a+k] -= 0.25*lib.einsum('lmef,ifde,mlad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 310 t1_3[a:a+k] -= 0.25*lib.einsum('mlef,ifde,lmad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 311 t1_3[a:a+k] += 0.25*lib.einsum('mlef,ifde,mlad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 312 313 t1_3 += 0.5*lib.einsum('ilaf,mefd,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 314 t1_3 -= 0.5*lib.einsum('ilaf,mefd,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 315 316 t1_3 -= 0.5*lib.einsum('ilaf,mdfe,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 317 t1_3 += 0.5*lib.einsum('ilaf,mdfe,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 318 319 t1_3 -= lib.einsum('ildf,mafe,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 320 t1_3 += lib.einsum('ilaf,mefd,mled->ia',t2_1[:],eris_ovvv,t2_1[a:a+k],optimize=True) 321 322 t1_3[a:a+k] += 0.5*lib.einsum('lmdf,iaef,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 323 t1_3[a:a+k] -= 0.5*lib.einsum('lmdf,iaef,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 324 t1_3[a:a+k] -= 0.5*lib.einsum('mldf,iaef,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 325 t1_3[a:a+k] += 0.5*lib.einsum('mldf,iaef,mlde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 326 327 t1_3[a:a+k] += lib.einsum('lmdf,iaef,lmde->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 328 t1_3[a:a+k] -= lib.einsum('lmef,iedf,lmad->ia',t2_1[:],eris_ovvv,t2_1[:],optimize=True) 329 330 t1_3 += lib.einsum('ilde,lead->ia',t2_2[:,a:a+k],eris_ovvv,optimize=True) 331 332 t1_3 -= lib.einsum('ildf,mefa,lmde->ia',t2_1[:],eris_ovvv, t2_1[:,a:a+k],optimize=True) 333 t1_3 += lib.einsum('lidf,mefa,lmde->ia',t2_1[:],eris_ovvv, t2_1[:,a:a+k],optimize=True) 334 335 t1_3 += lib.einsum('ilfd,mefa,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k] ,optimize=True) 336 t1_3 -= lib.einsum('ilfd,mefa,mlde->ia',t2_1[:],eris_ovvv,t2_1[a:a+k] ,optimize=True) 337 338 t1_3 += lib.einsum('ilaf,mefd,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 339 t1_3 -= lib.einsum('liaf,mefd,lmde->ia',t2_1[:],eris_ovvv,t2_1[:,a:a+k],optimize=True) 340 341 del eris_ovvv 342 a += k 343 344 t1_3 += 0.25*lib.einsum('inde,lamn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 345 t1_3 -= 0.25*lib.einsum('inde,lamn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 346 t1_3 -= 0.25*lib.einsum('nide,lamn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 347 t1_3 += 0.25*lib.einsum('nide,lamn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 348 t1_3 -= 0.25*lib.einsum('inde,maln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 349 t1_3 += 0.25*lib.einsum('inde,maln,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 350 t1_3 += 0.25*lib.einsum('nide,maln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 351 t1_3 -= 0.25*lib.einsum('nide,maln,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 352 t1_3 += lib.einsum('inde,lamn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 353 354 t1_3 += 0.5 * lib.einsum('inad,lemn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 355 t1_3 -= 0.5 * lib.einsum('inad,lemn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 356 t1_3 -= 0.5 * lib.einsum('niad,lemn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 357 t1_3 += 0.5 * lib.einsum('niad,lemn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 358 359 t1_3 -= 0.5 * lib.einsum('inad,meln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 360 t1_3 += 0.5 * lib.einsum('inad,meln,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 361 t1_3 += 0.5 * lib.einsum('niad,meln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 362 t1_3 -= 0.5 * lib.einsum('niad,meln,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 363 364 t1_3 -= 0.5 * lib.einsum('inad,lemn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 365 t1_3 += 0.5 * lib.einsum('niad,lemn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 366 367 t1_3 -= 0.5 * lib.einsum('inad,meln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 368 t1_3 += 0.5 * lib.einsum('niad,meln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 369 370 t1_3 -= 0.5 * lib.einsum('inad,lemn,lmed->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 371 t1_3 -= 0.5 * lib.einsum('inad,meln,mled->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 372 373 t1_3 += 0.5 * lib.einsum('inad,lemn,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 374 t1_3 -= 0.5 * lib.einsum('inad,lemn,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 375 376 t1_3 -= 0.5 * lib.einsum('inad,meln,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 377 t1_3 += 0.5 * lib.einsum('inad,meln,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 378 379 t1_3 -= 0.5 * lib.einsum('lnde,ianm,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 380 t1_3 += 0.5 * lib.einsum('lnde,ianm,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 381 t1_3 += 0.5 * lib.einsum('nlde,ianm,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 382 t1_3 -= 0.5 * lib.einsum('nlde,ianm,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 383 384 t1_3 += 0.5 * lib.einsum('lnde,naim,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 385 t1_3 -= 0.5 * lib.einsum('lnde,naim,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 386 t1_3 -= 0.5 * lib.einsum('nlde,naim,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 387 t1_3 += 0.5 * lib.einsum('nlde,naim,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 388 389 t1_3 -= lib.einsum('nled,ianm,mled->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 390 t1_3 += lib.einsum('nled,naim,mled->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 391 392 t1_3 -= 0.5*lib.einsum('lnde,ianm,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 393 t1_3 += 0.5*lib.einsum('lnde,ianm,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 394 t1_3 += 0.5*lib.einsum('nlde,ianm,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 395 t1_3 -= 0.5*lib.einsum('nlde,ianm,mlde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 396 397 t1_3 -= lib.einsum('lnde,ianm,lmde->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 398 399 t1_3 -= lib.einsum('lnde,ienm,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 400 t1_3 += lib.einsum('lnde,ienm,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 401 t1_3 += lib.einsum('nlde,ienm,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 402 t1_3 -= lib.einsum('nlde,ienm,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 403 404 t1_3 += lib.einsum('lnde,neim,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 405 t1_3 -= lib.einsum('lnde,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 406 t1_3 -= lib.einsum('nlde,neim,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 407 t1_3 += lib.einsum('nlde,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 408 409 t1_3 += lib.einsum('lnde,neim,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 410 t1_3 -= lib.einsum('lnde,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 411 412 t1_3 += lib.einsum('nled,ienm,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 413 t1_3 -= lib.einsum('nled,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 414 t1_3 += lib.einsum('lned,ienm,lmad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 415 416 t1_3 -= lib.einsum('lnde,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 417 t1_3 += lib.einsum('nlde,neim,mlad->ia',t2_1[:],eris_ovoo,t2_1[:],optimize=True) 418 419 t1_3 = t1_3/D1 420 421 cput0 = log.timer_debug1("Completed amplitude calculation", *cput0) 422 423 t1 = (t1_2, t1_3) 424 t2 = (t2_1, t2_2) 425 426 return t1, t2, t2_1_vvvv 427 428 429def compute_energy(myadc, t2, eris): 430 431 if myadc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 432 raise NotImplementedError(myadc.method) 433 434 eris_ovvo = eris.ovvo 435 436 t2_new = t2[0][:].copy() 437 438 if (myadc.method == "adc(3)"): 439 t2_new += t2[1][:] 440 441 #Compute MP2 correlation energy 442 443 e_mp = 0.5 * lib.einsum('ijab,iabj', t2_new, eris_ovvo,optimize=True) 444 e_mp -= 0.5 * lib.einsum('ijab,ibaj', t2_new, eris_ovvo,optimize=True) 445 e_mp -= 0.5 * lib.einsum('jiab,iabj', t2_new, eris_ovvo,optimize=True) 446 e_mp += 0.5 * lib.einsum('jiab,ibaj', t2_new, eris_ovvo,optimize=True) 447 e_mp += lib.einsum('ijab,iabj', t2_new, eris_ovvo,optimize=True) 448 449 del t2_new 450 return e_mp 451 452 453def contract_ladder(myadc,t_amp,vvvv): 454 455 nocc = myadc._nocc 456 nvir = myadc._nvir 457 458 t_amp = np.ascontiguousarray(t_amp.reshape(nocc*nocc,nvir*nvir).T) 459 t = np.zeros((nvir,nvir, nocc*nocc)) 460 chnk_size = radc_ao2mo.calculate_chunk_size(myadc) 461 462 a = 0 463 if isinstance(vvvv, list): 464 for dataset in vvvv: 465 k = dataset.shape[0] 466 dataset = dataset[:].reshape(-1,nvir*nvir) 467 t[a:a+k] = np.dot(dataset,t_amp).reshape(-1,nvir,nocc*nocc) 468 a += k 469 elif getattr(myadc, 'with_df', None): 470 for p in range(0,nvir,chnk_size): 471 vvvv_p = dfadc.get_vvvv_df(myadc, vvvv, p, chnk_size) 472 k = vvvv_p.shape[0] 473 vvvv_p = vvvv_p.reshape(-1,nvir*nvir) 474 t[a:a+k] = np.dot(vvvv_p,t_amp).reshape(-1,nvir,nocc*nocc) 475 del vvvv_p 476 a += k 477 else: 478 raise Exception("Unknown vvvv type") 479 480 del t_amp 481 t = np.ascontiguousarray(t.transpose(2,0,1)).reshape(nocc, nocc, nvir, nvir) 482 483 return t 484 485 486def density_matrix(myadc, T=None): 487 488 if T is None: 489 T = RADCIP(myadc).get_trans_moments() 490 491 nocc = myadc._nocc 492 nvir = myadc._nvir 493 494 n_singles = nocc 495 n_doubles = nvir * nocc * nocc 496 497 s1 = 0 498 f1 = n_singles 499 500 T_doubles = T[:,n_singles:] 501 T_doubles = T_doubles.reshape(-1,nvir,nocc,nocc) 502 T_doubles_transpose = T_doubles.transpose(0,1,3,2).copy() 503 T_bab = (2/3)*T_doubles + (1/3)*T_doubles_transpose 504 505 T_aaa = T_bab - T_bab.transpose(0,1,3,2) 506 507 T_a = T[:,s1:f1] 508 T_bab = T_bab.reshape(-1,n_doubles) 509 T_aaa = T_aaa.reshape(-1,n_doubles) 510 511 dm = 2 * np.dot(T_a,T_a.T) + np.dot(T_aaa, T_aaa.T) + 2 * np.dot(T_bab, T_bab.T) 512 513 return dm 514 515 516def analyze(myadc): 517 518 header = ("\n*************************************************************" 519 "\n Eigenvector analysis summary" 520 "\n*************************************************************") 521 logger.info(myadc, header) 522 523 myadc.analyze_eigenvector() 524 525 if myadc.compute_properties: 526 527 header = ("\n*************************************************************" 528 "\n Spectroscopic factors analysis summary" 529 "\n*************************************************************") 530 logger.info(myadc, header) 531 532 myadc.analyze_spec_factor() 533 534 535def compute_dyson_mo(myadc): 536 537 X = myadc.X 538 539 if X is None: 540 nroots = myadc.U.shape[1] 541 P,X = myadc.get_properties(nroots) 542 543 nroots = X.shape[1] 544 dyson_mo = np.dot(myadc.mo_coeff,X) 545 546 return dyson_mo 547 548 549class RADC(lib.StreamObject): 550 '''Ground state calculations 551 552 Attributes: 553 verbose : int 554 Print level. Default value equals to :class:`Mole.verbose` 555 max_memory : float or int 556 Allowed memory in MB. Default value equals to :class:`Mole.max_memory` 557 incore_complete : bool 558 Avoid all I/O. Default is False. 559 method : string 560 nth-order ADC method. Options are : ADC(2), ADC(2)-X, ADC(3). Default is ADC(2). 561 562 >>> mol = gto.M(atom = 'H 0 0 0; F 0 0 1.1', basis = 'ccpvdz') 563 >>> mf = scf.RHF(mol).run() 564 >>> myadc = adc.RADC(mf).run() 565 566 Saved results 567 568 e_corr : float 569 MPn correlation correction 570 e_tot : float 571 Total energy (HF + correlation) 572 t1, t2 : 573 T amplitudes t1[i,a], t2[i,j,a,b] (i,j in occ, a,b in virt) 574 ''' 575 incore_complete = getattr(__config__, 'adc_radc_RADC_incore_complete', False) 576 async_io = getattr(__config__, 'adc_radc_RADC_async_io', True) 577 blkmin = getattr(__config__, 'adc_radc_RADC_blkmin', 4) 578 memorymin = getattr(__config__, 'adc_radc_RADC_memorymin', 2000) 579 580 def __init__(self, mf, frozen=0, mo_coeff=None, mo_occ=None): 581 from pyscf import gto 582 583 if 'dft' in str(mf.__module__): 584 raise NotImplementedError('DFT reference for UADC') 585 586 if mo_coeff is None: mo_coeff = mf.mo_coeff 587 if mo_occ is None: mo_occ = mf.mo_occ 588 589 self.mol = mf.mol 590 self._scf = mf 591 self.verbose = self.mol.verbose 592 self.stdout = self.mol.stdout 593 self.max_memory = mf.max_memory 594 595 self.max_space = getattr(__config__, 'adc_radc_RADC_max_space', 12) 596 self.max_cycle = getattr(__config__, 'adc_radc_RADC_max_cycle', 50) 597 self.conv_tol = getattr(__config__, 'adc_radc_RADC_conv_tol', 1e-12) 598 self.tol_residual = getattr(__config__, 'adc_radc_RADC_tol_res', 1e-6) 599 self.scf_energy = mf.e_tot 600 601 self.frozen = frozen 602 self.incore_complete = self.incore_complete or self.mol.incore_anyway 603 604 self.mo_coeff = mo_coeff 605 self.mo_occ = mo_occ 606 self.e_corr = None 607 self.t1 = None 608 self.t2 = None 609 self.imds = lambda:None 610 self._nocc = mf.mol.nelectron//2 611 self._nmo = mo_coeff.shape[1] 612 self._nvir = self._nmo - self._nocc 613 self.mo_energy = mf.mo_energy 614 self.chkfile = mf.chkfile 615 self.method = "adc(2)" 616 self.method_type = "ip" 617 self.with_df = None 618 self.compute_properties = True 619 self.evec_print_tol = 0.1 620 self.spec_factor_print_tol = 0.1 621 622 self.E = None 623 self.U = None 624 self.P = None 625 self.X = None 626 627 keys = set(('tol_residual','conv_tol', 'e_corr', 'method', 'mo_coeff', 628 'mol', 'mo_energy', 'max_memory', 'incore_complete', 629 'scf_energy', 'e_tot', 't1', 'frozen', 'chkfile', 630 'max_space', 't2', 'mo_occ', 'max_cycle')) 631 632 self._keys = set(self.__dict__.keys()).union(keys) 633 634 compute_amplitudes = compute_amplitudes 635 compute_energy = compute_energy 636 transform_integrals = radc_ao2mo.transform_integrals_incore 637 make_rdm1 = density_matrix 638 639 def dump_flags(self, verbose=None): 640 logger.info(self, '') 641 logger.info(self, '******** %s ********', self.__class__) 642 logger.info(self, 'max_space = %d', self.max_space) 643 logger.info(self, 'max_cycle = %d', self.max_cycle) 644 logger.info(self, 'conv_tol = %s', self.conv_tol) 645 logger.info(self, 'max_memory %d MB (current use %d MB)', 646 self.max_memory, lib.current_memory()[0]) 647 return self 648 649 def dump_flags_gs(self, verbose=None): 650 logger.info(self, '') 651 logger.info(self, '******** %s ********', self.__class__) 652 logger.info(self, 'max_memory %d MB (current use %d MB)', 653 self.max_memory, lib.current_memory()[0]) 654 return self 655 656 def kernel_gs(self): 657 assert(self.mo_coeff is not None) 658 assert(self.mo_occ is not None) 659 660 self.method = self.method.lower() 661 if self.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 662 raise NotImplementedError(self.method) 663 664 if self.verbose >= logger.WARN: 665 self.check_sanity() 666 self.dump_flags_gs() 667 668 nmo = self._nmo 669 nao = self.mo_coeff.shape[0] 670 nmo_pair = nmo * (nmo+1) // 2 671 nao_pair = nao * (nao+1) // 2 672 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6 673 mem_now = lib.current_memory()[0] 674 675 if getattr(self, 'with_df', None) or getattr(self._scf, 'with_df', None): 676 if getattr(self, 'with_df', None): 677 self.with_df = self.with_df 678 else: 679 self.with_df = self._scf.with_df 680 681 def df_transform(): 682 return radc_ao2mo.transform_integrals_df(self) 683 self.transform_integrals = df_transform 684 elif (self._scf._eri is None or 685 (mem_incore+mem_now >= self.max_memory and not self.incore_complete)): 686 def outcore_transform(): 687 return radc_ao2mo.transform_integrals_outcore(self) 688 self.transform_integrals = outcore_transform 689 690 eris = self.transform_integrals() 691 692 self.e_corr, self.t1, self.t2 = compute_amplitudes_energy(self, eris=eris, verbose=self.verbose) 693 self._finalize() 694 695 return self.e_corr, self.t1, self.t2 696 697 def kernel(self, nroots=1, guess=None, eris=None): 698 assert(self.mo_coeff is not None) 699 assert(self.mo_occ is not None) 700 701 self.method = self.method.lower() 702 if self.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 703 raise NotImplementedError(self.method) 704 705 if self.verbose >= logger.WARN: 706 self.check_sanity() 707 self.dump_flags_gs() 708 709 nmo = self._nmo 710 nao = self.mo_coeff.shape[0] 711 nmo_pair = nmo * (nmo+1) // 2 712 nao_pair = nao * (nao+1) // 2 713 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6 714 mem_now = lib.current_memory()[0] 715 716 if getattr(self, 'with_df', None) or getattr(self._scf, 'with_df', None): 717 if getattr(self, 'with_df', None): 718 self.with_df = self.with_df 719 else: 720 self.with_df = self._scf.with_df 721 722 def df_transform(): 723 return radc_ao2mo.transform_integrals_df(self) 724 self.transform_integrals = df_transform 725 elif (self._scf._eri is None or 726 (mem_incore+mem_now >= self.max_memory and not self.incore_complete)): 727 def outcore_transform(): 728 return radc_ao2mo.transform_integrals_outcore(self) 729 self.transform_integrals = outcore_transform 730 731 eris = self.transform_integrals() 732 733 self.e_corr, self.t1, self.t2 = compute_amplitudes_energy(self, eris=eris, verbose=self.verbose) 734 self._finalize() 735 736 self.method_type = self.method_type.lower() 737 if(self.method_type == "ea"): 738 e_exc, v_exc, spec_fac, x, adc_es = self.ea_adc(nroots=nroots, guess=guess, eris=eris) 739 740 elif(self.method_type == "ip"): 741 e_exc, v_exc, spec_fac, x, adc_es = self.ip_adc(nroots=nroots, guess=guess, eris=eris) 742 743 else: 744 raise NotImplementedError(self.method_type) 745 self._adc_es = adc_es 746 return e_exc, v_exc, spec_fac, x 747 748 def _finalize(self): 749 '''Hook for dumping results and clearing up the object.''' 750 logger.note(self, 'E_corr = %.8f', 751 self.e_corr) 752 return self 753 754 def ea_adc(self, nroots=1, guess=None, eris=None): 755 adc_es = RADCEA(self) 756 e_exc, v_exc, spec_fac, x = adc_es.kernel(nroots, guess, eris) 757 return e_exc, v_exc, spec_fac, x, adc_es 758 759 def ip_adc(self, nroots=1, guess=None, eris=None): 760 adc_es = RADCIP(self) 761 e_exc, v_exc, spec_fac, x = adc_es.kernel(nroots, guess, eris) 762 return e_exc, v_exc, spec_fac, x, adc_es 763 764 def density_fit(self, auxbasis=None, with_df = None): 765 if with_df is None: 766 self.with_df = df.DF(self._scf.mol) 767 self.with_df.max_memory = self.max_memory 768 self.with_df.stdout = self.stdout 769 self.with_df.verbose = self.verbose 770 if auxbasis is None: 771 self.with_df.auxbasis = self._scf.with_df.auxbasis 772 else: 773 self.with_df.auxbasis = auxbasis 774 else: 775 self.with_df = with_df 776 return self 777 778 def analyze(self): 779 self._adc_es.analyze() 780 781 def compute_dyson_mo(self): 782 return self._adc_es.compute_dyson_mo() 783 784 785def get_imds_ea(adc, eris=None): 786 787 cput0 = (logger.process_clock(), logger.perf_counter()) 788 log = logger.Logger(adc.stdout, adc.verbose) 789 790 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 791 raise NotImplementedError(adc.method) 792 793 method = adc.method 794 795 t1 = adc.t1 796 t2 = adc.t2 797 798 t1_2 = t1[0] 799 800 eris_ovvo = eris.ovvo 801 nocc = adc._nocc 802 nvir = adc._nvir 803 804 e_occ = adc.mo_energy[:nocc].copy() 805 e_vir = adc.mo_energy[nocc:].copy() 806 807 idn_vir = np.identity(nvir) 808 809 if eris is None: 810 eris = adc.transform_integrals() 811 812 # a-b block 813 # Zeroth-order terms 814 815 M_ab = lib.einsum('ab,a->ab', idn_vir, e_vir) 816 817 # Second-order terms 818 t2_1 = t2[0][:] 819 820 M_ab += lib.einsum('l,lmad,lmbd->ab',e_occ ,t2_1, t2_1,optimize=True) 821 M_ab -= lib.einsum('l,lmad,mlbd->ab',e_occ ,t2_1, t2_1,optimize=True) 822 M_ab -= lib.einsum('l,mlad,lmbd->ab',e_occ ,t2_1, t2_1,optimize=True) 823 M_ab += lib.einsum('l,mlad,mlbd->ab',e_occ ,t2_1, t2_1,optimize=True) 824 M_ab += lib.einsum('l,lmad,lmbd->ab',e_occ,t2_1, t2_1,optimize=True) 825 M_ab += lib.einsum('l,mlad,mlbd->ab',e_occ,t2_1, t2_1,optimize=True) 826 827 M_ab -= 0.5 * lib.einsum('d,lmad,lmbd->ab',e_vir,t2_1, t2_1,optimize=True) 828 M_ab += 0.5 * lib.einsum('d,lmad,mlbd->ab',e_vir,t2_1, t2_1,optimize=True) 829 M_ab += 0.5 * lib.einsum('d,mlad,lmbd->ab',e_vir,t2_1, t2_1,optimize=True) 830 M_ab -= 0.5 * lib.einsum('d,mlad,mlbd->ab',e_vir,t2_1, t2_1,optimize=True) 831 M_ab -= 0.5 * lib.einsum('d,lmad,lmbd->ab',e_vir,t2_1, t2_1,optimize=True) 832 M_ab -= 0.5 * lib.einsum('d,mlad,mlbd->ab',e_vir,t2_1, t2_1,optimize=True) 833 834 M_ab_t = lib.einsum('lmad,lmbd->ab', t2_1,t2_1, optimize=True) 835 M_ab -= 1 * lib.einsum('a,ab->ab',e_vir,M_ab_t,optimize=True) 836 M_ab -= 1 * lib.einsum('b,ab->ab',e_vir,M_ab_t,optimize=True) 837 838 M_ab_t = lib.einsum('lmad,mlbd->ab', t2_1,t2_1, optimize=True) 839 M_ab += 0.5 * lib.einsum('a,ab->ab',e_vir,M_ab_t,optimize=True) 840 M_ab += 0.5 * lib.einsum('b,ab->ab',e_vir,M_ab_t,optimize=True) 841 del M_ab_t 842 843 M_ab -= 0.5 * lib.einsum('lmad,lbdm->ab',t2_1, eris_ovvo,optimize=True) 844 M_ab += 0.5 * lib.einsum('mlad,lbdm->ab',t2_1, eris_ovvo,optimize=True) 845 M_ab += 0.5 * lib.einsum('lmad,ldbm->ab',t2_1, eris_ovvo,optimize=True) 846 M_ab -= 0.5 * lib.einsum('mlad,ldbm->ab',t2_1, eris_ovvo,optimize=True) 847 M_ab -= lib.einsum('lmad,lbdm->ab',t2_1, eris_ovvo,optimize=True) 848 849 M_ab -= 0.5 * lib.einsum('lmbd,ladm->ab',t2_1, eris_ovvo,optimize=True) 850 M_ab += 0.5 * lib.einsum('mlbd,ladm->ab',t2_1, eris_ovvo,optimize=True) 851 M_ab += 0.5 * lib.einsum('lmbd,ldam->ab',t2_1, eris_ovvo,optimize=True) 852 M_ab -= 0.5 * lib.einsum('mlbd,ldam->ab',t2_1, eris_ovvo,optimize=True) 853 M_ab -= lib.einsum('lmbd,ladm->ab',t2_1, eris_ovvo,optimize=True) 854 855 del t2_1 856 cput0 = log.timer_debug1("Completed M_ab second-order terms ADC(2) calculation", *cput0) 857 858 #Third-order terms 859 860 if(method =='adc(3)'): 861 862 eris_oovv = eris.oovv 863 eris_oooo = eris.oooo 864 865 if isinstance(eris.ovvv, type(None)): 866 chnk_size = radc_ao2mo.calculate_chunk_size(adc) 867 else: 868 chnk_size = nocc 869 a = 0 870 for p in range(0,nocc,chnk_size): 871 if getattr(adc, 'with_df', None): 872 eris_ovvv = dfadc.get_ovvv_df(adc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 873 else: 874 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 875 k = eris_ovvv.shape[0] 876 M_ab += 4. * lib.einsum('ld,ldab->ab',t1_2[a:a+k], eris_ovvv,optimize=True) 877 M_ab -= lib.einsum('ld,lbad->ab',t1_2[a:a+k], eris_ovvv,optimize=True) 878 M_ab -= lib.einsum('ld,ladb->ab',t1_2[a:a+k], eris_ovvv,optimize=True) 879 del eris_ovvv 880 a += k 881 882 cput0 = log.timer_debug1("Completed M_ab ovvv ADC(3) calculation", *cput0) 883 t2_2 = t2[1][:] 884 885 M_ab -= 0.5 * lib.einsum('lmad,lbdm->ab',t2_2, eris_ovvo,optimize=True) 886 M_ab += 0.5 * lib.einsum('mlad,lbdm->ab',t2_2, eris_ovvo,optimize=True) 887 M_ab += 0.5 * lib.einsum('lmad,ldbm->ab',t2_2, eris_ovvo,optimize=True) 888 M_ab -= 0.5 * lib.einsum('mlad,ldbm->ab',t2_2, eris_ovvo,optimize=True) 889 M_ab -= lib.einsum('lmad,lbdm->ab',t2_2, eris_ovvo,optimize=True) 890 891 M_ab -= 0.5 * lib.einsum('lmbd,ladm->ab',t2_2,eris_ovvo,optimize=True) 892 M_ab += 0.5 * lib.einsum('mlbd,ladm->ab',t2_2,eris_ovvo,optimize=True) 893 M_ab += 0.5 * lib.einsum('lmbd,ldam->ab',t2_2, eris_ovvo,optimize=True) 894 M_ab -= 0.5 * lib.einsum('mlbd,ldam->ab',t2_2, eris_ovvo,optimize=True) 895 M_ab -= lib.einsum('lmbd,ladm->ab',t2_2,eris_ovvo,optimize=True) 896 t2_1 = t2[0][:] 897 898 M_ab += lib.einsum('l,lmbd,lmad->ab',e_occ, t2_1, t2_2, optimize=True) 899 M_ab -= lib.einsum('l,lmbd,mlad->ab',e_occ, t2_1, t2_2, optimize=True) 900 M_ab -= lib.einsum('l,mlbd,lmad->ab',e_occ, t2_1, t2_2, optimize=True) 901 M_ab += lib.einsum('l,mlbd,mlad->ab',e_occ, t2_1, t2_2, optimize=True) 902 M_ab += lib.einsum('l,lmbd,lmad->ab',e_occ, t2_1, t2_2, optimize=True) 903 M_ab += lib.einsum('l,mlbd,mlad->ab',e_occ, t2_1, t2_2, optimize=True) 904 905 M_ab += lib.einsum('l,lmad,lmbd->ab',e_occ, t2_1, t2_2, optimize=True) 906 M_ab -= lib.einsum('l,lmad,mlbd->ab',e_occ, t2_1, t2_2, optimize=True) 907 M_ab -= lib.einsum('l,mlad,lmbd->ab',e_occ, t2_1, t2_2, optimize=True) 908 M_ab += lib.einsum('l,mlad,mlbd->ab',e_occ, t2_1, t2_2, optimize=True) 909 M_ab += lib.einsum('l,lmad,lmbd->ab',e_occ, t2_1, t2_2, optimize=True) 910 M_ab += lib.einsum('l,mlad,mlbd->ab',e_occ, t2_1, t2_2, optimize=True) 911 912 M_ab -= 0.5*lib.einsum('d,lmbd,lmad->ab', e_vir, t2_1 ,t2_2, optimize=True) 913 M_ab += 0.5*lib.einsum('d,lmbd,mlad->ab', e_vir, t2_1 ,t2_2, optimize=True) 914 M_ab += 0.5*lib.einsum('d,mlbd,lmad->ab', e_vir, t2_1 ,t2_2, optimize=True) 915 M_ab -= 0.5*lib.einsum('d,mlbd,mlad->ab', e_vir, t2_1 ,t2_2, optimize=True) 916 M_ab -= 0.5*lib.einsum('d,lmbd,lmad->ab', e_vir, t2_1 ,t2_2, optimize=True) 917 M_ab -= 0.5*lib.einsum('d,mlbd,mlad->ab', e_vir, t2_1 ,t2_2, optimize=True) 918 919 M_ab -= 0.5*lib.einsum('d,lmad,lmbd->ab', e_vir, t2_1, t2_2, optimize=True) 920 M_ab += 0.5*lib.einsum('d,lmad,mlbd->ab', e_vir, t2_1, t2_2, optimize=True) 921 M_ab += 0.5*lib.einsum('d,mlad,lmbd->ab', e_vir, t2_1, t2_2, optimize=True) 922 M_ab -= 0.5*lib.einsum('d,mlad,mlbd->ab', e_vir, t2_1, t2_2, optimize=True) 923 M_ab -= 0.5*lib.einsum('d,lmad,lmbd->ab', e_vir, t2_1, t2_2, optimize=True) 924 M_ab -= 0.5*lib.einsum('d,mlad,mlbd->ab', e_vir, t2_1, t2_2, optimize=True) 925 926 M_ab_t = lib.einsum('lmbd,lmad->ab', t2_1,t2_2, optimize=True) 927 M_ab -= 1. * lib.einsum('a,ab->ab',e_vir, M_ab_t, optimize=True) 928 M_ab -= 1. * lib.einsum('a,ba->ab',e_vir, M_ab_t, optimize=True) 929 M_ab -= 1. * lib.einsum('b,ab->ab',e_vir, M_ab_t, optimize=True) 930 M_ab -= 1. * lib.einsum('b,ba->ab',e_vir, M_ab_t, optimize=True) 931 del M_ab_t 932 933 M_ab_t_1 = lib.einsum('lmbd,mlad->ab', t2_1,t2_2, optimize=True) 934 del t2_2 935 M_ab += 0.5 * lib.einsum('a,ab->ab',e_vir, M_ab_t_1, optimize=True) 936 M_ab += 0.5 * lib.einsum('a,ba->ab',e_vir, M_ab_t_1, optimize=True) 937 M_ab += 0.5 * lib.einsum('b,ab->ab',e_vir, M_ab_t_1, optimize=True) 938 M_ab += 0.5 * lib.einsum('b,ba->ab',e_vir, M_ab_t_1, optimize=True) 939 del M_ab_t_1 940 941 log.timer_debug1("Starting the small integrals calculation") 942 temp_t2_v_1 = lib.einsum('lned,mlbd->nemb',t2_1, t2_1,optimize=True) 943 M_ab -= lib.einsum('nemb,nmae->ab',temp_t2_v_1, eris_oovv, optimize=True) 944 M_ab -= lib.einsum('mbne,nmae->ab',temp_t2_v_1, eris_oovv, optimize=True) 945 M_ab += lib.einsum('nemb,maen->ab',temp_t2_v_1, eris_ovvo, optimize=True) 946 M_ab += lib.einsum('mbne,maen->ab',temp_t2_v_1, eris_ovvo, optimize=True) 947 M_ab += lib.einsum('nemb,neam->ab',temp_t2_v_1, eris_ovvo, optimize=True) 948 M_ab -= lib.einsum('name,nmeb->ab',temp_t2_v_1, eris_oovv, optimize=True) 949 M_ab -= lib.einsum('mena,nmeb->ab',temp_t2_v_1, eris_oovv, optimize=True) 950 M_ab += 2. * lib.einsum('name,nbem->ab',temp_t2_v_1, eris_ovvo, optimize=True) 951 M_ab += 2. * lib.einsum('mena,nbem->ab',temp_t2_v_1, eris_ovvo, optimize=True) 952 M_ab += lib.einsum('nbme,mean->ab',temp_t2_v_1, eris_ovvo, optimize=True) 953 del temp_t2_v_1 954 955 temp_t2_v_2 = lib.einsum('nled,mlbd->nemb',t2_1, t2_1,optimize=True) 956 M_ab += 2. * lib.einsum('nemb,nmae->ab',temp_t2_v_2, eris_oovv, optimize=True) 957 M_ab -= 2. * lib.einsum('nemb,maen->ab',temp_t2_v_2, eris_ovvo, optimize=True) 958 M_ab -= lib.einsum('nemb,neam->ab',temp_t2_v_2, eris_ovvo, optimize=True) 959 M_ab += 2. * lib.einsum('mena,nmeb->ab',temp_t2_v_2, eris_oovv, optimize=True) 960 M_ab -= 4. * lib.einsum('mena,nbem->ab',temp_t2_v_2, eris_ovvo, optimize=True) 961 M_ab -= lib.einsum('nemb,neam->ab',temp_t2_v_2, eris_ovvo, optimize=True) 962 del temp_t2_v_2 963 964 temp_t2_v_3 = lib.einsum('lned,lmbd->nemb',t2_1, t2_1,optimize=True) 965 M_ab -= lib.einsum('nemb,maen->ab',temp_t2_v_3, eris_ovvo, optimize=True) 966 M_ab += 2. * lib.einsum('nemb,nmae->ab',temp_t2_v_3, eris_oovv, optimize=True) 967 M_ab += 2. * lib.einsum('mena,nmeb->ab',temp_t2_v_3, eris_oovv, optimize=True) 968 M_ab -= lib.einsum('mena,nbem->ab',temp_t2_v_3, eris_ovvo, optimize=True) 969 del temp_t2_v_3 970 971 temp_t2_v_4 = lib.einsum('lnae,nmde->lmad',t2_1, eris_oovv,optimize=True) 972 M_ab -= lib.einsum('mlbd,lmad->ab',t2_1, temp_t2_v_4,optimize=True) 973 M_ab += 2. * lib.einsum('lmbd,lmad->ab',t2_1, temp_t2_v_4,optimize=True) 974 del temp_t2_v_4 975 976 temp_t2_v_5 = lib.einsum('nlae,nmde->lamd',t2_1, eris_oovv,optimize=True) 977 M_ab += 2. * lib.einsum('mlbd,lamd->ab',t2_1, temp_t2_v_5, optimize=True) 978 M_ab -= lib.einsum('lmbd,lamd->ab',t2_1, temp_t2_v_5, optimize=True) 979 del temp_t2_v_5 980 981 temp_t2_v_6 = lib.einsum('lnae,nedm->ladm',t2_1, eris_ovvo,optimize=True) 982 M_ab += 2. * lib.einsum('mlbd,ladm->ab',t2_1, temp_t2_v_6, optimize=True) 983 M_ab -= 4. * lib.einsum('lmbd,ladm->ab',t2_1, temp_t2_v_6, optimize=True) 984 del temp_t2_v_6 985 986 temp_t2_v_7 = lib.einsum('nlae,nedm->ladm',t2_1, eris_ovvo,optimize=True) 987 M_ab -= lib.einsum('mlbd,ladm->ab',t2_1, temp_t2_v_7, optimize=True) 988 M_ab += 2. * lib.einsum('lmbd,ladm->ab',t2_1, temp_t2_v_7, optimize=True) 989 del temp_t2_v_7 990 991 temp_t2_v_8 = lib.einsum('lned,mled->mn',t2_1, t2_1,optimize=True) 992 M_ab += 2.* lib.einsum('mn,nmab->ab',temp_t2_v_8, eris_oovv, optimize=True) 993 M_ab -= lib.einsum('mn,nbam->ab', temp_t2_v_8, eris_ovvo, optimize=True) 994 del temp_t2_v_8 995 996 temp_t2_v_9 = lib.einsum('nled,mled->mn',t2_1, t2_1,optimize=True) 997 M_ab -= 4.* lib.einsum('mn,nmab->ab',temp_t2_v_9, eris_oovv, optimize=True) 998 M_ab += 2. * lib.einsum('mn,nbam->ab',temp_t2_v_9, eris_ovvo, optimize=True) 999 del temp_t2_v_9 1000 1001 temp_t2_v_10 = lib.einsum('noad,nmol->mlad',t2_1, eris_oooo,optimize=True) 1002 M_ab -= 0.25*lib.einsum('mlbd,mlad->ab',t2_1, temp_t2_v_10, optimize=True) 1003 M_ab += 0.25*lib.einsum('lmbd,mlad->ab',t2_1, temp_t2_v_10, optimize=True) 1004 M_ab += 0.25*lib.einsum('mlbd,lmad->ab',t2_1, temp_t2_v_10, optimize=True) 1005 M_ab -= 0.25*lib.einsum('lmbd,lmad->ab',t2_1, temp_t2_v_10, optimize=True) 1006 M_ab -= lib.einsum('mlbd,mlad->ab',t2_1, temp_t2_v_10, optimize=True) 1007 del temp_t2_v_10 1008 1009 temp_t2_v_11 = lib.einsum('onad,nmol->mlad',t2_1, eris_oooo,optimize=True) 1010 M_ab += 0.25*lib.einsum('mlbd,mlad->ab',t2_1, temp_t2_v_11, optimize=True) 1011 M_ab -= 0.25*lib.einsum('lmbd,mlad->ab',t2_1, temp_t2_v_11, optimize=True) 1012 M_ab -= 0.25*lib.einsum('mlbd,lmad->ab',t2_1, temp_t2_v_11, optimize=True) 1013 M_ab += 0.25*lib.einsum('lmbd,lmad->ab',t2_1, temp_t2_v_11, optimize=True) 1014 del temp_t2_v_11 1015 log.timer_debug1("Completed M_ab ADC(3) small integrals calculation") 1016 1017 log.timer_debug1("Starting M_ab vvvv ADC(3) calculation") 1018 1019 if isinstance(eris.vvvv, np.ndarray): 1020 temp_t2 = adc.imds.t2_1_vvvv 1021 M_ab -= 0.25*lib.einsum('mlaf,mlbf->ab',t2_1, temp_t2, optimize=True) 1022 M_ab += 0.25*lib.einsum('mlaf,lmbf->ab',t2_1, temp_t2, optimize=True) 1023 M_ab += 0.25*lib.einsum('lmaf,mlbf->ab',t2_1, temp_t2, optimize=True) 1024 M_ab -= 0.25*lib.einsum('lmaf,lmbf->ab',t2_1, temp_t2, optimize=True) 1025 M_ab += 0.25*lib.einsum('mlaf,mlfb->ab',t2_1, temp_t2, optimize=True) 1026 M_ab -= 0.25*lib.einsum('mlaf,lmfb->ab',t2_1, temp_t2, optimize=True) 1027 M_ab -= 0.25*lib.einsum('lmaf,mlfb->ab',t2_1, temp_t2, optimize=True) 1028 M_ab += 0.25*lib.einsum('lmaf,lmfb->ab',t2_1, temp_t2, optimize=True) 1029 M_ab -= lib.einsum('mlaf,mlbf->ab',t2_1, temp_t2, optimize=True) 1030 1031 M_ab -= 0.25*lib.einsum('mlad,mlbd->ab', temp_t2, t2_1, optimize=True) 1032 M_ab += 0.25*lib.einsum('mlad,lmbd->ab', temp_t2, t2_1, optimize=True) 1033 M_ab += 0.25*lib.einsum('lmad,mlbd->ab', temp_t2, t2_1, optimize=True) 1034 M_ab -= 0.25*lib.einsum('lmad,lmbd->ab', temp_t2, t2_1, optimize=True) 1035 M_ab -= lib.einsum('mlad,mlbd->ab', temp_t2, t2_1, optimize=True) 1036 1037 M_ab += 0.25*lib.einsum('lmad,mlbd->ab',temp_t2, t2_1, optimize=True) 1038 M_ab -= 0.25*lib.einsum('lmad,lmbd->ab',temp_t2, t2_1, optimize=True) 1039 M_ab -= 0.25*lib.einsum('mlad,mlbd->ab',temp_t2, t2_1, optimize=True) 1040 M_ab += 0.25*lib.einsum('mlad,lmbd->ab',temp_t2, t2_1, optimize=True) 1041 del temp_t2 1042 1043 eris_vvvv = eris.vvvv 1044 eris_vvvv = eris_vvvv.reshape(nvir,nvir,nvir,nvir) 1045 M_ab -= lib.einsum('mldf,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1046 M_ab += lib.einsum('mldf,lmed,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1047 M_ab += lib.einsum('lmdf,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1048 M_ab -= lib.einsum('lmdf,lmed,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1049 M_ab += 0.5*lib.einsum('mldf,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1050 M_ab -= 0.5*lib.einsum('mldf,lmed,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1051 M_ab -= 0.5*lib.einsum('lmdf,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1052 M_ab += 0.5*lib.einsum('lmdf,lmed,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1053 M_ab += 2.*lib.einsum('mlfd,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1054 M_ab -= lib.einsum('mlfd,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1055 eris_vvvv = eris_vvvv.reshape(nvir*nvir,nvir*nvir) 1056 1057 else: 1058 temp_t2_vvvv = adc.imds.t2_1_vvvv[:] 1059 M_ab -= 0.25*lib.einsum('mlaf,mlbf->ab',t2_1, temp_t2_vvvv, optimize=True) 1060 M_ab += 0.25*lib.einsum('mlaf,lmbf->ab',t2_1, temp_t2_vvvv, optimize=True) 1061 M_ab += 0.25*lib.einsum('lmaf,mlbf->ab',t2_1, temp_t2_vvvv, optimize=True) 1062 M_ab -= 0.25*lib.einsum('lmaf,lmbf->ab',t2_1, temp_t2_vvvv, optimize=True) 1063 1064 M_ab += 0.25*lib.einsum('mlaf,mlfb->ab',t2_1, temp_t2_vvvv, optimize=True) 1065 M_ab -= 0.25*lib.einsum('mlaf,lmfb->ab',t2_1, temp_t2_vvvv, optimize=True) 1066 M_ab -= 0.25*lib.einsum('lmaf,mlfb->ab',t2_1, temp_t2_vvvv, optimize=True) 1067 M_ab += 0.25*lib.einsum('lmaf,lmfb->ab',t2_1, temp_t2_vvvv, optimize=True) 1068 1069 M_ab -= lib.einsum('mlaf,mlbf->ab',t2_1, temp_t2_vvvv, optimize=True) 1070 1071 M_ab += 0.25*lib.einsum('lmad,mlbd->ab',temp_t2_vvvv, t2_1, optimize=True) 1072 M_ab -= 0.25*lib.einsum('lmad,lmbd->ab',temp_t2_vvvv, t2_1, optimize=True) 1073 M_ab -= 0.25*lib.einsum('mlad,mlbd->ab',temp_t2_vvvv, t2_1, optimize=True) 1074 M_ab += 0.25*lib.einsum('mlad,lmbd->ab',temp_t2_vvvv, t2_1, optimize=True) 1075 1076 M_ab -= 0.25*lib.einsum('mlad,mlbd->ab', temp_t2_vvvv, t2_1, optimize=True) 1077 M_ab += 0.25*lib.einsum('mlad,lmbd->ab', temp_t2_vvvv, t2_1, optimize=True) 1078 M_ab += 0.25*lib.einsum('lmad,mlbd->ab', temp_t2_vvvv, t2_1, optimize=True) 1079 M_ab -= 0.25*lib.einsum('lmad,lmbd->ab', temp_t2_vvvv, t2_1, optimize=True) 1080 M_ab -= lib.einsum('mlad,mlbd->ab', temp_t2_vvvv, t2_1, optimize=True) 1081 del temp_t2_vvvv 1082 1083 chnk_size = radc_ao2mo.calculate_chunk_size(adc) 1084 a = 0 1085 temp = np.zeros((nvir,nvir)) 1086 1087 if isinstance(eris.vvvv, list): 1088 for dataset in eris.vvvv: 1089 k = dataset.shape[0] 1090 eris_vvvv = dataset[:].reshape(-1,nvir,nvir,nvir) 1091 temp[a:a+k] -= lib.einsum('mldf,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1092 temp[a:a+k] += lib.einsum('mldf,lmed,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1093 temp[a:a+k] += lib.einsum('lmdf,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1094 temp[a:a+k] -= lib.einsum('lmdf,lmed,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1095 temp[a:a+k] += 0.5*lib.einsum('mldf,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1096 temp[a:a+k] -= 0.5*lib.einsum('mldf,lmed,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1097 temp[a:a+k] -= 0.5*lib.einsum('lmdf,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1098 temp[a:a+k] += 0.5*lib.einsum('lmdf,lmed,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1099 temp[a:a+k] += 2.*lib.einsum('mlfd,mled,aebf->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1100 temp[a:a+k] -= lib.einsum('mlfd,mled,aefb->ab',t2_1, t2_1, eris_vvvv, optimize=True) 1101 del eris_vvvv 1102 a += k 1103 else: 1104 for p in range(0,nvir,chnk_size): 1105 1106 vvvv = dfadc.get_vvvv_df(adc, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 1107 k = vvvv.shape[0] 1108 temp[a:a+k] -= lib.einsum('mldf,mled,aebf->ab',t2_1, t2_1, vvvv, optimize=True) 1109 temp[a:a+k] += lib.einsum('mldf,lmed,aebf->ab',t2_1, t2_1, vvvv, optimize=True) 1110 temp[a:a+k] += lib.einsum('lmdf,mled,aebf->ab',t2_1, t2_1, vvvv, optimize=True) 1111 temp[a:a+k] -= lib.einsum('lmdf,lmed,aebf->ab',t2_1, t2_1, vvvv, optimize=True) 1112 temp[a:a+k] += 0.5*lib.einsum('mldf,mled,aefb->ab',t2_1, t2_1, vvvv, optimize=True) 1113 temp[a:a+k] -= 0.5*lib.einsum('mldf,lmed,aefb->ab',t2_1, t2_1, vvvv, optimize=True) 1114 temp[a:a+k] -= 0.5*lib.einsum('lmdf,mled,aefb->ab',t2_1, t2_1, vvvv, optimize=True) 1115 temp[a:a+k] += 0.5*lib.einsum('lmdf,lmed,aefb->ab',t2_1, t2_1, vvvv, optimize=True) 1116 temp[a:a+k] += 2.*lib.einsum('mlfd,mled,aebf->ab',t2_1, t2_1, vvvv, optimize=True) 1117 temp[a:a+k] -= lib.einsum('mlfd,mled,aefb->ab',t2_1, t2_1, vvvv, optimize=True) 1118 del vvvv 1119 a += k 1120 1121 M_ab += temp 1122 del temp 1123 del t2_1 1124 1125 cput0 = log.timer_debug1("Completed M_ab ADC(3) calculation", *cput0) 1126 return M_ab 1127 1128 1129def get_imds_ip(adc, eris=None): 1130 1131 cput0 = (logger.process_clock(), logger.perf_counter()) 1132 log = logger.Logger(adc.stdout, adc.verbose) 1133 1134 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 1135 raise NotImplementedError(adc.method) 1136 1137 method = adc.method 1138 1139 t1 = adc.t1 1140 t2 = adc.t2 1141 1142 t1_2 = t1[0] 1143 1144 nocc = adc._nocc 1145 1146 e_occ = adc.mo_energy[:nocc] 1147 e_vir = adc.mo_energy[nocc:] 1148 1149 idn_occ = np.identity(nocc) 1150 1151 if eris is None: 1152 eris = adc.transform_integrals() 1153 1154 eris_ovvo = eris.ovvo 1155 1156 # i-j block 1157 # Zeroth-order terms 1158 1159 M_ij = lib.einsum('ij,j->ij', idn_occ ,e_occ) 1160 1161 # Second-order terms 1162 t2_1 = t2[0][:] 1163 1164 M_ij += lib.einsum('d,ilde,jlde->ij',e_vir,t2_1, t2_1, optimize=True) 1165 M_ij -= lib.einsum('d,ilde,ljde->ij',e_vir,t2_1, t2_1, optimize=True) 1166 M_ij -= lib.einsum('d,lide,jlde->ij',e_vir,t2_1, t2_1, optimize=True) 1167 M_ij += lib.einsum('d,lide,ljde->ij',e_vir,t2_1, t2_1, optimize=True) 1168 M_ij += lib.einsum('d,ilde,jlde->ij',e_vir,t2_1, t2_1, optimize=True) 1169 M_ij += lib.einsum('d,iled,jled->ij',e_vir,t2_1, t2_1, optimize=True) 1170 1171 M_ij -= 0.5 * lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_1, optimize=True) 1172 M_ij += 0.5 * lib.einsum('l,ilde,ljde->ij',e_occ,t2_1, t2_1, optimize=True) 1173 M_ij += 0.5 * lib.einsum('l,lide,jlde->ij',e_occ,t2_1, t2_1, optimize=True) 1174 M_ij -= 0.5 * lib.einsum('l,lide,ljde->ij',e_occ,t2_1, t2_1, optimize=True) 1175 M_ij -= 0.5*lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_1, optimize=True) 1176 M_ij -= 0.5*lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_1, optimize=True) 1177 1178 M_ij_t = lib.einsum('ilde,jlde->ij', t2_1,t2_1, optimize=True) 1179 M_ij -= lib.einsum('i,ij->ij',e_occ,M_ij_t, optimize=True) 1180 M_ij -= lib.einsum('j,ij->ij',e_occ,M_ij_t, optimize=True) 1181 1182 M_ij_t = lib.einsum('ilde,ljde->ij', t2_1,t2_1, optimize=True) 1183 M_ij += 0.5 * lib.einsum('i,ij->ij',e_occ,M_ij_t, optimize=True) 1184 M_ij += 0.5 * lib.einsum('j,ij->ij',e_occ,M_ij_t, optimize=True) 1185 del M_ij_t 1186 1187 M_ij += 0.5 * lib.einsum('ilde,jdel->ij',t2_1, eris_ovvo,optimize=True) 1188 M_ij -= 0.5 * lib.einsum('lide,jdel->ij',t2_1, eris_ovvo,optimize=True) 1189 M_ij -= 0.5 * lib.einsum('ilde,jedl->ij',t2_1, eris_ovvo,optimize=True) 1190 M_ij += 0.5 * lib.einsum('lide,jedl->ij',t2_1, eris_ovvo,optimize=True) 1191 M_ij += lib.einsum('ilde,jdel->ij',t2_1, eris_ovvo,optimize=True) 1192 1193 M_ij += 0.5 * lib.einsum('jlde,idel->ij',t2_1, eris_ovvo,optimize=True) 1194 M_ij -= 0.5 * lib.einsum('ljde,idel->ij',t2_1, eris_ovvo,optimize=True) 1195 M_ij -= 0.5 * lib.einsum('jlde,ldei->ij',t2_1, eris_ovvo,optimize=True) 1196 M_ij += 0.5 * lib.einsum('ljde,ldei->ij',t2_1, eris_ovvo,optimize=True) 1197 M_ij += lib.einsum('jlde,idel->ij',t2_1, eris_ovvo,optimize=True) 1198 1199 del t2_1 1200 cput0 = log.timer_debug1("Completed M_ij second-order terms ADC(2) calculation", *cput0) 1201 # Third-order terms 1202 1203 if (method == "adc(3)"): 1204 1205 eris_oovv = eris.oovv 1206 eris_ovoo = eris.ovoo 1207 eris_oooo = eris.oooo 1208 1209 M_ij += lib.einsum('ld,ldji->ij',t1_2, eris_ovoo,optimize=True) 1210 M_ij -= lib.einsum('ld,jdli->ij',t1_2, eris_ovoo,optimize=True) 1211 M_ij += lib.einsum('ld,ldji->ij',t1_2, eris_ovoo,optimize=True) 1212 1213 M_ij += lib.einsum('ld,ldij->ij',t1_2, eris_ovoo,optimize=True) 1214 M_ij -= lib.einsum('ld,idlj->ij',t1_2, eris_ovoo,optimize=True) 1215 M_ij += lib.einsum('ld,ldij->ij',t1_2, eris_ovoo,optimize=True) 1216 t2_2 = t2[1][:] 1217 1218 M_ij += 0.5* lib.einsum('ilde,jdel->ij',t2_2, eris_ovvo,optimize=True) 1219 M_ij -= 0.5* lib.einsum('lide,jdel->ij',t2_2, eris_ovvo,optimize=True) 1220 M_ij -= 0.5* lib.einsum('ilde,jedl->ij',t2_2, eris_ovvo,optimize=True) 1221 M_ij += 0.5* lib.einsum('lide,jedl->ij',t2_2, eris_ovvo,optimize=True) 1222 M_ij += lib.einsum('ilde,jdel->ij',t2_2, eris_ovvo,optimize=True) 1223 1224 M_ij += 0.5* lib.einsum('jlde,ledi->ij',t2_2, eris_ovvo,optimize=True) 1225 M_ij -= 0.5* lib.einsum('ljde,ledi->ij',t2_2, eris_ovvo,optimize=True) 1226 M_ij -= 0.5* lib.einsum('jlde,iedl->ij',t2_2, eris_ovvo,optimize=True) 1227 M_ij += 0.5* lib.einsum('ljde,iedl->ij',t2_2, eris_ovvo,optimize=True) 1228 M_ij += lib.einsum('jlde,ledi->ij',t2_2, eris_ovvo,optimize=True) 1229 t2_1 = t2[0][:] 1230 1231 M_ij += lib.einsum('d,ilde,jlde->ij',e_vir,t2_1, t2_2,optimize=True) 1232 M_ij -= lib.einsum('d,ilde,ljde->ij',e_vir,t2_1, t2_2,optimize=True) 1233 M_ij -= lib.einsum('d,lide,jlde->ij',e_vir,t2_1, t2_2,optimize=True) 1234 M_ij += lib.einsum('d,lide,ljde->ij',e_vir,t2_1, t2_2,optimize=True) 1235 M_ij += lib.einsum('d,ilde,jlde->ij',e_vir,t2_1, t2_2,optimize=True) 1236 M_ij += lib.einsum('d,iled,jled->ij',e_vir,t2_1, t2_2,optimize=True) 1237 1238 M_ij += lib.einsum('d,jlde,ilde->ij',e_vir,t2_1, t2_2,optimize=True) 1239 M_ij -= lib.einsum('d,jlde,lide->ij',e_vir,t2_1, t2_2,optimize=True) 1240 M_ij -= lib.einsum('d,ljde,ilde->ij',e_vir,t2_1, t2_2,optimize=True) 1241 M_ij += lib.einsum('d,ljde,lide->ij',e_vir,t2_1, t2_2,optimize=True) 1242 M_ij += lib.einsum('d,jlde,ilde->ij',e_vir,t2_1, t2_2,optimize=True) 1243 M_ij += lib.einsum('d,jled,iled->ij',e_vir,t2_1, t2_2,optimize=True) 1244 1245 M_ij -= 0.5 * lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_2,optimize=True) 1246 M_ij += 0.5 * lib.einsum('l,ilde,ljde->ij',e_occ,t2_1, t2_2,optimize=True) 1247 M_ij += 0.5 * lib.einsum('l,lide,jlde->ij',e_occ,t2_1, t2_2,optimize=True) 1248 M_ij -= 0.5 * lib.einsum('l,lide,ljde->ij',e_occ,t2_1, t2_2,optimize=True) 1249 M_ij -= 0.5*lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_2,optimize=True) 1250 M_ij -= 0.5*lib.einsum('l,ilde,jlde->ij',e_occ,t2_1, t2_2,optimize=True) 1251 1252 M_ij -= 0.5 * lib.einsum('l,jlde,ilde->ij',e_occ,t2_1, t2_2,optimize=True) 1253 M_ij += 0.5 * lib.einsum('l,jlde,lide->ij',e_occ,t2_1, t2_2,optimize=True) 1254 M_ij += 0.5 * lib.einsum('l,ljde,ilde->ij',e_occ,t2_1, t2_2,optimize=True) 1255 M_ij -= 0.5 * lib.einsum('l,ljde,lide->ij',e_occ,t2_1, t2_2,optimize=True) 1256 M_ij -= 0.5*lib.einsum('l,jlde,ilde->ij',e_occ,t2_1, t2_2,optimize=True) 1257 M_ij -= 0.5*lib.einsum('l,jlde,ilde->ij',e_occ,t2_1, t2_2,optimize=True) 1258 1259 M_ij_t = lib.einsum('ilde,jlde->ij', t2_1,t2_2, optimize=True) 1260 M_ij -= 1. * lib.einsum('i,ij->ij',e_occ, M_ij_t, optimize=True) 1261 M_ij -= 1. * lib.einsum('i,ji->ij',e_occ, M_ij_t, optimize=True) 1262 M_ij -= 1. * lib.einsum('j,ij->ij',e_occ, M_ij_t, optimize=True) 1263 M_ij -= 1. * lib.einsum('j,ji->ij',e_occ, M_ij_t, optimize=True) 1264 del M_ij_t 1265 1266 M_ij_t_1 = lib.einsum('ilde,ljde->ij', t2_1,t2_2, optimize=True) 1267 del t2_2 1268 M_ij += 0.5 * lib.einsum('i,ij->ij',e_occ, M_ij_t_1, optimize=True) 1269 M_ij += 0.5 * lib.einsum('i,ji->ij',e_occ, M_ij_t_1, optimize=True) 1270 M_ij += 0.5 * lib.einsum('j,ij->ij',e_occ, M_ij_t_1, optimize=True) 1271 M_ij += 0.5 * lib.einsum('j,ji->ij',e_occ, M_ij_t_1, optimize=True) 1272 del M_ij_t_1 1273 1274 temp_t2_vvvv = adc.imds.t2_1_vvvv[:] 1275 M_ij += 0.25*lib.einsum('ilde,jlde->ij',t2_1, temp_t2_vvvv, optimize = True) 1276 M_ij -= 0.25*lib.einsum('ilde,ljde->ij',t2_1, temp_t2_vvvv, optimize = True) 1277 M_ij -= 0.25*lib.einsum('lide,jlde->ij',t2_1, temp_t2_vvvv, optimize = True) 1278 M_ij += 0.25*lib.einsum('lide,ljde->ij',t2_1, temp_t2_vvvv, optimize = True) 1279 M_ij -= 0.25*lib.einsum('ilde,jled->ij',t2_1, temp_t2_vvvv, optimize = True) 1280 M_ij += 0.25*lib.einsum('ilde,ljed->ij',t2_1, temp_t2_vvvv, optimize = True) 1281 M_ij += 0.25*lib.einsum('lide,jled->ij',t2_1, temp_t2_vvvv, optimize = True) 1282 M_ij -= 0.25*lib.einsum('lide,ljed->ij',t2_1, temp_t2_vvvv, optimize = True) 1283 M_ij +=lib.einsum('ilde,jlde->ij',t2_1, temp_t2_vvvv, optimize = True) 1284 del temp_t2_vvvv 1285 1286 log.timer_debug1("Starting the small integrals calculation") 1287 temp_t2_v_1 = lib.einsum('lmde,jldf->mejf',t2_1, t2_1,optimize=True) 1288 M_ij -= 2 * lib.einsum('mejf,mefi->ij',temp_t2_v_1, eris_ovvo,optimize = True) 1289 M_ij -= 2 * lib.einsum('jfme,mefi->ij',temp_t2_v_1, eris_ovvo,optimize = True) 1290 M_ij += lib.einsum('mejf,mife->ij',temp_t2_v_1, eris_oovv,optimize = True) 1291 M_ij += lib.einsum('jfme,mife->ij',temp_t2_v_1, eris_oovv,optimize = True) 1292 M_ij -= 2 * lib.einsum('meif,mefj->ij',temp_t2_v_1, eris_ovvo ,optimize = True) 1293 M_ij -= 2 * lib.einsum('ifme,mefj->ij',temp_t2_v_1, eris_ovvo ,optimize = True) 1294 M_ij += lib.einsum('meif,mjfe->ij',temp_t2_v_1, eris_oovv ,optimize = True) 1295 M_ij += lib.einsum('ifme,mjfe->ij',temp_t2_v_1, eris_oovv ,optimize = True) 1296 del temp_t2_v_1 1297 1298 temp_t2_v_2 = lib.einsum('lmde,ljdf->mejf',t2_1, t2_1,optimize=True) 1299 M_ij += 4 * lib.einsum('mejf,mefi->ij',temp_t2_v_2, eris_ovvo,optimize = True) 1300 M_ij += 4 * lib.einsum('meif,mefj->ij',temp_t2_v_2, eris_ovvo,optimize = True) 1301 M_ij -= 2 * lib.einsum('meif,mjfe->ij',temp_t2_v_2, eris_oovv,optimize = True) 1302 M_ij -= 2 * lib.einsum('mejf,mife->ij',temp_t2_v_2, eris_oovv,optimize = True) 1303 del temp_t2_v_2 1304 1305 temp_t2_v_3 = lib.einsum('mlde,jldf->mejf',t2_1, t2_1,optimize=True) 1306 M_ij += lib.einsum('mejf,mefi->ij',temp_t2_v_3, eris_ovvo,optimize = True) 1307 M_ij += lib.einsum('meif,mefj->ij',temp_t2_v_3, eris_ovvo,optimize = True) 1308 M_ij -= 2 *lib.einsum('meif,mjfe->ij',temp_t2_v_3, eris_oovv,optimize = True) 1309 M_ij -= 2 * lib.einsum('mejf,mife->ij',temp_t2_v_3, eris_oovv,optimize = True) 1310 del temp_t2_v_3 1311 1312 temp_t2_v_4 = lib.einsum('ilde,lmfe->idmf',t2_1, eris_oovv,optimize=True) 1313 M_ij -= 2 * lib.einsum('idmf,jmdf->ij',temp_t2_v_4, t2_1, optimize = True) 1314 M_ij += lib.einsum('idmf,mjdf->ij',temp_t2_v_4, t2_1, optimize = True) 1315 del temp_t2_v_4 1316 1317 temp_t2_v_5 = lib.einsum('lide,lmfe->idmf',t2_1, eris_oovv,optimize=True) 1318 M_ij += lib.einsum('idmf,jmdf->ij',temp_t2_v_5, t2_1, optimize = True) 1319 M_ij -= 2 * lib.einsum('idmf,mjdf->ij',temp_t2_v_5, t2_1, optimize = True) 1320 del temp_t2_v_5 1321 1322 temp_t2_v_6 = lib.einsum('ilde,lefm->idfm',t2_1, eris_ovvo,optimize=True) 1323 M_ij += 4 * lib.einsum('idfm,jmdf->ij',temp_t2_v_6, t2_1,optimize = True) 1324 M_ij -= 2 * lib.einsum('idfm,mjdf->ij',temp_t2_v_6, t2_1,optimize = True) 1325 del temp_t2_v_6 1326 1327 temp_t2_v_7 = lib.einsum('lide,lefm->idfm',t2_1, eris_ovvo,optimize=True) 1328 M_ij -= 2 * lib.einsum('idfm,jmdf->ij',temp_t2_v_7, t2_1,optimize = True) 1329 M_ij += lib.einsum('idfm,mjdf->ij',temp_t2_v_7, t2_1,optimize = True) 1330 del temp_t2_v_7 1331 1332 temp_t2_v_8 = lib.einsum('lmdf,lmde->fe',t2_1, t2_1,optimize=True) 1333 M_ij += 3 *lib.einsum('fe,jief->ij',temp_t2_v_8, eris_oovv, optimize = True) 1334 M_ij -= 1.5 *lib.einsum('fe,jfei->ij',temp_t2_v_8, eris_ovvo, optimize = True) 1335 M_ij += lib.einsum('ef,jief->ij',temp_t2_v_8, eris_oovv, optimize = True) 1336 M_ij -= 0.5 * lib.einsum('ef,jfei->ij',temp_t2_v_8, eris_ovvo, optimize = True) 1337 del temp_t2_v_8 1338 1339 temp_t2_v_9 = lib.einsum('lmdf,mlde->fe',t2_1, t2_1,optimize=True) 1340 M_ij -= 1.0 * lib.einsum('fe,jief->ij',temp_t2_v_9, eris_oovv, optimize = True) 1341 M_ij -= 1.0 * lib.einsum('ef,jief->ij',temp_t2_v_9, eris_oovv, optimize = True) 1342 M_ij += 0.5 * lib.einsum('fe,jfei->ij',temp_t2_v_9, eris_ovvo, optimize = True) 1343 M_ij += 0.5 * lib.einsum('ef,jfei->ij',temp_t2_v_9, eris_ovvo, optimize = True) 1344 del temp_t2_v_9 1345 1346 temp_t2_v_10 = lib.einsum('lnde,lmde->nm',t2_1, t2_1,optimize=True) 1347 M_ij -= 3.0 * lib.einsum('nm,jinm->ij',temp_t2_v_10, eris_oooo, optimize = True) 1348 M_ij -= 1.0 * lib.einsum('mn,jinm->ij',temp_t2_v_10, eris_oooo, optimize = True) 1349 M_ij += 1.5 * lib.einsum('nm,jmni->ij',temp_t2_v_10, eris_oooo, optimize = True) 1350 M_ij += 0.5 * lib.einsum('mn,jmni->ij',temp_t2_v_10, eris_oooo, optimize = True) 1351 del temp_t2_v_10 1352 1353 temp_t2_v_11 = lib.einsum('lnde,mlde->nm',t2_1, t2_1,optimize=True) 1354 M_ij += 1.0 * lib.einsum('nm,jinm->ij',temp_t2_v_11, eris_oooo, optimize = True) 1355 M_ij -= 0.5 * lib.einsum('nm,jmni->ij',temp_t2_v_11, eris_oooo, optimize = True) 1356 M_ij -= 0.5 * lib.einsum('mn,jmni->ij',temp_t2_v_11, eris_oooo, optimize = True) 1357 M_ij += 1.0 * lib.einsum('mn,jinm->ij',temp_t2_v_11, eris_oooo, optimize = True) 1358 del temp_t2_v_11 1359 1360 temp_t2_v_12 = lib.einsum('inde,lmde->inlm',t2_1, t2_1,optimize=True) 1361 M_ij += 1.25 * lib.einsum('inlm,jlnm->ij',temp_t2_v_12, eris_oooo, optimize = True) 1362 M_ij += 0.25 * lib.einsum('lmin,jlnm->ij',temp_t2_v_12, eris_oooo, optimize = True) 1363 M_ij -= 0.25 * lib.einsum('inlm,jmnl->ij',temp_t2_v_12, eris_oooo, optimize = True) 1364 M_ij -= 0.25 * lib.einsum('lmin,jmnl->ij',temp_t2_v_12, eris_oooo, optimize = True) 1365 1366 M_ij += 0.25 * lib.einsum('inlm,jlnm->ji',temp_t2_v_12, eris_oooo, optimize = True) 1367 M_ij -= 0.25 * lib.einsum('inlm,lnmj->ji',temp_t2_v_12, eris_oooo, optimize = True) 1368 M_ij += 1.00 * lib.einsum('inlm,ljmn->ji',temp_t2_v_12, eris_oooo, optimize = True) 1369 M_ij -= 0.25 * lib.einsum('lmin,lnmj->ji',temp_t2_v_12, eris_oooo, optimize = True) 1370 M_ij += 0.25 * lib.einsum('lmin,ljmn->ji',temp_t2_v_12, eris_oooo, optimize = True) 1371 del temp_t2_v_12 1372 1373 temp_t2_v_13 = lib.einsum('inde,mlde->inml',t2_1, t2_1,optimize=True) 1374 M_ij -= 0.25 * lib.einsum('inml,jlnm->ij',temp_t2_v_13, eris_oooo, optimize = True) 1375 M_ij -= 0.25 * lib.einsum('mlin,jlnm->ij',temp_t2_v_13, eris_oooo, optimize = True) 1376 M_ij += 0.25 * lib.einsum('inml,jmnl->ij',temp_t2_v_13, eris_oooo, optimize = True) 1377 M_ij += 0.25 * lib.einsum('mlin,jmnl->ij',temp_t2_v_13, eris_oooo, optimize = True) 1378 1379 M_ij -= 0.25 * lib.einsum('inml,jlnm->ji',temp_t2_v_13, eris_oooo, optimize = True) 1380 M_ij += 0.25 * lib.einsum('inml,lnmj->ji',temp_t2_v_13, eris_oooo, optimize = True) 1381 1382 M_ij -= 0.25 * lib.einsum('inml,ljmn->ji',temp_t2_v_13, eris_oooo, optimize = True) 1383 M_ij += 0.25 * lib.einsum('inml,lnmj->ji',temp_t2_v_13, eris_oooo, optimize = True) 1384 del temp_t2_v_13 1385 del t2_1 1386 1387 cput0 = log.timer_debug1("Completed M_ij ADC(n) calculation", *cput0) 1388 return M_ij 1389 1390 1391def ea_adc_diag(adc,M_ab=None,eris=None): 1392 1393 log = logger.Logger(adc.stdout, adc.verbose) 1394 1395 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 1396 raise NotImplementedError(adc.method) 1397 1398 if M_ab is None: 1399 M_ab = adc.get_imds() 1400 1401 nocc = adc._nocc 1402 nvir = adc._nvir 1403 1404 n_singles = nvir 1405 n_doubles = nocc * nvir * nvir 1406 1407 dim = n_singles + n_doubles 1408 1409 e_occ = adc.mo_energy[:nocc] 1410 e_vir = adc.mo_energy[nocc:] 1411 1412 s1 = 0 1413 f1 = n_singles 1414 s2 = f1 1415 f2 = s2 + n_doubles 1416 1417 d_ab = e_vir[:,None] + e_vir 1418 d_i = e_occ[:,None] 1419 D_n = -d_i + d_ab.reshape(-1) 1420 D_iab = D_n.reshape(-1) 1421 1422 diag = np.zeros(dim) 1423 1424 # Compute precond in p1-p1 block 1425 1426 M_ab_diag = np.diagonal(M_ab) 1427 diag[s1:f1] = M_ab_diag.copy() 1428 1429 # Compute precond in 2p1h-2p1h block 1430 1431 diag[s2:f2] = D_iab 1432 del D_iab 1433 1434# ###### Additional terms for the preconditioner #### 1435# 1436# if (method == "adc(2)-x" or method == "adc(3)"): 1437# 1438# if eris is None: 1439# eris = adc.transform_integrals() 1440# 1441# #TODO Implement this for out-of-core and density-fitted algorithms 1442# if isinstance(eris.vvvv, np.ndarray): 1443# 1444# eris_oovv = eris.oovv 1445# eris_ovvo = eris.ovvo 1446# eris_vvvv = eris.vvvv 1447# 1448# temp = np.zeros((nocc, eris_vvvv.shape[0])) 1449# temp[:] += np.diag(eris_vvvv) 1450# diag[s2:f2] += temp.reshape(-1) 1451# 1452# eris_ovov_p = np.ascontiguousarray(eris_oovv[:].transpose(0,2,1,3)) 1453# eris_ovov_p = eris_ovov_p.reshape(nocc*nvir, nocc*nvir) 1454# 1455# temp = np.zeros((nvir, nocc, nvir)) 1456# temp[:] += np.diagonal(eris_ovov_p).reshape(nocc, nvir) 1457# temp = np.ascontiguousarray(temp.transpose(1,0,2)) 1458# diag[s2:f2] += -temp.reshape(-1) 1459# 1460# eris_ovov_p = np.ascontiguousarray(eris_oovv[:].transpose(0,2,1,3)) 1461# eris_ovov_p = eris_ovov_p.reshape(nocc*nvir, nocc*nvir) 1462# 1463# temp = np.zeros((nvir, nocc, nvir)) 1464# temp[:] += np.diagonal(eris_ovov_p).reshape(nocc, nvir) 1465# temp = np.ascontiguousarray(temp.transpose(1,2,0)) 1466# diag[s2:f2] += -temp.reshape(-1) 1467# else: 1468# raise Exception("Precond not available for out-of-core and density-fitted algo") 1469 1470 log.timer_debug1("Completed ea_diag calculation") 1471 return diag 1472 1473 1474def ip_adc_diag(adc,M_ij=None,eris=None): 1475 1476 log = logger.Logger(adc.stdout, adc.verbose) 1477 1478 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 1479 raise NotImplementedError(adc.method) 1480 1481 if M_ij is None: 1482 M_ij = adc.get_imds() 1483 1484 nocc = adc._nocc 1485 nvir = adc._nvir 1486 1487 n_singles = nocc 1488 n_doubles = nvir * nocc * nocc 1489 1490 dim = n_singles + n_doubles 1491 1492 e_occ = adc.mo_energy[:nocc] 1493 e_vir = adc.mo_energy[nocc:] 1494 1495 s1 = 0 1496 f1 = n_singles 1497 s2 = f1 1498 f2 = s2 + n_doubles 1499 1500 d_ij = e_occ[:,None] + e_occ 1501 d_a = e_vir[:,None] 1502 D_n = -d_a + d_ij.reshape(-1) 1503 D_aij = D_n.reshape(-1) 1504 1505 diag = np.zeros(dim) 1506 1507 # Compute precond in h1-h1 block 1508 M_ij_diag = np.diagonal(M_ij) 1509 diag[s1:f1] = M_ij_diag.copy() 1510 1511 # Compute precond in 2p1h-2p1h block 1512 1513 diag[s2:f2] = D_aij.copy() 1514 1515# ###### Additional terms for the preconditioner #### 1516# if (method == "adc(2)-x" or method == "adc(3)"): 1517# 1518# if eris is None: 1519# eris = adc.transform_integrals() 1520# 1521# if isinstance(eris.vvvv, np.ndarray): 1522# 1523# eris_oooo = eris.oooo 1524# eris_oovv = eris.oovv 1525# eris_ovvo = eris.ovvo 1526# 1527# eris_oooo_p = np.ascontiguousarray(eris_oooo.transpose(0,2,1,3)) 1528# eris_oooo_p = eris_oooo_p.reshape(nocc*nocc, nocc*nocc) 1529# 1530# temp = np.zeros((nvir, eris_oooo_p.shape[0])) 1531# temp[:] += np.diag(eris_oooo_p) 1532# diag[s2:f2] += -temp.reshape(-1) 1533# 1534# eris_ovov_p = np.ascontiguousarray(eris_oovv.transpose(0,2,1,3)) 1535# eris_ovov_p = eris_ovov_p.reshape(nocc*nvir, nocc*nvir) 1536# 1537# temp = np.zeros((nocc, nocc, nvir)) 1538# temp[:] += np.diagonal(eris_ovov_p).reshape(nocc, nvir) 1539# temp = np.ascontiguousarray(temp.transpose(2,1,0)) 1540# diag[s2:f2] += temp.reshape(-1) 1541# 1542# eris_ovov_p = np.ascontiguousarray(eris_oovv.transpose(0,2,1,3)) 1543# eris_ovov_p = eris_ovov_p.reshape(nocc*nvir, nocc*nvir) 1544# 1545# temp = np.zeros((nocc, nocc, nvir)) 1546# temp[:] += np.diagonal(eris_ovov_p).reshape(nocc, nvir) 1547# temp = np.ascontiguousarray(temp.transpose(2,0,1)) 1548# diag[s2:f2] += temp.reshape(-1) 1549# else: 1550# raise Exception("Precond not available for out-of-core and density-fitted algo") 1551 1552 diag = -diag 1553 log.timer_debug1("Completed ea_diag calculation") 1554 1555 return diag 1556 1557 1558def ea_contract_r_vvvv(myadc,r2,vvvv): 1559 1560 nocc = myadc._nocc 1561 nvir = myadc._nvir 1562 1563 r2_vvvv = np.zeros((nocc,nvir,nvir)) 1564 r2 = np.ascontiguousarray(r2.reshape(nocc,-1)) 1565 chnk_size = radc_ao2mo.calculate_chunk_size(myadc) 1566 1567 a = 0 1568 if isinstance(vvvv, list): 1569 for dataset in vvvv: 1570 k = dataset.shape[0] 1571 dataset = dataset[:].reshape(-1,nvir*nvir) 1572 r2_vvvv[:,a:a+k] = np.dot(r2,dataset.T).reshape(nocc,-1,nvir) 1573 del dataset 1574 a += k 1575 elif getattr(myadc, 'with_df', None): 1576 for p in range(0,nvir,chnk_size): 1577 vvvv_p = dfadc.get_vvvv_df(myadc, vvvv, p, chnk_size) 1578 k = vvvv_p.shape[0] 1579 vvvv_p = vvvv_p.reshape(-1,nvir*nvir) 1580 r2_vvvv[:,a:a+k] = np.dot(r2,vvvv_p.T).reshape(nocc,-1,nvir) 1581 del vvvv_p 1582 a += k 1583 else: 1584 raise Exception("Unknown vvvv type") 1585 1586 r2_vvvv = r2_vvvv.reshape(-1) 1587 1588 return r2_vvvv 1589 1590 1591def ea_adc_matvec(adc, M_ab=None, eris=None): 1592 1593 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 1594 raise NotImplementedError(adc.method) 1595 1596 method = adc.method 1597 1598 1599 nocc = adc._nocc 1600 nvir = adc._nvir 1601 1602 n_singles = nvir 1603 n_doubles = nocc * nvir * nvir 1604 1605 dim = n_singles + n_doubles 1606 1607 e_occ = adc.mo_energy[:nocc] 1608 e_vir = adc.mo_energy[nocc:] 1609 1610 if eris is None: 1611 eris = adc.transform_integrals() 1612 1613 s1 = 0 1614 f1 = n_singles 1615 s2 = f1 1616 f2 = s2 + n_doubles 1617 1618 d_ab = e_vir[:,None] + e_vir 1619 d_i = e_occ[:,None] 1620 D_n = -d_i + d_ab.reshape(-1) 1621 D_iab = D_n.reshape(-1) 1622 1623 if M_ab is None: 1624 M_ab = adc.get_imds() 1625 1626 #Calculate sigma vector 1627 def sigma_(r): 1628 cput0 = (logger.process_clock(), logger.perf_counter()) 1629 log = logger.Logger(adc.stdout, adc.verbose) 1630 1631 s = np.zeros((dim)) 1632 1633 r1 = r[s1:f1] 1634 r2 = r[s2:f2] 1635 1636 r2 = r2.reshape(nocc,nvir,nvir) 1637 1638############ ADC(2) ab block ############################ 1639 1640 s[s1:f1] = lib.einsum('ab,b->a',M_ab,r1) 1641 1642############# ADC(2) a - ibc and ibc - a coupling blocks ######################### 1643 1644 if isinstance(eris.ovvv, type(None)): 1645 chnk_size = radc_ao2mo.calculate_chunk_size(adc) 1646 else: 1647 chnk_size = nocc 1648 a = 0 1649 temp_doubles = np.zeros((nocc,nvir,nvir)) 1650 for p in range(0,nocc,chnk_size): 1651 if getattr(adc, 'with_df', None): 1652 eris_ovvv = dfadc.get_ovvv_df(adc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 1653 else: 1654 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 1655 k = eris_ovvv.shape[0] 1656 1657 s[s1:f1] += 2. * lib.einsum('icab,ibc->a', eris_ovvv, r2[a:a+k], optimize = True) 1658 s[s1:f1] -= lib.einsum('ibac,ibc->a', eris_ovvv, r2[a:a+k], optimize = True) 1659 1660 temp_doubles[a:a+k] += lib.einsum('icab,a->ibc', eris_ovvv, r1, optimize = True) 1661 del eris_ovvv 1662 a += k 1663 1664 s[s2:f2] += temp_doubles.reshape(-1) 1665################ ADC(2) iab - jcd block ############################ 1666 1667 s[s2:f2] += D_iab * r2.reshape(-1) 1668 1669############### ADC(3) iab - jcd block ############################ 1670 1671 if (method == "adc(2)-x" or method == "adc(3)"): 1672 1673 eris_oovv = eris.oovv 1674 eris_ovvo = eris.ovvo 1675 1676 r2 = r2.reshape(nocc, nvir, nvir) 1677 1678 if isinstance(eris.vvvv, np.ndarray): 1679 r_bab_t = r2.reshape(nocc,-1) 1680 eris_vvvv = eris.vvvv 1681 s[s2:f2] += np.dot(r_bab_t,eris_vvvv.T).reshape(-1) 1682 elif isinstance(eris.vvvv, list): 1683 s[s2:f2] += ea_contract_r_vvvv(adc,r2,eris.vvvv) 1684 else: 1685 s[s2:f2] += ea_contract_r_vvvv(adc,r2,eris.Lvv) 1686 1687 s[s2:f2] -= 0.5*lib.einsum('jzyi,jzx->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1688 s[s2:f2] += 0.5*lib.einsum('jzyi,jxz->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1689 1690 s[s2:f2] -= 0.5*lib.einsum('jiyz,jxz->ixy',eris_oovv,r2,optimize = True).reshape(-1) 1691 s[s2:f2] += 0.5*lib.einsum('jzyi,jxz->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1692 s[s2:f2] -= 0.5*lib.einsum('jixz,jzy->ixy',eris_oovv,r2,optimize = True).reshape(-1) 1693 s[s2:f2] -= 0.5*lib.einsum('jixw,jwy->ixy',eris_oovv,r2,optimize = True).reshape(-1) 1694 s[s2:f2] -= 0.5*lib.einsum('jiyw,jxw->ixy',eris_oovv,r2,optimize = True).reshape(-1) 1695 s[s2:f2] += 0.5*lib.einsum('jwyi,jxw->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1696 s[s2:f2] += 0.5*lib.einsum('jwyi,jxw->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1697 s[s2:f2] -= 0.5*lib.einsum('jwyi,jwx->ixy',eris_ovvo,r2,optimize = True).reshape(-1) 1698 1699 #print("Calculating additional terms for adc(3)") 1700 1701 if (method == "adc(3)"): 1702 1703 eris_ovoo = eris.ovoo 1704 1705############### ADC(3) a - ibc block and ibc-a coupling blocks ######################## 1706 1707 t2_1 = adc.t2[0][:] 1708 1709 temp = 0.25 * lib.einsum('lmab,jab->lmj',t2_1,r2) 1710 temp -= 0.25 * lib.einsum('lmab,jba->lmj',t2_1,r2) 1711 temp -= 0.25 * lib.einsum('mlab,jab->lmj',t2_1,r2) 1712 temp += 0.25 * lib.einsum('mlab,jba->lmj',t2_1,r2) 1713 1714 s[s1:f1] += lib.einsum('lmj,lamj->a',temp, eris_ovoo, optimize=True) 1715 s[s1:f1] -= lib.einsum('lmj,malj->a',temp, eris_ovoo, optimize=True) 1716 del temp 1717 1718 temp_1 = -lib.einsum('lmzw,jzw->jlm',t2_1,r2) 1719 s[s1:f1] -= lib.einsum('jlm,lamj->a',temp_1, eris_ovoo, optimize=True) 1720 1721 temp_s_a = lib.einsum('jlwd,jzw->lzd',t2_1,r2,optimize=True) 1722 temp_s_a -= lib.einsum('jlwd,jwz->lzd',t2_1,r2,optimize=True) 1723 temp_s_a -= lib.einsum('ljwd,jzw->lzd',t2_1,r2,optimize=True) 1724 temp_s_a += lib.einsum('ljwd,jwz->lzd',t2_1,r2,optimize=True) 1725 temp_s_a += lib.einsum('ljdw,jzw->lzd',t2_1,r2,optimize=True) 1726 1727 temp_s_a_1 = -lib.einsum('jlzd,jwz->lwd',t2_1,r2,optimize=True) 1728 temp_s_a_1 += lib.einsum('jlzd,jzw->lwd',t2_1,r2,optimize=True) 1729 temp_s_a_1 += lib.einsum('ljzd,jwz->lwd',t2_1,r2,optimize=True) 1730 temp_s_a_1 -= lib.einsum('ljzd,jzw->lwd',t2_1,r2,optimize=True) 1731 temp_s_a_1 += -lib.einsum('ljdz,jwz->lwd',t2_1,r2,optimize=True) 1732 1733 temp_t2_r2_1 = lib.einsum('jlwd,jzw->lzd',t2_1,r2,optimize=True) 1734 temp_t2_r2_1 -= lib.einsum('jlwd,jwz->lzd',t2_1,r2,optimize=True) 1735 temp_t2_r2_1 += lib.einsum('jlwd,jzw->lzd',t2_1,r2,optimize=True) 1736 temp_t2_r2_1 -= lib.einsum('ljwd,jzw->lzd',t2_1,r2,optimize=True) 1737 1738 temp_t2_r2_2 = -lib.einsum('jlzd,jwz->lwd',t2_1,r2,optimize=True) 1739 temp_t2_r2_2 += lib.einsum('jlzd,jzw->lwd',t2_1,r2,optimize=True) 1740 temp_t2_r2_2 -= lib.einsum('jlzd,jwz->lwd',t2_1,r2,optimize=True) 1741 temp_t2_r2_2 += lib.einsum('ljzd,jwz->lwd',t2_1,r2,optimize=True) 1742 1743 temp_t2_r2_3 = -lib.einsum('ljzd,jzw->lwd',t2_1,r2,optimize=True) 1744 1745 temp_a = t2_1.transpose(0,3,1,2).copy() 1746 temp_b = temp_a.reshape(nocc*nvir,nocc*nvir) 1747 r2_t = r2.reshape(nocc*nvir,-1) 1748 temp_c = np.dot(temp_b,r2_t).reshape(nocc,nvir,nvir) 1749 temp_t2_r2_4 = temp_c.transpose(0,2,1).copy() 1750 1751 del t2_1 1752 1753 if isinstance(eris.ovvv, type(None)): 1754 chnk_size = radc_ao2mo.calculate_chunk_size(adc) 1755 else: 1756 chnk_size = nocc 1757 a = 0 1758 temp = np.zeros((nocc,nvir,nvir)) 1759 temp_1_1 = np.zeros((nocc,nvir,nvir)) 1760 temp_2_1 = np.zeros((nocc,nvir,nvir)) 1761 for p in range(0,nocc,chnk_size): 1762 if getattr(adc, 'with_df', None): 1763 eris_ovvv = dfadc.get_ovvv_df(adc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 1764 else: 1765 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 1766 k = eris_ovvv.shape[0] 1767 1768 temp_1_1[a:a+k] = lib.einsum('ldxb,b->lxd', eris_ovvv,r1,optimize=True) 1769 temp_1_1[a:a+k] -= lib.einsum('lbxd,b->lxd', eris_ovvv,r1,optimize=True) 1770 temp_2_1[a:a+k] = lib.einsum('ldxb,b->lxd', eris_ovvv,r1,optimize=True) 1771 1772 s[s1:f1] += 0.5*lib.einsum('lzd,ldza->a',temp_s_a[a:a+k],eris_ovvv,optimize=True) 1773 s[s1:f1] -= 0.5*lib.einsum('lzd,lazd->a',temp_s_a[a:a+k],eris_ovvv,optimize=True) 1774 s[s1:f1] -= 0.5*lib.einsum('lwd,ldwa->a',temp_s_a_1[a:a+k],eris_ovvv,optimize=True) 1775 s[s1:f1] += 0.5*lib.einsum('lwd,lawd->a',temp_s_a_1[a:a+k],eris_ovvv,optimize=True) 1776 1777 s[s1:f1] += 0.5*lib.einsum('lzd,ldza->a',temp_t2_r2_1[a:a+k],eris_ovvv,optimize=True) 1778 1779 s[s1:f1] -= 0.5*lib.einsum('lwd,ldwa->a',temp_t2_r2_2[a:a+k],eris_ovvv,optimize=True) 1780 1781 s[s1:f1] += 0.5*lib.einsum('lwd,lawd->a',temp_t2_r2_3[a:a+k],eris_ovvv,optimize=True) 1782 1783 s[s1:f1] -= 0.5*lib.einsum('lzd,lazd->a',temp_t2_r2_4[a:a+k],eris_ovvv,optimize=True) 1784 1785 temp[a:a+k] -= lib.einsum('lbyd,b->lyd',eris_ovvv,r1,optimize=True) 1786 1787 del eris_ovvv 1788 a += k 1789 1790 t2_1 = adc.t2[0][:] 1791 temp_1 = -lib.einsum('lyd,lixd->ixy',temp,t2_1,optimize=True) 1792 s[s2:f2] -= temp_1.reshape(-1) 1793 1794 del temp_s_a 1795 del temp_s_a_1 1796 del temp_t2_r2_1 1797 del temp_t2_r2_2 1798 del temp_t2_r2_3 1799 del temp_t2_r2_4 1800 1801 temp_1 = lib.einsum('b,lbmi->lmi',r1,eris_ovoo) 1802 s[s2:f2] += lib.einsum('lmi,lmxy->ixy',temp_1, t2_1, optimize=True).reshape(-1) 1803 1804 temp = lib.einsum('lxd,lidy->ixy',temp_1_1,t2_1,optimize=True) 1805 temp += lib.einsum('lxd,ilyd->ixy',temp_2_1,t2_1,optimize=True) 1806 temp -= lib.einsum('lxd,ildy->ixy',temp_2_1,t2_1,optimize=True) 1807 s[s2:f2] += temp.reshape(-1) 1808 1809 del t2_1 1810 del temp 1811 del temp_1 1812 del temp_1_1 1813 del temp_2_1 1814 1815 cput0 = log.timer_debug1("completed sigma vector calculation", *cput0) 1816 return s 1817 1818 return sigma_ 1819 1820 1821def ip_adc_matvec(adc, M_ij=None, eris=None): 1822 1823 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 1824 raise NotImplementedError(adc.method) 1825 1826 method = adc.method 1827 1828 1829 nocc = adc._nocc 1830 nvir = adc._nvir 1831 1832 n_singles = nocc 1833 n_doubles = nvir * nocc * nocc 1834 1835 dim = n_singles + n_doubles 1836 1837 e_occ = adc.mo_energy[:nocc] 1838 e_vir = adc.mo_energy[nocc:] 1839 1840 if eris is None: 1841 eris = adc.transform_integrals() 1842 1843 s1 = 0 1844 f1 = n_singles 1845 s2 = f1 1846 f2 = s2 + n_doubles 1847 1848 d_ij = e_occ[:,None] + e_occ 1849 d_a = e_vir[:,None] 1850 D_n = -d_a + d_ij.reshape(-1) 1851 D_aij = D_n.reshape(-1) 1852 1853 if M_ij is None: 1854 M_ij = adc.get_imds() 1855 1856 #Calculate sigma vector 1857 def sigma_(r): 1858 cput0 = (logger.process_clock(), logger.perf_counter()) 1859 log = logger.Logger(adc.stdout, adc.verbose) 1860 1861 s = np.zeros((dim)) 1862 1863 r1 = r[s1:f1] 1864 r2 = r[s2:f2] 1865 1866 r2 = r2.reshape(nvir,nocc,nocc) 1867 1868 eris_ovoo = eris.ovoo 1869 1870############ ADC(2) ij block ############################ 1871 1872 s[s1:f1] = lib.einsum('ij,j->i',M_ij,r1) 1873 1874############ ADC(2) i - kja block ######################### 1875 1876 s[s1:f1] += 2. * lib.einsum('jaki,ajk->i', eris_ovoo, r2, optimize = True) 1877 s[s1:f1] -= lib.einsum('kaji,ajk->i', eris_ovoo, r2, optimize = True) 1878 1879############## ADC(2) ajk - i block ############################ 1880 1881 temp = lib.einsum('jaki,i->ajk', eris_ovoo, r1, optimize = True).reshape(-1) 1882 s[s2:f2] += temp.reshape(-1) 1883 1884################ ADC(2) ajk - bil block ############################ 1885 1886 s[s2:f2] += D_aij * r2.reshape(-1) 1887 1888############### ADC(3) ajk - bil block ############################ 1889 1890 if (method == "adc(2)-x" or method == "adc(3)"): 1891 1892 eris_oooo = eris.oooo 1893 eris_oovv = eris.oovv 1894 eris_ovvo = eris.ovvo 1895 1896 s[s2:f2] -= 0.5*lib.einsum('kijl,ali->ajk',eris_oooo, r2, optimize = True).reshape(-1) 1897 s[s2:f2] -= 0.5*lib.einsum('klji,ail->ajk',eris_oooo ,r2, optimize = True).reshape(-1) 1898 1899 s[s2:f2] += 0.5*lib.einsum('klba,bjl->ajk',eris_oovv,r2,optimize = True).reshape(-1) 1900 1901 s[s2:f2] += 0.5*lib.einsum('jabl,bkl->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1902 s[s2:f2] -= 0.5*lib.einsum('jabl,blk->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1903 s[s2:f2] += 0.5*lib.einsum('jlba,blk->ajk',eris_oovv,r2,optimize = True).reshape(-1) 1904 s[s2:f2] -= 0.5*lib.einsum('jabl,blk->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1905 1906 s[s2:f2] += 0.5*lib.einsum('kiba,bji->ajk',eris_oovv,r2,optimize = True).reshape(-1) 1907 1908 s[s2:f2] += 0.5*lib.einsum('jiba,bik->ajk',eris_oovv,r2,optimize = True).reshape(-1) 1909 s[s2:f2] -= 0.5*lib.einsum('jabi,bik->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1910 s[s2:f2] -= 0.5*lib.einsum('jabi,bik->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1911 s[s2:f2] += 0.5*lib.einsum('jabi,bki->ajk',eris_ovvo,r2,optimize = True).reshape(-1) 1912 1913 if (method == "adc(3)"): 1914 1915 eris_ovoo = eris.ovoo 1916 t2_1 = adc.t2[0] 1917 1918################ ADC(3) i - kja block and ajk - i ############################ 1919 1920 temp = 0.25 * lib.einsum('ijbc,aij->abc',t2_1, r2, optimize=True) 1921 temp -= 0.25 * lib.einsum('ijbc,aji->abc',t2_1, r2, optimize=True) 1922 temp -= 0.25 * lib.einsum('jibc,aij->abc',t2_1, r2, optimize=True) 1923 temp += 0.25 * lib.einsum('jibc,aji->abc',t2_1, r2, optimize=True) 1924 1925 temp_1 = lib.einsum('kjcb,ajk->abc',t2_1,r2, optimize=True) 1926 1927 if isinstance(eris.ovvv, type(None)): 1928 chnk_size = radc_ao2mo.calculate_chunk_size(adc) 1929 else: 1930 chnk_size = nocc 1931 a = 0 1932 temp_singles = np.zeros((nocc)) 1933 temp_doubles = np.zeros((nvir,nvir,nvir)) 1934 for p in range(0,nocc,chnk_size): 1935 if getattr(adc, 'with_df', None): 1936 eris_ovvv = dfadc.get_ovvv_df(adc, eris.Lov, eris.Lvv, p, chnk_size).reshape(-1,nvir,nvir,nvir) 1937 else: 1938 eris_ovvv = radc_ao2mo.unpack_eri_1(eris.ovvv, nvir) 1939 k = eris_ovvv.shape[0] 1940 1941 temp_singles[a:a+k] += lib.einsum('abc,icab->i',temp, eris_ovvv, optimize=True) 1942 temp_singles[a:a+k] -= lib.einsum('abc,ibac->i',temp, eris_ovvv, optimize=True) 1943 temp_singles[a:a+k] += lib.einsum('abc,icab->i',temp_1, eris_ovvv, optimize=True) 1944 temp_doubles = lib.einsum('i,icab->cba',r1[a:a+k],eris_ovvv,optimize=True) 1945 s[s2:f2] += lib.einsum('cba,kjcb->ajk',temp_doubles, t2_1, optimize=True).reshape(-1) 1946 del eris_ovvv 1947 del temp_doubles 1948 a += k 1949 1950 s[s1:f1] += temp_singles 1951 temp = np.zeros_like(r2) 1952 temp = lib.einsum('jlab,ajk->blk',t2_1,r2,optimize=True) 1953 temp -= lib.einsum('jlab,akj->blk',t2_1,r2,optimize=True) 1954 temp -= lib.einsum('ljab,ajk->blk',t2_1,r2,optimize=True) 1955 temp += lib.einsum('ljab,akj->blk',t2_1,r2,optimize=True) 1956 temp += lib.einsum('ljba,ajk->blk',t2_1,r2,optimize=True) 1957 1958 temp_1 = np.zeros_like(r2) 1959 temp_1 = lib.einsum('jlab,ajk->blk',t2_1,r2,optimize=True) 1960 temp_1 -= lib.einsum('jlab,akj->blk',t2_1,r2,optimize=True) 1961 temp_1 += lib.einsum('jlab,ajk->blk',t2_1,r2,optimize=True) 1962 temp_1 -= lib.einsum('ljab,ajk->blk',t2_1,r2,optimize=True) 1963 1964 temp_2 = lib.einsum('jlba,akj->blk',t2_1,r2, optimize=True) 1965 1966 s[s1:f1] += 0.5*lib.einsum('blk,lbik->i',temp,eris_ovoo,optimize=True) 1967 s[s1:f1] -= 0.5*lib.einsum('blk,iblk->i',temp,eris_ovoo,optimize=True) 1968 s[s1:f1] += 0.5*lib.einsum('blk,lbik->i',temp_1,eris_ovoo,optimize=True) 1969 s[s1:f1] -= 0.5*lib.einsum('blk,iblk->i',temp_2,eris_ovoo,optimize=True) 1970 del temp 1971 del temp_1 1972 del temp_2 1973 1974 temp = np.zeros_like(r2) 1975 temp = -lib.einsum('klab,akj->blj',t2_1,r2,optimize=True) 1976 temp += lib.einsum('klab,ajk->blj',t2_1,r2,optimize=True) 1977 temp += lib.einsum('lkab,akj->blj',t2_1,r2,optimize=True) 1978 temp -= lib.einsum('lkab,ajk->blj',t2_1,r2,optimize=True) 1979 temp -= lib.einsum('lkba,akj->blj',t2_1,r2,optimize=True) 1980 1981 temp_1 = np.zeros_like(r2) 1982 temp_1 = -lib.einsum('klab,akj->blj',t2_1,r2,optimize=True) 1983 temp_1 += lib.einsum('klab,ajk->blj',t2_1,r2,optimize=True) 1984 temp_1 -= lib.einsum('klab,akj->blj',t2_1,r2,optimize=True) 1985 temp_1 += lib.einsum('lkab,akj->blj',t2_1,r2,optimize=True) 1986 1987 temp_2 = -lib.einsum('klba,ajk->blj',t2_1,r2,optimize=True) 1988 1989 s[s1:f1] -= 0.5*lib.einsum('blj,lbij->i',temp,eris_ovoo,optimize=True) 1990 s[s1:f1] += 0.5*lib.einsum('blj,iblj->i',temp,eris_ovoo,optimize=True) 1991 s[s1:f1] -= 0.5*lib.einsum('blj,lbij->i',temp_1,eris_ovoo,optimize=True) 1992 s[s1:f1] += 0.5*lib.einsum('blj,iblj->i',temp_2,eris_ovoo,optimize=True) 1993 1994 del temp 1995 del temp_1 1996 del temp_2 1997 1998 temp_1 = lib.einsum('i,lbik->kbl',r1,eris_ovoo) 1999 temp_1 -= lib.einsum('i,iblk->kbl',r1,eris_ovoo) 2000 temp_2 = lib.einsum('i,lbik->kbl',r1,eris_ovoo) 2001 2002 temp = lib.einsum('kbl,ljba->ajk',temp_1,t2_1,optimize=True) 2003 temp += lib.einsum('kbl,jlab->ajk',temp_2,t2_1,optimize=True) 2004 temp -= lib.einsum('kbl,ljab->ajk',temp_2,t2_1,optimize=True) 2005 s[s2:f2] += temp.reshape(-1) 2006 2007 temp = -lib.einsum('i,iblj->jbl',r1,eris_ovoo,optimize=True) 2008 temp_1 = -lib.einsum('jbl,klba->ajk',temp,t2_1,optimize=True) 2009 s[s2:f2] -= temp_1.reshape(-1) 2010 2011 del temp 2012 del temp_1 2013 del temp_2 2014 del t2_1 2015 2016 cput0 = log.timer_debug1("completed sigma vector calculation", *cput0) 2017 s *= -1.0 2018 2019 return s 2020 2021 return sigma_ 2022 2023 2024def ea_compute_trans_moments(adc, orb): 2025 2026 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 2027 raise NotImplementedError(adc.method) 2028 2029 method = adc.method 2030 2031 t2_1 = adc.t2[0][:] 2032 t1_2 = adc.t1[0][:] 2033 2034 nocc = adc._nocc 2035 nvir = adc._nvir 2036 2037 n_singles = nvir 2038 n_doubles = nocc * nvir * nvir 2039 2040 dim = n_singles + n_doubles 2041 2042 idn_vir = np.identity(nvir) 2043 2044 s1 = 0 2045 f1 = n_singles 2046 s2 = f1 2047 f2 = s2 + n_doubles 2048 2049 T = np.zeros((dim)) 2050 2051######## ADC(2) part ############################################ 2052 2053 if orb < nocc: 2054 2055 T[s1:f1] = -t1_2[orb,:] 2056 2057 t2_1_t = -t2_1.transpose(1,0,2,3) 2058 2059 T[s2:f2] += t2_1_t[:,orb,:,:].reshape(-1) 2060 2061 else: 2062 2063 T[s1:f1] += idn_vir[(orb-nocc), :] 2064 T[s1:f1] -= 0.25*lib.einsum('klc,klac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2065 T[s1:f1] -= 0.25*lib.einsum('lkc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2066 2067 T[s1:f1] -= 0.25*lib.einsum('klc,klac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2068 T[s1:f1] += 0.25*lib.einsum('lkc,klac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2069 T[s1:f1] += 0.25*lib.einsum('klc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2070 T[s1:f1] -= 0.25*lib.einsum('lkc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_1, optimize = True) 2071 2072######### ADC(3) 2p-1h part ############################################ 2073 2074 if(method=="adc(2)-x"or adc.method=="adc(3)"): 2075 2076 t2_2 = adc.t2[1][:] 2077 2078 if orb < nocc: 2079 2080 t2_2_t = -t2_2.transpose(1,0,2,3) 2081 2082 T[s2:f2] += t2_2_t[:,orb,:,:].reshape(-1) 2083 2084########### ADC(3) 1p part ############################################ 2085 2086 if(adc.method=="adc(3)"): 2087 2088 t1_3 = adc.t1[1] 2089 2090 if orb < nocc: 2091 T[s1:f1] += 0.5*lib.einsum('kac,ck->a',t2_1[:,orb,:,:], t1_2.T,optimize = True) 2092 T[s1:f1] -= 0.5*lib.einsum('kac,ck->a',t2_1[orb,:,:,:], t1_2.T,optimize = True) 2093 T[s1:f1] -= 0.5*lib.einsum('kac,ck->a',t2_1[orb,:,:,:], t1_2.T,optimize = True) 2094 T[s1:f1] -= t1_3[orb,:] 2095 2096 else: 2097 2098 T[s1:f1] -= 0.25*lib.einsum('klc,klac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2099 T[s1:f1] -= 0.25*lib.einsum('lkc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2100 2101 T[s1:f1] -= 0.25*lib.einsum('klac,klc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2102 T[s1:f1] -= 0.25*lib.einsum('lkac,lkc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2103 2104 T[s1:f1] -= 0.25*lib.einsum('klc,klac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2105 T[s1:f1] += 0.25*lib.einsum('klc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2106 T[s1:f1] += 0.25*lib.einsum('lkc,klac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2107 T[s1:f1] -= 0.25*lib.einsum('lkc,lkac->a',t2_1[:,:,(orb-nocc),:], t2_2, optimize = True) 2108 2109 T[s1:f1] -= 0.25*lib.einsum('klac,klc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2110 T[s1:f1] += 0.25*lib.einsum('klac,lkc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2111 T[s1:f1] += 0.25*lib.einsum('lkac,klc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2112 T[s1:f1] -= 0.25*lib.einsum('lkac,lkc->a',t2_1, t2_2[:,:,(orb-nocc),:],optimize = True) 2113 2114 del t2_2 2115 del t2_1 2116 2117 T_aaa = T[n_singles:].reshape(nocc,nvir,nvir).copy() 2118 T_aaa = T_aaa - T_aaa.transpose(0,2,1) 2119 T[n_singles:] += T_aaa.reshape(-1) 2120 2121 return T 2122 2123 2124def ip_compute_trans_moments(adc, orb): 2125 2126 if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): 2127 raise NotImplementedError(adc.method) 2128 2129 method = adc.method 2130 2131 t2_1 = adc.t2[0][:] 2132 t1_2 = adc.t1[0][:] 2133 2134 nocc = adc._nocc 2135 nvir = adc._nvir 2136 2137 n_singles = nocc 2138 n_doubles = nvir * nocc * nocc 2139 2140 dim = n_singles + n_doubles 2141 2142 idn_occ = np.identity(nocc) 2143 2144 s1 = 0 2145 f1 = n_singles 2146 s2 = f1 2147 f2 = s2 + n_doubles 2148 2149 T = np.zeros((dim)) 2150 2151######## ADC(2) 1h part ############################################ 2152 2153 if orb < nocc: 2154 T[s1:f1] = idn_occ[orb, :] 2155 T[s1:f1] += 0.25*lib.einsum('kdc,ikdc->i',t2_1[:,orb,:,:], t2_1, optimize = True) 2156 T[s1:f1] -= 0.25*lib.einsum('kcd,ikdc->i',t2_1[:,orb,:,:], t2_1, optimize = True) 2157 T[s1:f1] -= 0.25*lib.einsum('kdc,ikcd->i',t2_1[:,orb,:,:], t2_1, optimize = True) 2158 T[s1:f1] += 0.25*lib.einsum('kcd,ikcd->i',t2_1[:,orb,:,:], t2_1, optimize = True) 2159 T[s1:f1] -= 0.25*lib.einsum('kdc,ikdc->i',t2_1[orb,:,:,:], t2_1, optimize = True) 2160 T[s1:f1] -= 0.25*lib.einsum('kcd,ikcd->i',t2_1[orb,:,:,:], t2_1, optimize = True) 2161 else: 2162 T[s1:f1] += t1_2[:,(orb-nocc)] 2163 2164######## ADC(2) 2h-1p part ############################################ 2165 2166 t2_1_t = t2_1.transpose(2,3,1,0) 2167 2168 T[s2:f2] = t2_1_t[(orb-nocc),:,:,:].reshape(-1) 2169 2170######## ADC(3) 2h-1p part ############################################ 2171 2172 if(method=='adc(2)-x'or method=='adc(3)'): 2173 2174 t2_2 = adc.t2[1][:] 2175 2176 if orb >= nocc: 2177 t2_2_t = t2_2.transpose(2,3,1,0) 2178 2179 T[s2:f2] += t2_2_t[(orb-nocc),:,:,:].reshape(-1) 2180 2181######### ADC(3) 1h part ############################################ 2182 2183 if(method=='adc(3)'): 2184 2185 t1_3 = adc.t1[1] 2186 2187 if orb < nocc: 2188 T[s1:f1] += 0.25*lib.einsum('kdc,ikdc->i',t2_1[:,orb,:,:], t2_2, optimize = True) 2189 T[s1:f1] -= 0.25*lib.einsum('kcd,ikdc->i',t2_1[:,orb,:,:], t2_2, optimize = True) 2190 T[s1:f1] -= 0.25*lib.einsum('kdc,ikcd->i',t2_1[:,orb,:,:], t2_2, optimize = True) 2191 T[s1:f1] += 0.25*lib.einsum('kcd,ikcd->i',t2_1[:,orb,:,:], t2_2, optimize = True) 2192 T[s1:f1] -= 0.25*lib.einsum('kdc,ikdc->i',t2_1[orb,:,:,:], t2_2, optimize = True) 2193 T[s1:f1] -= 0.25*lib.einsum('kcd,ikcd->i',t2_1[orb,:,:,:], t2_2, optimize = True) 2194 2195 T[s1:f1] += 0.25*lib.einsum('ikdc,kdc->i',t2_1, t2_2[:,orb,:,:],optimize = True) 2196 T[s1:f1] -= 0.25*lib.einsum('ikcd,kdc->i',t2_1, t2_2[:,orb,:,:],optimize = True) 2197 T[s1:f1] -= 0.25*lib.einsum('ikdc,kcd->i',t2_1, t2_2[:,orb,:,:],optimize = True) 2198 T[s1:f1] += 0.25*lib.einsum('ikcd,kcd->i',t2_1, t2_2[:,orb,:,:],optimize = True) 2199 T[s1:f1] -= 0.25*lib.einsum('ikcd,kcd->i',t2_1, t2_2[orb,:,:,:],optimize = True) 2200 T[s1:f1] -= 0.25*lib.einsum('ikdc,kdc->i',t2_1, t2_2[orb,:,:,:],optimize = True) 2201 else: 2202 T[s1:f1] += 0.5*lib.einsum('ikc,kc->i',t2_1[:,:,(orb-nocc),:], t1_2,optimize = True) 2203 T[s1:f1] -= 0.5*lib.einsum('kic,kc->i',t2_1[:,:,(orb-nocc),:], t1_2,optimize = True) 2204 T[s1:f1] += 0.5*lib.einsum('ikc,kc->i',t2_1[:,:,(orb-nocc),:], t1_2,optimize = True) 2205 T[s1:f1] += t1_3[:,(orb-nocc)] 2206 2207 del t2_2 2208 del t2_1 2209 2210 T_aaa = T[n_singles:].reshape(nvir,nocc,nocc).copy() 2211 T_aaa = T_aaa - T_aaa.transpose(0,2,1) 2212 T[n_singles:] += T_aaa.reshape(-1) 2213 2214 return T 2215 2216 2217def get_trans_moments(adc): 2218 2219 nmo = adc.nmo 2220 T = [] 2221 for orb in range(nmo): 2222 2223 T_a = adc.compute_trans_moments(orb) 2224 T.append(T_a) 2225 2226 T = np.array(T) 2227 return T 2228 2229 2230def analyze_eigenvector_ea(adc): 2231 2232 nocc = adc._nocc 2233 nvir = adc._nvir 2234 evec_print_tol = adc.evec_print_tol 2235 2236 logger.info(adc, "Number of occupied orbitals = %d", nocc) 2237 logger.info(adc, "Number of virtual orbitals = %d", nvir) 2238 logger.info(adc, "Print eigenvector elements > %f\n", evec_print_tol) 2239 2240 n_singles = nvir 2241 U = adc.U 2242 2243 for I in range(U.shape[1]): 2244 U1 = U[:n_singles,I] 2245 U2 = U[n_singles:,I].reshape(nocc,nvir,nvir) 2246 U1dotU1 = np.dot(U1, U1) 2247 U2dotU2 = 2.*np.dot(U2.ravel(), U2.ravel()) - np.dot(U2.ravel(), U2.transpose(0,2,1).ravel()) 2248 2249 U_sq = U[:,I].copy()**2 2250 ind_idx = np.argsort(-U_sq) 2251 U_sq = U_sq[ind_idx] 2252 U_sorted = U[ind_idx,I].copy() 2253 2254 U_sorted = U_sorted[U_sq > evec_print_tol**2] 2255 ind_idx = ind_idx[U_sq > evec_print_tol**2] 2256 2257 singles_idx = [] 2258 doubles_idx = [] 2259 singles_val = [] 2260 doubles_val = [] 2261 iter_num = 0 2262 2263 for orb_idx in ind_idx: 2264 2265 if orb_idx < n_singles: 2266 a_idx = orb_idx + 1 + nocc 2267 singles_idx.append(a_idx) 2268 singles_val.append(U_sorted[iter_num]) 2269 2270 if orb_idx >= n_singles: 2271 iab_idx = orb_idx - n_singles 2272 ab_rem = iab_idx % (nvir*nvir) 2273 i_idx = iab_idx //(nvir*nvir) 2274 a_idx = ab_rem//nvir 2275 b_idx = ab_rem % nvir 2276 doubles_idx.append((i_idx + 1, a_idx + 1 + nocc, b_idx + 1 + nocc)) 2277 doubles_val.append(U_sorted[iter_num]) 2278 2279 iter_num += 1 2280 2281 logger.info(adc, '%s | root %d | norm(1p) = %6.4f | norm(1h2p) = %6.4f ', 2282 adc.method ,I, U1dotU1, U2dotU2) 2283 2284 if singles_val: 2285 logger.info(adc, "\n1p block: ") 2286 logger.info(adc, " a U(a)") 2287 logger.info(adc, "------------------") 2288 for idx, print_singles in enumerate(singles_idx): 2289 logger.info(adc, ' %4d %7.4f', print_singles, singles_val[idx]) 2290 2291 if doubles_val: 2292 logger.info(adc, "\n1h2p block: ") 2293 logger.info(adc, " i a b U(i,a,b)") 2294 logger.info(adc, "-------------------------------") 2295 for idx, print_doubles in enumerate(doubles_idx): 2296 logger.info(adc, ' %4d %4d %4d %7.4f', 2297 print_doubles[0], print_doubles[1], print_doubles[2], doubles_val[idx]) 2298 2299 logger.info(adc, "\n*************************************************************\n") 2300 2301 2302def analyze_eigenvector_ip(adc): 2303 2304 nocc = adc._nocc 2305 nvir = adc._nvir 2306 2307 n_singles = nocc 2308 evec_print_tol = adc.evec_print_tol 2309 U = adc.U 2310 2311 logger.info(adc, "Number of occupied orbitals = %d", nocc) 2312 logger.info(adc, "Number of virtual orbitals = %d", nvir) 2313 logger.info(adc, "Print eigenvector elements > %f\n", evec_print_tol) 2314 2315 for I in range(U.shape[1]): 2316 U1 = U[:n_singles,I] 2317 U2 = U[n_singles:,I].reshape(nvir,nocc,nocc) 2318 U1dotU1 = np.dot(U1, U1) 2319 U2dotU2 = 2.*np.dot(U2.ravel(), U2.ravel()) - np.dot(U2.ravel(), U2.transpose(0,2,1).ravel()) 2320 2321 U_sq = U[:,I].copy()**2 2322 ind_idx = np.argsort(-U_sq) 2323 U_sq = U_sq[ind_idx] 2324 U_sorted = U[ind_idx,I].copy() 2325 2326 U_sorted = U_sorted[U_sq > evec_print_tol**2] 2327 ind_idx = ind_idx[U_sq > evec_print_tol**2] 2328 2329 singles_idx = [] 2330 doubles_idx = [] 2331 singles_val = [] 2332 doubles_val = [] 2333 iter_num = 0 2334 2335 for orb_idx in ind_idx: 2336 2337 if orb_idx < n_singles: 2338 i_idx = orb_idx + 1 2339 singles_idx.append(i_idx) 2340 singles_val.append(U_sorted[iter_num]) 2341 2342 if orb_idx >= n_singles: 2343 aij_idx = orb_idx - n_singles 2344 ij_rem = aij_idx % (nocc*nocc) 2345 a_idx = aij_idx//(nocc*nocc) 2346 i_idx = ij_rem//nocc 2347 j_idx = ij_rem % nocc 2348 doubles_idx.append((a_idx + 1 + n_singles, i_idx + 1, j_idx + 1)) 2349 doubles_val.append(U_sorted[iter_num]) 2350 2351 iter_num += 1 2352 2353 logger.info(adc, '%s | root %d | norm(1h) = %6.4f | norm(2h1p) = %6.4f ', 2354 adc.method ,I, U1dotU1, U2dotU2) 2355 2356 if singles_val: 2357 logger.info(adc, "\n1h block: ") 2358 logger.info(adc, " i U(i)") 2359 logger.info(adc, "------------------") 2360 for idx, print_singles in enumerate(singles_idx): 2361 logger.info(adc, ' %4d %7.4f', print_singles, singles_val[idx]) 2362 2363 if doubles_val: 2364 logger.info(adc, "\n2h1p block: ") 2365 logger.info(adc, " i j a U(i,j,a)") 2366 logger.info(adc, "-------------------------------") 2367 for idx, print_doubles in enumerate(doubles_idx): 2368 logger.info(adc, ' %4d %4d %4d %7.4f', 2369 print_doubles[1], print_doubles[2], print_doubles[0], doubles_val[idx]) 2370 2371 logger.info(adc, "\n*************************************************************\n") 2372 2373 2374def analyze_spec_factor(adc): 2375 2376 X = adc.X 2377 X_2 = (X.copy()**2)*2 2378 thresh = adc.spec_factor_print_tol 2379 2380 logger.info(adc, "Print spectroscopic factors > %E\n", adc.spec_factor_print_tol) 2381 2382 for i in range(X_2.shape[1]): 2383 2384 sort = np.argsort(-X_2[:,i]) 2385 X_2_row = X_2[:,i] 2386 X_2_row = X_2_row[sort] 2387 2388 if not adc.mol.symmetry: 2389 sym = np.repeat(['A'], X_2_row.shape[0]) 2390 else: 2391 sym = [symm.irrep_id2name(adc.mol.groupname, x) for x in adc._scf.mo_coeff.orbsym] 2392 sym = np.array(sym) 2393 2394 sym = sym[sort] 2395 2396 spec_Contribution = X_2_row[X_2_row > thresh] 2397 index_mo = sort[X_2_row > thresh]+1 2398 2399 if np.sum(spec_Contribution) == 0.0: 2400 continue 2401 2402 logger.info(adc,'%s | root %d \n',adc.method ,i) 2403 logger.info(adc, " HF MO Spec. Contribution Orbital symmetry") 2404 logger.info(adc, "-----------------------------------------------------------") 2405 2406 for c in range(index_mo.shape[0]): 2407 logger.info(adc, ' %3.d %10.8f %s', index_mo[c], spec_Contribution[c], sym[c]) 2408 2409 logger.info(adc, '\nPartial spec. factor sum = %10.8f', np.sum(spec_Contribution)) 2410 logger.info(adc, "\n*************************************************************\n") 2411 2412 2413def renormalize_eigenvectors_ea(adc, nroots=1): 2414 2415 nocc = adc._nocc 2416 nvir = adc._nvir 2417 2418 n_singles = nvir 2419 2420 U = adc.U 2421 2422 for I in range(U.shape[1]): 2423 U1 = U[:n_singles,I] 2424 U2 = U[n_singles:,I].reshape(nocc,nvir,nvir) 2425 UdotU = np.dot(U1, U1) + 2.*np.dot(U2.ravel(), U2.ravel()) - np.dot(U2.ravel(), U2.transpose(0,2,1).ravel()) 2426 U[:,I] /= np.sqrt(UdotU) 2427 2428 return U 2429 2430 2431def renormalize_eigenvectors_ip(adc, nroots=1): 2432 2433 nocc = adc._nocc 2434 nvir = adc._nvir 2435 2436 n_singles = nocc 2437 2438 U = adc.U 2439 2440 for I in range(U.shape[1]): 2441 U1 = U[:n_singles,I] 2442 U2 = U[n_singles:,I].reshape(nvir,nocc,nocc) 2443 UdotU = np.dot(U1, U1) + 2.*np.dot(U2.ravel(), U2.ravel()) - np.dot(U2.ravel(), U2.transpose(0,2,1).ravel()) 2444 U[:,I] /= np.sqrt(UdotU) 2445 2446 return U 2447 2448 2449def get_properties(adc, nroots=1): 2450 2451 #Transition moments 2452 T = adc.get_trans_moments() 2453 2454 #Spectroscopic amplitudes 2455 U = adc.renormalize_eigenvectors(nroots) 2456 X = np.dot(T, U).reshape(-1, nroots) 2457 2458 #Spectroscopic factors 2459 P = 2.0*lib.einsum("pi,pi->i", X, X) 2460 2461 return P,X 2462 2463 2464class RADCEA(RADC): 2465 '''restricted ADC for EA energies and spectroscopic amplitudes 2466 2467 Attributes: 2468 verbose : int 2469 Print level. Default value equals to :class:`Mole.verbose` 2470 max_memory : float or int 2471 Allowed memory in MB. Default value equals to :class:`Mole.max_memory` 2472 incore_complete : bool 2473 Avoid all I/O. Default is False. 2474 method : string 2475 nth-order ADC method. Options are : ADC(2), ADC(2)-X, ADC(3). Default is ADC(2). 2476 conv_tol : float 2477 Convergence threshold for Davidson iterations. Default is 1e-12. 2478 max_cycle : int 2479 Number of Davidson iterations. Default is 50. 2480 max_space : int 2481 Space size to hold trial vectors for Davidson iterative 2482 diagonalization. Default is 12. 2483 2484 Kwargs: 2485 nroots : int 2486 Number of roots (eigenvalues) requested. Default value is 1. 2487 2488 >>> myadc = adc.RADC(mf).run() 2489 >>> myadcea = adc.RADC(myadc).run() 2490 2491 Saved results 2492 2493 e_ea : float or list of floats 2494 EA energy (eigenvalue). For nroots = 1, it is a single float 2495 number. If nroots > 1, it is a list of floats for the lowest 2496 nroots eigenvalues. 2497 v_ip : array 2498 Eigenvectors for each EA transition. 2499 p_ea : float 2500 Spectroscopic amplitudes for each EA transition. 2501 ''' 2502 def __init__(self, adc): 2503 self.mol = adc.mol 2504 self.verbose = adc.verbose 2505 self.stdout = adc.stdout 2506 self.max_memory = adc.max_memory 2507 self.max_space = adc.max_space 2508 self.max_cycle = adc.max_cycle 2509 self.conv_tol = adc.conv_tol 2510 self.tol_residual = adc.tol_residual 2511 self.t1 = adc.t1 2512 self.t2 = adc.t2 2513 self.imds = adc.imds 2514 self.e_corr = adc.e_corr 2515 self.method = adc.method 2516 self.method_type = adc.method_type 2517 self._scf = adc._scf 2518 self._nocc = adc._nocc 2519 self._nvir = adc._nvir 2520 self._nmo = adc._nmo 2521 self.mo_coeff = adc.mo_coeff 2522 self.mo_energy = adc.mo_energy 2523 self.nmo = adc._nmo 2524 self.transform_integrals = adc.transform_integrals 2525 self.with_df = adc.with_df 2526 self.compute_properties = adc.compute_properties 2527 self.E = None 2528 self.U = None 2529 self.P = None 2530 self.X = None 2531 self.evec_print_tol = adc.evec_print_tol 2532 self.spec_factor_print_tol = adc.spec_factor_print_tol 2533 2534 keys = set(('tol_residual','conv_tol', 'e_corr', 'method', 'mo_coeff', 2535 'mo_energy', 'max_memory', 't1', 'max_space', 't2', 2536 'max_cycle')) 2537 2538 self._keys = set(self.__dict__.keys()).union(keys) 2539 2540 kernel = kernel 2541 get_imds = get_imds_ea 2542 matvec = ea_adc_matvec 2543 get_diag = ea_adc_diag 2544 compute_trans_moments = ea_compute_trans_moments 2545 get_trans_moments = get_trans_moments 2546 renormalize_eigenvectors = renormalize_eigenvectors_ea 2547 get_properties = get_properties 2548 analyze_spec_factor = analyze_spec_factor 2549 analyze_eigenvector = analyze_eigenvector_ea 2550 analyze = analyze 2551 compute_dyson_mo = compute_dyson_mo 2552 2553 def get_init_guess(self, nroots=1, diag=None, ascending = True): 2554 if diag is None : 2555 diag = self.ea_adc_diag() 2556 idx = None 2557 if ascending: 2558 idx = np.argsort(diag) 2559 else: 2560 idx = np.argsort(diag)[::-1] 2561 guess = np.zeros((diag.shape[0], nroots)) 2562 min_shape = min(diag.shape[0], nroots) 2563 guess[:min_shape,:min_shape] = np.identity(min_shape) 2564 g = np.zeros((diag.shape[0], nroots)) 2565 g[idx] = guess.copy() 2566 guess = [] 2567 for p in range(g.shape[1]): 2568 guess.append(g[:,p]) 2569 return guess 2570 2571 2572 def gen_matvec(self, imds=None, eris=None): 2573 if imds is None: imds = self.get_imds(eris) 2574 diag = self.get_diag(imds, eris) 2575 matvec = self.matvec(imds, eris) 2576 return matvec, diag 2577 2578 2579class RADCIP(RADC): 2580 '''restricted ADC for IP energies and spectroscopic amplitudes 2581 2582 Attributes: 2583 verbose : int 2584 Print level. Default value equals to :class:`Mole.verbose` 2585 max_memory : float or int 2586 Allowed memory in MB. Default value equals to :class:`Mole.max_memory` 2587 incore_complete : bool 2588 Avoid all I/O. Default is False. 2589 method : string 2590 nth-order ADC method. Options are : ADC(2), ADC(2)-X, ADC(3). Default is ADC(2). 2591 conv_tol : float 2592 Convergence threshold for Davidson iterations. Default is 1e-12. 2593 max_cycle : int 2594 Number of Davidson iterations. Default is 50. 2595 max_space : int 2596 Space size to hold trial vectors for Davidson iterative diagonalization. Default is 12. 2597 2598 Kwargs: 2599 nroots : int 2600 Number of roots (eigenvalues) requested. Default value is 1. 2601 2602 >>> myadc = adc.RADC(mf).run() 2603 >>> myadcip = adc.RADC(myadc).run() 2604 2605 Saved results 2606 2607 e_ip : float or list of floats 2608 IP energy (eigenvalue). For nroots = 1, it is a single float 2609 number. If nroots > 1, it is a list of floats for the lowest 2610 nroots eigenvalues. 2611 v_ip : array 2612 Eigenvectors for each IP transition. 2613 p_ip : float 2614 Spectroscopic amplitudes for each IP transition. 2615 ''' 2616 def __init__(self, adc): 2617 self.mol = adc.mol 2618 self.verbose = adc.verbose 2619 self.stdout = adc.stdout 2620 self.max_memory = adc.max_memory 2621 self.max_space = adc.max_space 2622 self.max_cycle = adc.max_cycle 2623 self.conv_tol = adc.conv_tol 2624 self.tol_residual = adc.tol_residual 2625 self.t1 = adc.t1 2626 self.t2 = adc.t2 2627 self.imds = adc.imds 2628 self.e_corr = adc.e_corr 2629 self.method = adc.method 2630 self.method_type = adc.method_type 2631 self._scf = adc._scf 2632 self._nocc = adc._nocc 2633 self._nvir = adc._nvir 2634 self._nmo = adc._nmo 2635 self.mo_coeff = adc.mo_coeff 2636 self.mo_energy = adc.mo_energy 2637 self.nmo = adc._nmo 2638 self.transform_integrals = adc.transform_integrals 2639 self.with_df = adc.with_df 2640 self.compute_properties = adc.compute_properties 2641 self.E = None 2642 self.U = None 2643 self.P = None 2644 self.X = None 2645 self.evec_print_tol = adc.evec_print_tol 2646 self.spec_factor_print_tol = adc.spec_factor_print_tol 2647 2648 keys = set(('tol_residual','conv_tol', 'e_corr', 'method', 'mo_coeff', 2649 'mo_energy_b', 'max_memory', 't1', 'mo_energy_a', 2650 'max_space', 't2', 'max_cycle')) 2651 self._keys = set(self.__dict__.keys()).union(keys) 2652 2653 kernel = kernel 2654 get_imds = get_imds_ip 2655 get_diag = ip_adc_diag 2656 matvec = ip_adc_matvec 2657 compute_trans_moments = ip_compute_trans_moments 2658 get_trans_moments = get_trans_moments 2659 renormalize_eigenvectors = renormalize_eigenvectors_ip 2660 get_properties = get_properties 2661 analyze_spec_factor = analyze_spec_factor 2662 analyze_eigenvector = analyze_eigenvector_ip 2663 analyze = analyze 2664 compute_dyson_mo = compute_dyson_mo 2665 2666 def get_init_guess(self, nroots=1, diag=None, ascending = True): 2667 if diag is None : 2668 diag = self.ip_adc_diag() 2669 idx = None 2670 if ascending: 2671 idx = np.argsort(diag) 2672 else: 2673 idx = np.argsort(diag)[::-1] 2674 guess = np.zeros((diag.shape[0], nroots)) 2675 min_shape = min(diag.shape[0], nroots) 2676 guess[:min_shape,:min_shape] = np.identity(min_shape) 2677 g = np.zeros((diag.shape[0], nroots)) 2678 g[idx] = guess.copy() 2679 guess = [] 2680 for p in range(g.shape[1]): 2681 guess.append(g[:,p]) 2682 return guess 2683 2684 def gen_matvec(self, imds=None, eris=None): 2685 if imds is None: imds = self.get_imds(eris) 2686 diag = self.get_diag(imds, eris) 2687 matvec = self.matvec(imds, eris) 2688 return matvec, diag 2689 2690if __name__ == '__main__': 2691 from pyscf import scf 2692 from pyscf import gto 2693 from pyscf import adc 2694 2695 r = 1.098 2696 mol = gto.Mole() 2697 mol.atom = [ 2698 ['N', ( 0., 0. , -r/2 )], 2699 ['N', ( 0., 0. , r/2)],] 2700 mol.basis = {'N':'aug-cc-pvdz'} 2701 mol.verbose = 0 2702 mol.build() 2703 mf = scf.RHF(mol) 2704 mf.conv_tol = 1e-12 2705 mf.kernel() 2706 2707 myadc = adc.ADC(mf) 2708 ecorr, t_amp1, t_amp2 = myadc.kernel_gs() 2709 print(ecorr - -0.3220169236051954) 2710 2711 myadcip = RADCIP(myadc) 2712 e,v,p = kernel(myadcip,nroots=3) 2713 print("ADC(2) IP energies") 2714 print (e[0] - 0.5434389910483670) 2715 print (e[1] - 0.6240296243595950) 2716 print (e[2] - 0.6240296243595956) 2717 2718 print("ADC(2) IP spectroscopic factors") 2719 print (p[0] - 1.7688097076459075) 2720 print (p[1] - 1.8192921131700284) 2721 print (p[2] - 1.8192921131700293) 2722 2723 myadcea = RADCEA(myadc) 2724 e,v,p = kernel(myadcea,nroots=3) 2725 print("ADC(2) EA energies") 2726 print (e[0] - 0.0961781923822576) 2727 print (e[1] - 0.1258326916409743) 2728 print (e[2] - 0.1380779405750178) 2729 2730 print("ADC(2) EA spectroscopic factors") 2731 print (p[0] - 1.9832854445007961) 2732 print (p[1] - 1.9634368668786559) 2733 print (p[2] - 1.9783719593912672) 2734 2735 myadc = adc.ADC(mf) 2736 myadc.method = "adc(3)" 2737 ecorr, t_amp1, t_amp2 = myadc.kernel_gs() 2738 print(ecorr - -0.31694173142858517) 2739 2740 myadcip = RADCIP(myadc) 2741 e,v,p = kernel(myadcip,nroots=3) 2742 print("ADC(3) IP energies") 2743 print (e[0] - 0.5667526829981027) 2744 print (e[1] - 0.6099995170092525) 2745 print (e[2] - 0.6099995170092529) 2746 2747 print("ADC(3) IP spectroscopic factors") 2748 print (p[0] - 1.8173191958988848) 2749 print (p[1] - 1.8429224413853840) 2750 print (p[2] - 1.8429224413853851) 2751 2752 myadcea = RADCEA(myadc) 2753 e,v,p = kernel(myadcea,nroots=3) 2754 2755 print("ADC(3) EA energies") 2756 print (e[0] - 0.0936790850738445) 2757 print (e[1] - 0.0983654552141278) 2758 print (e[2] - 0.1295709313652367) 2759 2760 print("ADC(3) EA spectroscopic factors") 2761 print (p[0] - 1.8324175318668088) 2762 print (p[1] - 1.9840991060607487) 2763 print (p[2] - 1.9638550014980212) 2764 2765 myadc.method = "adc(2)-x" 2766 e,v,p = myadc.kernel(nroots=4) 2767 print("ADC(2)-x IP energies") 2768 print (e[0] - 0.5405255360673724) 2769 print (e[1] - 0.6208026698756577) 2770 print (e[2] - 0.6208026698756582) 2771 print (e[3] - 0.6465332771967947) 2772 2773 myadc.method_type = "ea" 2774 e,v,p = myadc.kernel(nroots=4) 2775 print("ADC(2)-x EA energies") 2776 print (e[0] - 0.0953065329985665) 2777 print (e[1] - 0.1238833070823509) 2778 print (e[2] - 0.1365693811939308) 2779 print (e[3] - 0.1365693811939316) 2780