1#!/usr/bin/env python 2# Copyright 2014-2021 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15 16from functools import reduce 17 18import numpy 19import numpy as np 20 21from pyscf import lib 22from pyscf import ao2mo 23from pyscf.lib import logger 24from pyscf.cc import ccsd 25from pyscf.cc import uccsd 26from pyscf.mp import ump2 27from pyscf.cc import gintermediates as imd 28 29#einsum = np.einsum 30einsum = lib.einsum 31 32# This is unrestricted (U)CCSD, i.e. spin-orbital form. 33 34 35def update_amps(cc, t1, t2, eris): 36 nocc, nvir = t1.shape 37 fock = eris.fock 38 39 fov = fock[:nocc,nocc:] 40 foo = fock[:nocc,:nocc] 41 fvv = fock[nocc:,nocc:] 42 43 tau = imd.make_tau(t2, t1, t1) 44 45 Fvv = imd.cc_Fvv(t1, t2, eris) 46 Foo = imd.cc_Foo(t1, t2, eris) 47 Fov = imd.cc_Fov(t1, t2, eris) 48 Woooo = imd.cc_Woooo(t1, t2, eris) 49 Wvvvv = imd.cc_Wvvvv(t1, t2, eris) 50 Wovvo = imd.cc_Wovvo(t1, t2, eris) 51 52 # Move energy terms to the other side 53 Fvv -= np.diag(np.diag(fvv)) 54 Foo -= np.diag(np.diag(foo)) 55 56 # T1 equation 57 t1new = np.array(fov).conj() 58 t1new += einsum('ie,ae->ia', t1, Fvv) 59 t1new += -einsum('ma,mi->ia', t1, Foo) 60 t1new += einsum('imae,me->ia', t2, Fov) 61 t1new += -einsum('nf,naif->ia', t1, eris.ovov) 62 t1new += -0.5*einsum('imef,maef->ia', t2, eris.ovvv) 63 t1new += -0.5*einsum('mnae,mnie->ia', t2, eris.ooov) 64 65 # T2 equation 66 t2new = np.array(eris.oovv).conj() 67 Ftmp = Fvv - 0.5*einsum('mb,me->be', t1, Fov) 68 tmp = einsum('ijae,be->ijab', t2, Ftmp) 69 t2new += (tmp - tmp.transpose(0,1,3,2)) 70 Ftmp = Foo + 0.5*einsum('je,me->mj', t1, Fov) 71 tmp = einsum('imab,mj->ijab', t2, Ftmp) 72 t2new -= (tmp - tmp.transpose(1,0,2,3)) 73 t2new += 0.5*einsum('mnab,mnij->ijab', tau, Woooo) 74 t2new += 0.5*einsum('ijef,abef->ijab', tau, Wvvvv) 75 tmp = einsum('imae,mbej->ijab', t2, Wovvo) 76 tmp -= -einsum('ie,ma,mbje->ijab', t1, t1, eris.ovov) 77 t2new += (tmp - tmp.transpose(0,1,3,2) 78 - tmp.transpose(1,0,2,3) + tmp.transpose(1,0,3,2) ) 79 tmp = einsum('ie,jeba->ijab', t1, np.array(eris.ovvv).conj()) 80 t2new += (tmp - tmp.transpose(1,0,2,3)) 81 tmp = einsum('ma,mbij->ijab', t1, eris.ovoo) 82 t2new -= (tmp - tmp.transpose(0,1,3,2)) 83 84 mo_e = eris.fock.diagonal() 85 eia = mo_e[:nocc,None] - mo_e[None,nocc:] - cc.level_shift 86 eijab = lib.direct_sum('ia,jb->ijab', eia, eia) 87 t1new /= eia 88 t2new /= eijab 89 90 return t1new, t2new 91 92 93def energy(cc, t1, t2, eris): 94 nocc, nvir = t1.shape 95 fock = eris.fock 96 e = einsum('ia,ia', fock[:nocc,nocc:], t1) 97 e += 0.25*np.einsum('ijab,ijab', t2, eris.oovv) 98 e += 0.5 *np.einsum('ia,jb,ijab', t1, t1, eris.oovv) 99 return e.real 100 101 102get_frozen_mask = ump2.get_frozen_mask 103 104 105class UCCSD(ccsd.CCSD): 106 107 def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): 108 ccsd.CCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) 109 # Spin-orbital CCSD needs a stricter tolerance than spatial-orbital 110 self.conv_tol_normt = 1e-6 111 112 @property 113 def nocc(self): 114 nocca, noccb = self.get_nocc() 115 return nocca + noccb 116 117 @property 118 def nmo(self): 119 nmoa, nmob = self.get_nmo() 120 return nmoa + nmob 121 122 get_nocc = uccsd.get_nocc 123 get_nmo = uccsd.get_nmo 124 get_frozen_mask = get_frozen_mask 125 126 def init_amps(self, eris): 127 time0 = logger.process_clock(), logger.perf_counter() 128 mo_e = eris.fock.diagonal() 129 nocc = self.nocc 130 eia = mo_e[:nocc,None] - mo_e[None,nocc:] 131 eijab = lib.direct_sum('ia,jb->ijab',eia,eia) 132 t1 = eris.fock[:nocc,nocc:] / eia 133 eris_oovv = np.array(eris.oovv) 134 t2 = eris_oovv/eijab 135 self.emp2 = 0.25*einsum('ijab,ijab',t2,eris_oovv.conj()).real 136 logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2) 137 logger.timer(self, 'init mp2', *time0) 138 return self.emp2, t1, t2 139 140 energy = energy 141 142 def kernel(self, t1=None, t2=None, eris=None, mbpt2=False): 143 return self.ccsd(t1, t2, eris, mbpt2) 144 def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False): 145 '''Ground-state unrestricted (U)CCSD. 146 147 Kwargs: 148 mbpt2 : bool 149 Use one-shot MBPT2 approximation to CCSD. 150 ''' 151 if mbpt2: 152 #cctyp = 'MBPT2' 153 #self.e_corr, self.t1, self.t2 = self.init_amps(eris) 154 raise NotImplementedError 155 156 if eris is None: 157 eris = self.ao2mo(self.mo_coeff) 158 self.eris = eris 159 return ccsd.CCSD.ccsd(self, t1, t2, eris) 160 161 def ao2mo(self, mo_coeff=None): 162 return _PhysicistsERIs(self, mo_coeff) 163 164 def update_amps(self, t1, t2, eris): 165 return update_amps(self, t1, t2, eris) 166 167 def nip(self): 168 nocc = self.nocc 169 nvir = self.nmo - nocc 170 self._nip = nocc + nocc*(nocc-1)/2*nvir 171 return self._nip 172 173 def nea(self): 174 nocc = self.nocc 175 nvir = self.nmo - nocc 176 self._nea = nvir + nocc*nvir*(nvir-1)/2 177 return self._nea 178 179 def nee(self): 180 nocc = self.nocc 181 nvir = self.nmo - nocc 182 self._nee = nocc*nvir + nocc*(nocc-1)/2*nvir*(nvir-1)/2 183 return self._nee 184 185 def ipccsd_matvec(self, vector): 186 # Ref: Tu, Wang, and Li, J. Chem. Phys. 136, 174102 (2012) Eqs.(8)-(9) 187 if not getattr(self, 'imds', None): 188 self.imds = _IMDS(self) 189 if not self.imds.made_ip_imds: 190 self.imds.make_ip() 191 imds = self.imds 192 193 r1,r2 = self.vector_to_amplitudes_ip(vector) 194 195 # Eq. (8) 196 Hr1 = -einsum('mi,m->i',imds.Foo,r1) 197 Hr1 += einsum('me,mie->i',imds.Fov,r2) 198 Hr1 += -0.5*einsum('nmie,mne->i',imds.Wooov,r2) 199 # Eq. (9) 200 Hr2 = einsum('ae,ije->ija',imds.Fvv,r2) 201 tmp1 = einsum('mi,mja->ija',imds.Foo,r2) 202 Hr2 += (-tmp1 + tmp1.transpose(1,0,2)) 203 Hr2 += -einsum('maji,m->ija',imds.Wovoo,r1) 204 Hr2 += 0.5*einsum('mnij,mna->ija',imds.Woooo,r2) 205 tmp2 = einsum('maei,mje->ija',imds.Wovvo,r2) 206 Hr2 += (tmp2 - tmp2.transpose(1,0,2)) 207 Hr2 += 0.5*einsum('mnef,ijae,mnf->ija',imds.Woovv,self.t2,r2) 208 209 vector = self.amplitudes_to_vector_ip(Hr1,Hr2) 210 return vector 211 212 def ipccsd_diag(self): 213 if not getattr(self, 'imds', None): 214 self.imds = _IMDS(self) 215 if not self.imds.made_ip_imds: 216 self.imds.make_ip() 217 imds = self.imds 218 219 t1, t2 = self.t1, self.t2 220 nocc, nvir = t1.shape 221 222 Hr1 = -np.diag(imds.Foo) 223 Hr2 = np.zeros((nocc,nocc,nvir),dtype=t1.dtype) 224 for i in range(nocc): 225 for j in range(nocc): 226 for a in range(nvir): 227 Hr2[i,j,a] += imds.Fvv[a,a] 228 Hr2[i,j,a] += -imds.Foo[i,i] 229 Hr2[i,j,a] += -imds.Foo[j,j] 230 Hr2[i,j,a] += 0.5*(imds.Woooo[i,j,i,j]-imds.Woooo[j,i,i,j]) 231 Hr2[i,j,a] += imds.Wovvo[i,a,a,i] 232 Hr2[i,j,a] += imds.Wovvo[j,a,a,j] 233 Hr2[i,j,a] += 0.5*(np.dot(imds.Woovv[i,j,:,a],t2[i,j,a,:]) 234 -np.dot(imds.Woovv[j,i,:,a],t2[i,j,a,:])) 235 236 vector = self.amplitudes_to_vector_ip(Hr1,Hr2) 237 return vector 238 239 def vector_to_amplitudes_ip(self,vector): 240 nocc = self.nocc 241 nvir = self.nmo - nocc 242 r1 = vector[:nocc].copy() 243 r2 = np.zeros((nocc,nocc,nvir), vector.dtype) 244 index = nocc 245 for i in range(nocc): 246 for j in range(i): 247 for a in range(nvir): 248 r2[i,j,a] = vector[index] 249 r2[j,i,a] = -vector[index] 250 index += 1 251 return [r1,r2] 252 253 def amplitudes_to_vector_ip(self,r1,r2): 254 nocc = self.nocc 255 nvir = self.nmo - nocc 256 size = nocc + nocc*(nocc-1)/2*nvir 257 vector = np.zeros(size, r1.dtype) 258 vector[:nocc] = r1.copy() 259 index = nocc 260 for i in range(nocc): 261 for j in range(i): 262 for a in range(nvir): 263 vector[index] = r2[i,j,a] 264 index += 1 265 return vector 266 267 def eaccsd_matvec(self,vector): 268 # Ref: Nooijen and Bartlett, J. Chem. Phys. 102, 3629 (1994) Eqs.(30)-(31) 269 if not getattr(self, 'imds', None): 270 self.imds = _IMDS(self) 271 if not self.imds.made_ea_imds: 272 self.imds.make_ea() 273 imds = self.imds 274 275 r1,r2 = self.vector_to_amplitudes_ea(vector) 276 277 # Eq. (30) 278 Hr1 = einsum('ac,c->a',imds.Fvv,r1) 279 Hr1 += einsum('ld,lad->a',imds.Fov,r2) 280 Hr1 += 0.5*einsum('alcd,lcd->a',imds.Wvovv,r2) 281 # Eq. (31) 282 Hr2 = einsum('abcj,c->jab',imds.Wvvvo,r1) 283 tmp1 = einsum('ac,jcb->jab',imds.Fvv,r2) 284 Hr2 += (tmp1 - tmp1.transpose(0,2,1)) 285 Hr2 += -einsum('lj,lab->jab',imds.Foo,r2) 286 tmp2 = einsum('lbdj,lad->jab',imds.Wovvo,r2) 287 Hr2 += (tmp2 - tmp2.transpose(0,2,1)) 288 nvir = self.nmo-self.nocc 289 for a in range(nvir): 290 Hr2[:,a,:] += 0.5*einsum('bcd,jcd->jb',imds.Wvvvv[a],r2) 291 Hr2 += -0.5*einsum('klcd,lcd,kjab->jab',imds.Woovv,r2,self.t2) 292 293 vector = self.amplitudes_to_vector_ea(Hr1,Hr2) 294 return vector 295 296 def eaccsd_diag(self): 297 if not getattr(self, 'imds', None): 298 self.imds = _IMDS(self) 299 if not self.imds.made_ea_imds: 300 self.imds.make_ea() 301 imds = self.imds 302 303 t1, t2 = self.t1, self.t2 304 nocc, nvir = t1.shape 305 306 Hr1 = np.diag(imds.Fvv) 307 Hr2 = np.zeros((nocc,nvir,nvir),dtype=t1.dtype) 308 for a in range(nvir): 309 _Wvvvva = np.array(imds.Wvvvv[a]) 310 for b in range(a): 311 for j in range(nocc): 312 Hr2[j,a,b] += imds.Fvv[a,a] 313 Hr2[j,a,b] += imds.Fvv[b,b] 314 Hr2[j,a,b] += -imds.Foo[j,j] 315 Hr2[j,a,b] += imds.Wovvo[j,b,b,j] 316 Hr2[j,a,b] += imds.Wovvo[j,a,a,j] 317 Hr2[j,a,b] += 0.5*(_Wvvvva[b,a,b]-_Wvvvva[b,b,a]) 318 Hr2[j,a,b] -= 0.5*(np.dot(imds.Woovv[:,j,a,b],t2[:,j,a,b]) - 319 np.dot(imds.Woovv[:,j,b,a],t2[:,j,a,b])) 320 321 vector = self.amplitudes_to_vector_ea(Hr1,Hr2) 322 return vector 323 324 def vector_to_amplitudes_ea(self,vector): 325 nocc = self.nocc 326 nvir = self.nmo - nocc 327 r1 = vector[:nvir].copy() 328 r2 = np.zeros((nocc,nvir,nvir), vector.dtype) 329 index = nvir 330 for i in range(nocc): 331 for a in range(nvir): 332 for b in range(a): 333 r2[i,a,b] = vector[index] 334 r2[i,b,a] = -vector[index] 335 index += 1 336 return [r1,r2] 337 338 def amplitudes_to_vector_ea(self,r1,r2): 339 nocc = self.nocc 340 nvir = self.nmo - nocc 341 size = nvir + nvir*(nvir-1)/2*nocc 342 vector = np.zeros(size, r1.dtype) 343 vector[:nvir] = r1.copy() 344 index = nvir 345 for i in range(nocc): 346 for a in range(nvir): 347 for b in range(a): 348 vector[index] = r2[i,a,b] 349 index += 1 350 return vector 351 352 def eeccsd_matvec(self,vector): 353 # Ref: Wang, Tu, and Wang, J. Chem. Theory Comput. 10, 5567 (2014) Eqs.(9)-(10) 354 # Note: Last line in Eq. (10) is superfluous. 355 # See, e.g. Gwaltney, Nooijen, and Barlett, Chem. Phys. Lett. 248, 189 (1996) 356 if not getattr(self, 'imds', None): 357 self.imds = _IMDS(self) 358 if not self.imds.made_ee_imds: 359 self.imds.make_ee() 360 imds = self.imds 361 362 r1,r2 = self.vector_to_amplitudes_ee(vector) 363 364 # Eq. (9) 365 Hr1 = einsum('ae,ie->ia',imds.Fvv,r1) 366 Hr1 += -einsum('mi,ma->ia',imds.Foo,r1) 367 Hr1 += einsum('me,imae->ia',imds.Fov,r2) 368 Hr1 += einsum('maei,me->ia',imds.Wovvo,r1) 369 Hr1 += -0.5*einsum('mnie,mnae->ia',imds.Wooov,r2) 370 Hr1 += 0.5*einsum('amef,imef->ia',imds.Wvovv,r2) 371 # Eq. (10) 372 tmpab = einsum('be,ijae->ijab',imds.Fvv,r2) 373 tmpab += -0.5*einsum('mnef,ijae,mnbf->ijab',imds.Woovv,self.t2,r2) 374 tmpab += -einsum('mbij,ma->ijab',imds.Wovoo,r1) 375 tmpab += -einsum('amef,ijfb,me->ijab',imds.Wvovv,self.t2,r1) 376 tmpij = -einsum('mj,imab->ijab',imds.Foo,r2) 377 tmpij += -0.5*einsum('mnef,imab,jnef->ijab',imds.Woovv,self.t2,r2) 378 tmpij += einsum('abej,ie->ijab',imds.Wvvvo,r1) 379 tmpij += einsum('mnie,njab,me->ijab',imds.Wooov,self.t2,r1) 380 381 tmpabij = einsum('mbej,imae->ijab',imds.Wovvo,r2) 382 383 Hr2 = (tmpab - tmpab.transpose(0,1,3,2) 384 + tmpij - tmpij.transpose(1,0,2,3) 385 + 0.5*einsum('mnij,mnab->ijab',imds.Woooo,r2) 386 + 0.5*einsum('abef,ijef->ijab',imds.Wvvvv,r2) 387 + tmpabij - tmpabij.transpose(0,1,3,2) 388 - tmpabij.transpose(1,0,2,3) + tmpabij.transpose(1,0,3,2) ) 389 390 vector = self.amplitudes_to_vector_ee(Hr1,Hr2) 391 return vector 392 393 def eeccsd_diag(self): 394 if not getattr(self, 'imds', None): 395 self.imds = _IMDS(self) 396 if not self.imds.made_ee_imds: 397 self.imds.make_ee() 398 imds = self.imds 399 400 t1, t2 = self.t1, self.t2 401 nocc, nvir = t1.shape 402 403 Hr1 = np.zeros((nocc,nvir), dtype=t1.dtype) 404 Hr2 = np.zeros((nocc,nocc,nvir,nvir), dtype=t1.dtype) 405 for i in range(nocc): 406 for a in range(nvir): 407 Hr1[i,a] = imds.Fvv[a,a] - imds.Foo[i,i] + imds.Wovvo[i,a,a,i] 408 for a in range(nvir): 409 tmp = 0.5*(np.einsum('ijeb,ijbe->ijb', imds.Woovv, t2) - 410 np.einsum('jieb,ijbe->ijb', imds.Woovv, t2)) 411 Hr2[:,:,:,a] += imds.Fvv[a,a] + tmp 412 Hr2[:,:,a,:] += imds.Fvv[a,a] + tmp 413 _Wvvvva = np.array(imds.Wvvvv[a]) 414 for b in range(a): 415 Hr2[:,:,a,b] += 0.5*(_Wvvvva[b,a,b]-_Wvvvva[b,b,a]) 416 for i in range(nocc): 417 tmp = imds.Wovvo[i,a,a,i] 418 Hr2[:,i,:,a] += tmp 419 Hr2[i,:,:,a] += tmp 420 Hr2[:,i,a,:] += tmp 421 Hr2[i,:,a,:] += tmp 422 for i in range(nocc): 423 tmp = 0.5*(np.einsum('kjab,jkab->jab', imds.Woovv, t2) - 424 np.einsum('kjba,jkab->jab', imds.Woovv, t2)) 425 Hr2[:,i,:,:] += -imds.Foo[i,i] + tmp 426 Hr2[i,:,:,:] += -imds.Foo[i,i] + tmp 427 for j in range(i): 428 Hr2[i,j,:,:] += 0.5*(imds.Woooo[i,j,i,j]-imds.Woooo[j,i,i,j]) 429 430 vector = self.amplitudes_to_vector_ee(Hr1,Hr2) 431 return vector 432 433 def vector_to_amplitudes_ee(self,vector): 434 nocc = self.nocc 435 nvir = self.nmo - nocc 436 r1 = vector[:nocc*nvir].copy().reshape((nocc,nvir)) 437 r2 = np.zeros((nocc,nocc,nvir,nvir), vector.dtype) 438 index = nocc*nvir 439 for i in range(nocc): 440 for j in range(i): 441 for a in range(nvir): 442 for b in range(a): 443 r2[i,j,a,b] = vector[index] 444 r2[j,i,a,b] = -vector[index] 445 r2[i,j,b,a] = -vector[index] 446 r2[j,i,b,a] = vector[index] 447 index += 1 448 return [r1,r2] 449 450 def amplitudes_to_vector_ee(self,r1,r2): 451 nocc = self.nocc 452 nvir = self.nmo - nocc 453 size = nocc*nvir + nocc*(nocc-1)/2*nvir*(nvir-1)/2 454 vector = np.zeros(size, r1.dtype) 455 vector[:nocc*nvir] = r1.copy().reshape(nocc*nvir) 456 index = nocc*nvir 457 for i in range(nocc): 458 for j in range(i): 459 for a in range(nvir): 460 for b in range(a): 461 vector[index] = r2[i,j,a,b] 462 index += 1 463 return vector 464 465 def amplitudes_from_rccsd(self, t1, t2): 466 '''Convert spatial orbital T1,T2 to spin-orbital T1,T2''' 467 nocc, nvir = t1.shape 468 nocc2 = nocc * 2 469 nvir2 = nvir * 2 470 t1s = np.zeros((nocc2,nvir2)) 471 t1s[:nocc,:nvir] = t1 472 t1s[nocc:,nvir:] = t1 473 474 t2s = np.zeros((nocc2,nocc2,nvir2,nvir2)) 475 t2s[:nocc,nocc:,:nvir,nvir:] = t2 476 t2s[nocc:,:nocc,nvir:,:nvir] = t2 477 t2s[:nocc,nocc:,nvir:,:nvir] =-t2.transpose(0,1,3,2) 478 t2s[nocc:,:nocc,:nvir,nvir:] =-t2.transpose(0,1,3,2) 479 t2s[:nocc,:nocc,:nvir,:nvir] = t2 - t2.transpose(0,1,3,2) 480 t2s[nocc:,nocc:,nvir:,nvir:] = t2 - t2.transpose(0,1,3,2) 481 return t1s, t2s 482 483 484class _PhysicistsERIs: 485 def __init__(self, cc, mo_coeff=None, method='incore', 486 ao2mofn=ao2mo.outcore.general_iofree): 487 cput0 = (logger.process_clock(), logger.perf_counter()) 488 moidx = cc.get_frozen_mask() 489 if mo_coeff is None: 490 self.mo_coeff = mo_coeff = [cc.mo_coeff[0][:,moidx[0]], 491 cc.mo_coeff[1][:,moidx[1]]] 492 else: 493 self.mo_coeff = mo_coeff = [mo_coeff[0][:,moidx[0]], 494 mo_coeff[1][:,moidx[1]]] 495 496 nocc = cc.nocc 497 nmo = cc.nmo 498 nvir = nmo - nocc 499 mem_incore = nmo**4*2 * 8/1e6 500 mem_now = lib.current_memory()[0] 501 502 self.fock, so_coeff, spin = uspatial2spin(cc, moidx, mo_coeff) 503 504 log = logger.Logger(cc.stdout, cc.verbose) 505 if (method == 'incore' and (mem_incore+mem_now < cc.max_memory) 506 or cc.mol.incore_anyway): 507 eri = ao2mofn(cc._scf.mol, (so_coeff,so_coeff,so_coeff,so_coeff), compact=0) 508 if mo_coeff[0].dtype == np.float: eri = eri.real 509 eri = eri.reshape((nmo,)*4) 510 for i in range(nmo): 511 for j in range(i): 512 if spin[i] != spin[j]: 513 eri[i,j,:,:] = eri[j,i,:,:] = 0. 514 eri[:,:,i,j] = eri[:,:,j,i] = 0. 515 eri1 = eri - eri.transpose(0,3,2,1) 516 eri1 = eri1.transpose(0,2,1,3) 517 518 self.dtype = eri1.dtype 519 self.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy() 520 self.ooov = eri1[:nocc,:nocc,:nocc,nocc:].copy() 521 self.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy() 522 self.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy() 523 self.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy() 524 self.ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy() 525 self.vvvv = eri1[nocc:,nocc:,nocc:,nocc:].copy() 526 else: 527 self.feri1 = lib.H5TmpFile() 528 orbo = so_coeff[:,:nocc] 529 orbv = so_coeff[:,nocc:] 530 if mo_coeff[0].dtype == np.complex128: ds_type = 'c16' 531 else: ds_type = 'f8' 532 self.oooo = self.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), ds_type) 533 self.ooov = self.feri1.create_dataset('ooov', (nocc,nocc,nocc,nvir), ds_type) 534 self.ovoo = self.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), ds_type) 535 self.oovv = self.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), ds_type) 536 self.ovov = self.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), ds_type) 537 self.ovvv = self.feri1.create_dataset('ovvv', (nocc,nvir,nvir,nvir), ds_type) 538 self.vvvv = self.feri1.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type) 539 540 cput1 = logger.process_clock(), logger.perf_counter() 541 # <ij||pq> = <ij|pq> - <ij|qp> = (ip|jq) - (iq|jp) 542 buf = ao2mofn(cc._scf.mol, (orbo,so_coeff,orbo,so_coeff), compact=0) 543 if mo_coeff[0].dtype == np.float: buf = buf.real 544 buf = buf.reshape((nocc,nmo,nocc,nmo)) 545 for i in range(nocc): 546 for p in range(nmo): 547 if spin[i] != spin[p]: 548 buf[i,p,:,:] = 0. 549 buf[:,:,i,p] = 0. 550 buf1 = buf - buf.transpose(0,3,2,1) 551 buf1 = buf1.transpose(0,2,1,3) 552 cput1 = log.timer_debug1('transforming oopq', *cput1) 553 self.dtype = buf1.dtype 554 self.oooo[:,:,:,:] = buf1[:,:,:nocc,:nocc] 555 self.ooov[:,:,:,:] = buf1[:,:,:nocc,nocc:] 556 self.oovv[:,:,:,:] = buf1[:,:,nocc:,nocc:] 557 558 cput1 = logger.process_clock(), logger.perf_counter() 559 # <ia||pq> = <ia|pq> - <ia|qp> = (ip|aq) - (iq|ap) 560 buf = ao2mofn(cc._scf.mol, (orbo,so_coeff,orbv,so_coeff), compact=0) 561 if mo_coeff[0].dtype == np.float: buf = buf.real 562 buf = buf.reshape((nocc,nmo,nvir,nmo)) 563 for p in range(nmo): 564 for i in range(nocc): 565 if spin[i] != spin[p]: 566 buf[i,p,:,:] = 0. 567 for a in range(nvir): 568 if spin[nocc+a] != spin[p]: 569 buf[:,:,a,p] = 0. 570 buf1 = buf - buf.transpose(0,3,2,1) 571 buf1 = buf1.transpose(0,2,1,3) 572 cput1 = log.timer_debug1('transforming ovpq', *cput1) 573 self.ovoo[:,:,:,:] = buf1[:,:,:nocc,:nocc] 574 self.ovov[:,:,:,:] = buf1[:,:,:nocc,nocc:] 575 self.ovvv[:,:,:,:] = buf1[:,:,nocc:,nocc:] 576 577 for a in range(nvir): 578 orbva = orbv[:,a].reshape(-1,1) 579 buf = ao2mofn(cc._scf.mol, (orbva,orbv,orbv,orbv), compact=0) 580 if mo_coeff[0].dtype == np.float: buf = buf.real 581 buf = buf.reshape((1,nvir,nvir,nvir)) 582 for b in range(nvir): 583 if spin[nocc+a] != spin[nocc+b]: 584 buf[0,b,:,:] = 0. 585 for c in range(nvir): 586 if spin[nocc+b] != spin[nocc+c]: 587 buf[:,:,b,c] = buf[:,:,c,b] = 0. 588 buf1 = buf - buf.transpose(0,3,2,1) 589 buf1 = buf1.transpose(0,2,1,3) 590 self.vvvv[a] = buf1[:] 591 592 cput1 = log.timer_debug1('transforming vvvv', *cput1) 593 594 log.timer('CCSD integral transformation', *cput0) 595 596 597def uspatial2spin(cc, moidx, mo_coeff): 598 '''Convert the results of an unrestricted mean-field calculation to spin-orbital form. 599 600 Spin-orbital ordering is determined by orbital energy without regard for spin. 601 602 Returns: 603 fock : (nso,nso) ndarray 604 The Fock matrix in the basis of spin-orbitals 605 so_coeff : (nao, nso) ndarray 606 The matrix of spin-orbital coefficients in the AO basis 607 spin : (nso,) ndarary 608 The spin (0 or 1) of each spin-orbital 609 ''' 610 611 dm = cc._scf.make_rdm1(cc.mo_coeff, cc.mo_occ) 612 fockao = cc._scf.get_hcore() + cc._scf.get_veff(cc.mol, dm) 613 fockab = list() 614 for a in range(2): 615 fockab.append( reduce(numpy.dot, (mo_coeff[a].T, fockao[a], mo_coeff[a])) ) 616 617 nocc = cc.nocc 618 nao = cc.mo_coeff[0].shape[0] 619 nmo = cc.nmo 620 so_coeff = np.zeros((nao, nmo), dtype=mo_coeff[0].dtype) 621 nocc_a = int(sum(cc.mo_occ[0]*moidx[0])) 622 nocc_b = int(sum(cc.mo_occ[1]*moidx[1])) 623 nmo_a = fockab[0].shape[0] 624 nmo_b = fockab[1].shape[0] 625 nvir_a = nmo_a - nocc_a 626 #nvir_b = nmo_b - nocc_b 627 oa = range(0,nocc_a) 628 ob = range(nocc_a,nocc) 629 va = range(nocc,nocc+nvir_a) 630 vb = range(nocc+nvir_a,nmo) 631 spin = np.zeros(nmo, dtype=int) 632 spin[oa] = 0 633 spin[ob] = 1 634 spin[va] = 0 635 spin[vb] = 1 636 so_coeff[:,oa] = mo_coeff[0][:,:nocc_a] 637 so_coeff[:,ob] = mo_coeff[1][:,:nocc_b] 638 so_coeff[:,va] = mo_coeff[0][:,nocc_a:nmo_a] 639 so_coeff[:,vb] = mo_coeff[1][:,nocc_b:nmo_b] 640 641 fock = np.zeros((nmo, nmo), dtype=fockab[0].dtype) 642 fock[np.ix_(oa,oa)] = fockab[0][:nocc_a,:nocc_a] 643 fock[np.ix_(oa,va)] = fockab[0][:nocc_a,nocc_a:] 644 fock[np.ix_(va,oa)] = fockab[0][nocc_a:,:nocc_a] 645 fock[np.ix_(va,va)] = fockab[0][nocc_a:,nocc_a:] 646 fock[np.ix_(ob,ob)] = fockab[1][:nocc_b,:nocc_b] 647 fock[np.ix_(ob,vb)] = fockab[1][:nocc_b,nocc_b:] 648 fock[np.ix_(vb,ob)] = fockab[1][nocc_b:,:nocc_b] 649 fock[np.ix_(vb,vb)] = fockab[1][nocc_b:,nocc_b:] 650 651# Do not sort because it's different to the orbital ordering generated by 652# get_frozen_mask function in AO-direct vvvv contraction 653# idxo = np.diagonal(fock[:nocc,:nocc]).argsort() 654# idxv = nocc + np.diagonal(fock[nocc:,nocc:]).argsort() 655# idx = np.concatenate((idxo,idxv)) 656# spin = spin[idx] 657# so_coeff = so_coeff[:,idx] 658# fock = fock[:, idx][idx] 659 660 return fock, so_coeff, spin 661 662 663class _IMDS: 664 # Exactly the same as RCCSD IMDS except 665 # -- rintermediates --> uintermediates 666 # -- Loo, Lvv, cc_Fov --> Foo, Fvv, Fov 667 # -- One less 2-virtual intermediate 668 def __init__(self, cc): 669 self.cc = cc 670 self._made_shared = False 671 self.made_ip_imds = False 672 self.made_ea_imds = False 673 self.made_ee_imds = False 674 675 def _make_shared(self): 676 cput0 = (logger.process_clock(), logger.perf_counter()) 677 log = logger.Logger(self.cc.stdout, self.cc.verbose) 678 679 t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris 680 self.Foo = imd.Foo(t1,t2,eris) 681 self.Fvv = imd.Fvv(t1,t2,eris) 682 self.Fov = imd.Fov(t1,t2,eris) 683 684 # 2 virtuals 685 self.Wovvo = imd.Wovvo(t1,t2,eris) 686 self.Woovv = eris.oovv 687 688 log.timer('EOM-CCSD shared intermediates', *cput0) 689 690 def make_ip(self): 691 if self._made_shared is False: 692 self._make_shared() 693 self._made_shared = True 694 695 cput0 = (logger.process_clock(), logger.perf_counter()) 696 log = logger.Logger(self.cc.stdout, self.cc.verbose) 697 698 t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris 699 700 # 0 or 1 virtuals 701 self.Woooo = imd.Woooo(t1,t2,eris) 702 self.Wooov = imd.Wooov(t1,t2,eris) 703 self.Wovoo = imd.Wovoo(t1,t2,eris) 704 705 self.made_ip_imds = True 706 log.timer('EOM-CCSD IP intermediates', *cput0) 707 708 def make_ea(self): 709 if self._made_shared is False: 710 self._make_shared() 711 self._made_shared = True 712 713 cput0 = (logger.process_clock(), logger.perf_counter()) 714 log = logger.Logger(self.cc.stdout, self.cc.verbose) 715 716 t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris 717 718 # 3 or 4 virtuals 719 self.Wvovv = imd.Wvovv(t1,t2,eris) 720 self.Wvvvv = imd.Wvvvv(t1,t2,eris) 721 self.Wvvvo = imd.Wvvvo(t1,t2,eris,self.Wvvvv) 722 723 self.made_ea_imds = True 724 log.timer('EOM-CCSD EA intermediates', *cput0) 725 726 def make_ee(self): 727 if self._made_shared is False: 728 self._make_shared() 729 self._made_shared = True 730 731 cput0 = (logger.process_clock(), logger.perf_counter()) 732 log = logger.Logger(self.cc.stdout, self.cc.verbose) 733 734 t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris 735 736 if self.made_ip_imds is False: 737 # 0 or 1 virtuals 738 self.Woooo = imd.Woooo(t1,t2,eris) 739 self.Wooov = imd.Wooov(t1,t2,eris) 740 self.Wovoo = imd.Wovoo(t1,t2,eris) 741 if self.made_ea_imds is False: 742 # 3 or 4 virtuals 743 self.Wvovv = imd.Wvovv(t1,t2,eris) 744 self.Wvvvv = imd.Wvvvv(t1,t2,eris) 745 self.Wvvvo = imd.Wvvvo(t1,t2,eris,self.Wvvvv) 746 747 self.made_ee_imds = True 748 log.timer('EOM-CCSD EE intermediates', *cput0) 749 750if __name__ == '__main__': 751 from pyscf import scf 752 from pyscf import gto 753 mol = gto.Mole() 754 mol.atom = [['O', (0., 0., 0.)], 755 ['O', (1.21, 0., 0.)]] 756 mol.basis = 'cc-pvdz' 757 mol.spin = 2 758 mol.build() 759 mf = scf.UHF(mol) 760 print(mf.scf()) 761 762 # Freeze 1s electrons 763 frozen = [[0,1], [0,1]] 764 # also acceptable 765 #frozen = 4 766 ucc = UCCSD(mf, frozen=frozen) 767 ecc, t1, t2 = ucc.kernel() 768 print(ecc - -0.3486987472235819) 769 exit() 770 771 mol = gto.Mole() 772 mol.atom = [ 773 [8 , (0. , 0. , 0.)], 774 [1 , (0. , -0.757 , 0.587)], 775 [1 , (0. , 0.757 , 0.587)]] 776 mol.basis = 'cc-pvdz' 777 mol.spin = 0 778 mol.build() 779 mf = scf.UHF(mol) 780 print(mf.scf()) 781 782 mycc = UCCSD(mf) 783 ecc, t1, t2 = mycc.kernel() 784 print(ecc - -0.2133432712431435) 785 e,v = mycc.ipccsd(nroots=8) 786 print(e[0] - 0.4335604332073799) 787 print(e[2] - 0.5187659896045407) 788 print(e[4] - 0.6782876002229172) 789 790 mycc.verbose = 5 791 e,v = mycc.eaccsd(nroots=8) 792 print(e[0] - 0.16737886338859731) 793 print(e[2] - 0.24027613852009164) 794 print(e[4] - 0.51006797826488071) 795 print("e=", e) 796 797 e,v = mycc.eeccsd(nroots=4) 798 print(e[0] - 0.2757159395886167) 799 print(e[1] - 0.2757159395886167) 800 print(e[2] - 0.2757159395886167) 801 print(e[3] - 0.3005716731825082) 802