#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from functools import reduce
import numpy
import numpy as np
from pyscf import lib
from pyscf import ao2mo
from pyscf.lib import logger
from pyscf.cc import ccsd
from pyscf.cc import uccsd
from pyscf.mp import ump2
from pyscf.cc import gintermediates as imd
#einsum = np.einsum
einsum = lib.einsum
# This is unrestricted (U)CCSD, i.e. spin-orbital form.
[docs]
def update_amps(cc, t1, t2, eris):
nocc, nvir = t1.shape
fock = eris.fock
fov = fock[:nocc,nocc:]
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
tau = imd.make_tau(t2, t1, t1)
Fvv = imd.cc_Fvv(t1, t2, eris)
Foo = imd.cc_Foo(t1, t2, eris)
Fov = imd.cc_Fov(t1, t2, eris)
Woooo = imd.cc_Woooo(t1, t2, eris)
Wvvvv = imd.cc_Wvvvv(t1, t2, eris)
Wovvo = imd.cc_Wovvo(t1, t2, eris)
# Move energy terms to the other side
Fvv -= np.diag(np.diag(fvv))
Foo -= np.diag(np.diag(foo))
# T1 equation
t1new = np.array(fov).conj()
t1new += einsum('ie,ae->ia', t1, Fvv)
t1new += -einsum('ma,mi->ia', t1, Foo)
t1new += einsum('imae,me->ia', t2, Fov)
t1new += -einsum('nf,naif->ia', t1, eris.ovov)
t1new += -0.5*einsum('imef,maef->ia', t2, eris.ovvv)
t1new += -0.5*einsum('mnae,mnie->ia', t2, eris.ooov)
# T2 equation
t2new = np.array(eris.oovv).conj()
Ftmp = Fvv - 0.5*einsum('mb,me->be', t1, Fov)
tmp = einsum('ijae,be->ijab', t2, Ftmp)
t2new += (tmp - tmp.transpose(0,1,3,2))
Ftmp = Foo + 0.5*einsum('je,me->mj', t1, Fov)
tmp = einsum('imab,mj->ijab', t2, Ftmp)
t2new -= (tmp - tmp.transpose(1,0,2,3))
t2new += 0.5*einsum('mnab,mnij->ijab', tau, Woooo)
t2new += 0.5*einsum('ijef,abef->ijab', tau, Wvvvv)
tmp = einsum('imae,mbej->ijab', t2, Wovvo)
tmp -= -einsum('ie,ma,mbje->ijab', t1, t1, eris.ovov)
t2new += (tmp - tmp.transpose(0,1,3,2)
- tmp.transpose(1,0,2,3) + tmp.transpose(1,0,3,2) )
tmp = einsum('ie,jeba->ijab', t1, np.array(eris.ovvv).conj())
t2new += (tmp - tmp.transpose(1,0,2,3))
tmp = einsum('ma,mbij->ijab', t1, eris.ovoo)
t2new -= (tmp - tmp.transpose(0,1,3,2))
mo_e = eris.fock.diagonal()
eia = mo_e[:nocc,None] - mo_e[None,nocc:] - cc.level_shift
eijab = lib.direct_sum('ia,jb->ijab', eia, eia)
t1new /= eia
t2new /= eijab
return t1new, t2new
[docs]
def energy(cc, t1, t2, eris):
nocc, nvir = t1.shape
fock = eris.fock
e = einsum('ia,ia', fock[:nocc,nocc:], t1)
e += 0.25*np.einsum('ijab,ijab', t2, eris.oovv)
e += 0.5 *np.einsum('ia,jb,ijab', t1, t1, eris.oovv)
return e.real
get_frozen_mask = ump2.get_frozen_mask
[docs]
class UCCSD(ccsd.CCSD):
def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None):
ccsd.CCSD.__init__(self, mf, frozen, mo_coeff, mo_occ)
# Spin-orbital CCSD needs a stricter tolerance than spatial-orbital
self.conv_tol_normt = 1e-6
@property
def nocc(self):
nocca, noccb = self.get_nocc()
return nocca + noccb
@property
def nmo(self):
nmoa, nmob = self.get_nmo()
return nmoa + nmob
get_nocc = uccsd.get_nocc
get_nmo = uccsd.get_nmo
get_frozen_mask = get_frozen_mask
[docs]
def init_amps(self, eris):
time0 = logger.process_clock(), logger.perf_counter()
mo_e = eris.fock.diagonal()
nocc = self.nocc
eia = mo_e[:nocc,None] - mo_e[None,nocc:]
eijab = lib.direct_sum('ia,jb->ijab',eia,eia)
t1 = eris.fock[:nocc,nocc:] / eia
eris_oovv = np.array(eris.oovv)
t2 = eris_oovv/eijab
self.emp2 = 0.25*einsum('ijab,ijab',t2,eris_oovv.conj()).real
logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2)
logger.timer(self, 'init mp2', *time0)
return self.emp2, t1, t2
energy = energy
[docs]
def kernel(self, t1=None, t2=None, eris=None, mbpt2=False):
return self.ccsd(t1, t2, eris, mbpt2)
[docs]
def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False):
'''Ground-state unrestricted (U)CCSD.
Kwargs:
mbpt2 : bool
Use one-shot MBPT2 approximation to CCSD.
'''
if mbpt2:
#cctyp = 'MBPT2'
#self.e_corr, self.t1, self.t2 = self.init_amps(eris)
raise NotImplementedError
if eris is None:
eris = self.ao2mo(self.mo_coeff)
self.eris = eris
return ccsd.CCSD.ccsd(self, t1, t2, eris)
[docs]
def ao2mo(self, mo_coeff=None):
return _PhysicistsERIs(self, mo_coeff)
[docs]
def update_amps(self, t1, t2, eris):
return update_amps(self, t1, t2, eris)
[docs]
def nip(self):
nocc = self.nocc
nvir = self.nmo - nocc
self._nip = nocc + nocc*(nocc-1)//2*nvir
return self._nip
[docs]
def nea(self):
nocc = self.nocc
nvir = self.nmo - nocc
self._nea = nvir + nocc*nvir*(nvir-1)//2
return self._nea
[docs]
def nee(self):
nocc = self.nocc
nvir = self.nmo - nocc
self._nee = nocc*nvir + nocc*(nocc-1)//2*nvir*(nvir-1)//2
return self._nee
[docs]
def ipccsd_matvec(self, vector):
# Ref: Tu, Wang, and Li, J. Chem. Phys. 136, 174102 (2012) Eqs.(8)-(9)
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ip_imds:
self.imds.make_ip()
imds = self.imds
r1,r2 = self.vector_to_amplitudes_ip(vector)
# Eq. (8)
Hr1 = -einsum('mi,m->i',imds.Foo,r1)
Hr1 += einsum('me,mie->i',imds.Fov,r2)
Hr1 += -0.5*einsum('nmie,mne->i',imds.Wooov,r2)
# Eq. (9)
Hr2 = einsum('ae,ije->ija',imds.Fvv,r2)
tmp1 = einsum('mi,mja->ija',imds.Foo,r2)
Hr2 += (-tmp1 + tmp1.transpose(1,0,2))
Hr2 += -einsum('maji,m->ija',imds.Wovoo,r1)
Hr2 += 0.5*einsum('mnij,mna->ija',imds.Woooo,r2)
tmp2 = einsum('maei,mje->ija',imds.Wovvo,r2)
Hr2 += (tmp2 - tmp2.transpose(1,0,2))
Hr2 += 0.5*einsum('mnef,ijae,mnf->ija',imds.Woovv,self.t2,r2)
vector = self.amplitudes_to_vector_ip(Hr1,Hr2)
return vector
[docs]
def ipccsd_diag(self):
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ip_imds:
self.imds.make_ip()
imds = self.imds
t1, t2 = self.t1, self.t2
nocc, nvir = t1.shape
Hr1 = -np.diag(imds.Foo)
Hr2 = np.zeros((nocc,nocc,nvir),dtype=t1.dtype)
for i in range(nocc):
for j in range(nocc):
for a in range(nvir):
Hr2[i,j,a] += imds.Fvv[a,a]
Hr2[i,j,a] += -imds.Foo[i,i]
Hr2[i,j,a] += -imds.Foo[j,j]
Hr2[i,j,a] += 0.5*(imds.Woooo[i,j,i,j]-imds.Woooo[j,i,i,j])
Hr2[i,j,a] += imds.Wovvo[i,a,a,i]
Hr2[i,j,a] += imds.Wovvo[j,a,a,j]
Hr2[i,j,a] += 0.5*(np.dot(imds.Woovv[i,j,:,a],t2[i,j,a,:])
-np.dot(imds.Woovv[j,i,:,a],t2[i,j,a,:]))
vector = self.amplitudes_to_vector_ip(Hr1,Hr2)
return vector
[docs]
def vector_to_amplitudes_ip(self,vector):
nocc = self.nocc
nvir = self.nmo - nocc
r1 = vector[:nocc].copy()
r2 = np.zeros((nocc,nocc,nvir), vector.dtype)
index = nocc
for i in range(nocc):
for j in range(i):
for a in range(nvir):
r2[i,j,a] = vector[index]
r2[j,i,a] = -vector[index]
index += 1
return [r1,r2]
[docs]
def amplitudes_to_vector_ip(self,r1,r2):
nocc = self.nocc
nvir = self.nmo - nocc
size = nocc + nocc*(nocc-1)//2*nvir
vector = np.zeros(size, r1.dtype)
vector[:nocc] = r1.copy()
index = nocc
for i in range(nocc):
for j in range(i):
for a in range(nvir):
vector[index] = r2[i,j,a]
index += 1
return vector
[docs]
def eaccsd_matvec(self,vector):
# Ref: Nooijen and Bartlett, J. Chem. Phys. 102, 3629 (1994) Eqs.(30)-(31)
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ea_imds:
self.imds.make_ea()
imds = self.imds
r1,r2 = self.vector_to_amplitudes_ea(vector)
# Eq. (30)
Hr1 = einsum('ac,c->a',imds.Fvv,r1)
Hr1 += einsum('ld,lad->a',imds.Fov,r2)
Hr1 += 0.5*einsum('alcd,lcd->a',imds.Wvovv,r2)
# Eq. (31)
Hr2 = einsum('abcj,c->jab',imds.Wvvvo,r1)
tmp1 = einsum('ac,jcb->jab',imds.Fvv,r2)
Hr2 += (tmp1 - tmp1.transpose(0,2,1))
Hr2 += -einsum('lj,lab->jab',imds.Foo,r2)
tmp2 = einsum('lbdj,lad->jab',imds.Wovvo,r2)
Hr2 += (tmp2 - tmp2.transpose(0,2,1))
nvir = self.nmo-self.nocc
for a in range(nvir):
Hr2[:,a,:] += 0.5*einsum('bcd,jcd->jb',imds.Wvvvv[a],r2)
Hr2 += -0.5*einsum('klcd,lcd,kjab->jab',imds.Woovv,r2,self.t2)
vector = self.amplitudes_to_vector_ea(Hr1,Hr2)
return vector
[docs]
def eaccsd_diag(self):
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ea_imds:
self.imds.make_ea()
imds = self.imds
t1, t2 = self.t1, self.t2
nocc, nvir = t1.shape
Hr1 = np.diag(imds.Fvv)
Hr2 = np.zeros((nocc,nvir,nvir),dtype=t1.dtype)
for a in range(nvir):
_Wvvvva = np.array(imds.Wvvvv[a])
for b in range(a):
for j in range(nocc):
Hr2[j,a,b] += imds.Fvv[a,a]
Hr2[j,a,b] += imds.Fvv[b,b]
Hr2[j,a,b] += -imds.Foo[j,j]
Hr2[j,a,b] += imds.Wovvo[j,b,b,j]
Hr2[j,a,b] += imds.Wovvo[j,a,a,j]
Hr2[j,a,b] += 0.5*(_Wvvvva[b,a,b]-_Wvvvva[b,b,a])
Hr2[j,a,b] -= 0.5*(np.dot(imds.Woovv[:,j,a,b],t2[:,j,a,b]) -
np.dot(imds.Woovv[:,j,b,a],t2[:,j,a,b]))
vector = self.amplitudes_to_vector_ea(Hr1,Hr2)
return vector
[docs]
def vector_to_amplitudes_ea(self,vector):
nocc = self.nocc
nvir = self.nmo - nocc
r1 = vector[:nvir].copy()
r2 = np.zeros((nocc,nvir,nvir), vector.dtype)
index = nvir
for i in range(nocc):
for a in range(nvir):
for b in range(a):
r2[i,a,b] = vector[index]
r2[i,b,a] = -vector[index]
index += 1
return [r1,r2]
[docs]
def amplitudes_to_vector_ea(self,r1,r2):
nocc = self.nocc
nvir = self.nmo - nocc
size = nvir + nvir*(nvir-1)//2*nocc
vector = np.zeros(size, r1.dtype)
vector[:nvir] = r1.copy()
index = nvir
for i in range(nocc):
for a in range(nvir):
for b in range(a):
vector[index] = r2[i,a,b]
index += 1
return vector
[docs]
def eeccsd_matvec(self,vector):
# Ref: Wang, Tu, and Wang, J. Chem. Theory Comput. 10, 5567 (2014) Eqs.(9)-(10)
# Note: Last line in Eq. (10) is superfluous.
# See, e.g. Gwaltney, Nooijen, and Barlett, Chem. Phys. Lett. 248, 189 (1996)
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ee_imds:
self.imds.make_ee()
imds = self.imds
r1,r2 = self.vector_to_amplitudes_ee(vector)
# Eq. (9)
Hr1 = einsum('ae,ie->ia',imds.Fvv,r1)
Hr1 += -einsum('mi,ma->ia',imds.Foo,r1)
Hr1 += einsum('me,imae->ia',imds.Fov,r2)
Hr1 += einsum('maei,me->ia',imds.Wovvo,r1)
Hr1 += -0.5*einsum('mnie,mnae->ia',imds.Wooov,r2)
Hr1 += 0.5*einsum('amef,imef->ia',imds.Wvovv,r2)
# Eq. (10)
tmpab = einsum('be,ijae->ijab',imds.Fvv,r2)
tmpab += -0.5*einsum('mnef,ijae,mnbf->ijab',imds.Woovv,self.t2,r2)
tmpab += -einsum('mbij,ma->ijab',imds.Wovoo,r1)
tmpab += -einsum('amef,ijfb,me->ijab',imds.Wvovv,self.t2,r1)
tmpij = -einsum('mj,imab->ijab',imds.Foo,r2)
tmpij += -0.5*einsum('mnef,imab,jnef->ijab',imds.Woovv,self.t2,r2)
tmpij += einsum('abej,ie->ijab',imds.Wvvvo,r1)
tmpij += einsum('mnie,njab,me->ijab',imds.Wooov,self.t2,r1)
tmpabij = einsum('mbej,imae->ijab',imds.Wovvo,r2)
Hr2 = (tmpab - tmpab.transpose(0,1,3,2)
+ tmpij - tmpij.transpose(1,0,2,3)
+ 0.5*einsum('mnij,mnab->ijab',imds.Woooo,r2)
+ 0.5*einsum('abef,ijef->ijab',imds.Wvvvv,r2)
+ tmpabij - tmpabij.transpose(0,1,3,2)
- tmpabij.transpose(1,0,2,3) + tmpabij.transpose(1,0,3,2) )
vector = self.amplitudes_to_vector_ee(Hr1,Hr2)
return vector
[docs]
def eeccsd_diag(self):
if not getattr(self, 'imds', None):
self.imds = _IMDS(self)
if not self.imds.made_ee_imds:
self.imds.make_ee()
imds = self.imds
t1, t2 = self.t1, self.t2
nocc, nvir = t1.shape
Hr1 = np.zeros((nocc,nvir), dtype=t1.dtype)
Hr2 = np.zeros((nocc,nocc,nvir,nvir), dtype=t1.dtype)
for i in range(nocc):
for a in range(nvir):
Hr1[i,a] = imds.Fvv[a,a] - imds.Foo[i,i] + imds.Wovvo[i,a,a,i]
for a in range(nvir):
tmp = 0.5*(np.einsum('ijeb,ijbe->ijb', imds.Woovv, t2) -
np.einsum('jieb,ijbe->ijb', imds.Woovv, t2))
Hr2[:,:,:,a] += imds.Fvv[a,a] + tmp
Hr2[:,:,a,:] += imds.Fvv[a,a] + tmp
_Wvvvva = np.array(imds.Wvvvv[a])
for b in range(a):
Hr2[:,:,a,b] += 0.5*(_Wvvvva[b,a,b]-_Wvvvva[b,b,a])
for i in range(nocc):
tmp = imds.Wovvo[i,a,a,i]
Hr2[:,i,:,a] += tmp
Hr2[i,:,:,a] += tmp
Hr2[:,i,a,:] += tmp
Hr2[i,:,a,:] += tmp
for i in range(nocc):
tmp = 0.5*(np.einsum('kjab,jkab->jab', imds.Woovv, t2) -
np.einsum('kjba,jkab->jab', imds.Woovv, t2))
Hr2[:,i,:,:] += -imds.Foo[i,i] + tmp
Hr2[i,:,:,:] += -imds.Foo[i,i] + tmp
for j in range(i):
Hr2[i,j,:,:] += 0.5*(imds.Woooo[i,j,i,j]-imds.Woooo[j,i,i,j])
vector = self.amplitudes_to_vector_ee(Hr1,Hr2)
return vector
[docs]
def vector_to_amplitudes_ee(self,vector):
nocc = self.nocc
nvir = self.nmo - nocc
r1 = vector[:nocc*nvir].copy().reshape((nocc,nvir))
r2 = np.zeros((nocc,nocc,nvir,nvir), vector.dtype)
index = nocc*nvir
for i in range(nocc):
for j in range(i):
for a in range(nvir):
for b in range(a):
r2[i,j,a,b] = vector[index]
r2[j,i,a,b] = -vector[index]
r2[i,j,b,a] = -vector[index]
r2[j,i,b,a] = vector[index]
index += 1
return [r1,r2]
[docs]
def amplitudes_to_vector_ee(self,r1,r2):
nocc = self.nocc
nvir = self.nmo - nocc
size = nocc*nvir + nocc*(nocc-1)//2*nvir*(nvir-1)//2
vector = np.zeros(size, r1.dtype)
vector[:nocc*nvir] = r1.copy().reshape(nocc*nvir)
index = nocc*nvir
for i in range(nocc):
for j in range(i):
for a in range(nvir):
for b in range(a):
vector[index] = r2[i,j,a,b]
index += 1
return vector
[docs]
def amplitudes_from_rccsd(self, t1, t2):
'''Convert spatial orbital T1,T2 to spin-orbital T1,T2'''
nocc, nvir = t1.shape
nocc2 = nocc * 2
nvir2 = nvir * 2
t1s = np.zeros((nocc2,nvir2))
t1s[:nocc,:nvir] = t1
t1s[nocc:,nvir:] = t1
t2s = np.zeros((nocc2,nocc2,nvir2,nvir2))
t2s[:nocc,nocc:,:nvir,nvir:] = t2
t2s[nocc:,:nocc,nvir:,:nvir] = t2
t2s[:nocc,nocc:,nvir:,:nvir] =-t2.transpose(0,1,3,2)
t2s[nocc:,:nocc,:nvir,nvir:] =-t2.transpose(0,1,3,2)
t2s[:nocc,:nocc,:nvir,:nvir] = t2 - t2.transpose(0,1,3,2)
t2s[nocc:,nocc:,nvir:,nvir:] = t2 - t2.transpose(0,1,3,2)
return t1s, t2s
class _PhysicistsERIs:
def __init__(self, cc, mo_coeff=None, method='incore',
ao2mofn=ao2mo.outcore.general_iofree):
cput0 = (logger.process_clock(), logger.perf_counter())
moidx = cc.get_frozen_mask()
if mo_coeff is None:
self.mo_coeff = mo_coeff = [cc.mo_coeff[0][:,moidx[0]],
cc.mo_coeff[1][:,moidx[1]]]
else:
self.mo_coeff = mo_coeff = [mo_coeff[0][:,moidx[0]],
mo_coeff[1][:,moidx[1]]]
nocc = cc.nocc
nmo = cc.nmo
nvir = nmo - nocc
mem_incore = nmo**4*2 * 8/1e6
mem_now = lib.current_memory()[0]
self.fock, so_coeff, spin = uspatial2spin(cc, moidx, mo_coeff)
log = logger.Logger(cc.stdout, cc.verbose)
if (method == 'incore' and (mem_incore+mem_now < cc.max_memory)
or cc.mol.incore_anyway):
eri = ao2mofn(cc._scf.mol, (so_coeff,so_coeff,so_coeff,so_coeff), compact=0)
if mo_coeff[0].dtype == np.double: eri = eri.real
eri = eri.reshape((nmo,)*4)
for i in range(nmo):
for j in range(i):
if spin[i] != spin[j]:
eri[i,j,:,:] = eri[j,i,:,:] = 0.
eri[:,:,i,j] = eri[:,:,j,i] = 0.
eri1 = eri - eri.transpose(0,3,2,1)
eri1 = eri1.transpose(0,2,1,3)
self.dtype = eri1.dtype
self.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy()
self.ooov = eri1[:nocc,:nocc,:nocc,nocc:].copy()
self.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy()
self.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy()
self.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy()
self.ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy()
self.vvvv = eri1[nocc:,nocc:,nocc:,nocc:].copy()
else:
self.feri1 = lib.H5TmpFile()
orbo = so_coeff[:,:nocc]
orbv = so_coeff[:,nocc:]
if mo_coeff[0].dtype == np.complex128: ds_type = 'c16'
else: ds_type = 'f8'
self.oooo = self.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), ds_type)
self.ooov = self.feri1.create_dataset('ooov', (nocc,nocc,nocc,nvir), ds_type)
self.ovoo = self.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), ds_type)
self.oovv = self.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), ds_type)
self.ovov = self.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), ds_type)
self.ovvv = self.feri1.create_dataset('ovvv', (nocc,nvir,nvir,nvir), ds_type)
self.vvvv = self.feri1.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type)
cput1 = logger.process_clock(), logger.perf_counter()
# <ij||pq> = <ij|pq> - <ij|qp> = (ip|jq) - (iq|jp)
buf = ao2mofn(cc._scf.mol, (orbo,so_coeff,orbo,so_coeff), compact=0)
if mo_coeff[0].dtype == np.double: buf = buf.real
buf = buf.reshape((nocc,nmo,nocc,nmo))
for i in range(nocc):
for p in range(nmo):
if spin[i] != spin[p]:
buf[i,p,:,:] = 0.
buf[:,:,i,p] = 0.
buf1 = buf - buf.transpose(0,3,2,1)
buf1 = buf1.transpose(0,2,1,3)
cput1 = log.timer_debug1('transforming oopq', *cput1)
self.dtype = buf1.dtype
self.oooo[:,:,:,:] = buf1[:,:,:nocc,:nocc]
self.ooov[:,:,:,:] = buf1[:,:,:nocc,nocc:]
self.oovv[:,:,:,:] = buf1[:,:,nocc:,nocc:]
cput1 = logger.process_clock(), logger.perf_counter()
# <ia||pq> = <ia|pq> - <ia|qp> = (ip|aq) - (iq|ap)
buf = ao2mofn(cc._scf.mol, (orbo,so_coeff,orbv,so_coeff), compact=0)
if mo_coeff[0].dtype == np.double: buf = buf.real
buf = buf.reshape((nocc,nmo,nvir,nmo))
for p in range(nmo):
for i in range(nocc):
if spin[i] != spin[p]:
buf[i,p,:,:] = 0.
for a in range(nvir):
if spin[nocc+a] != spin[p]:
buf[:,:,a,p] = 0.
buf1 = buf - buf.transpose(0,3,2,1)
buf1 = buf1.transpose(0,2,1,3)
cput1 = log.timer_debug1('transforming ovpq', *cput1)
self.ovoo[:,:,:,:] = buf1[:,:,:nocc,:nocc]
self.ovov[:,:,:,:] = buf1[:,:,:nocc,nocc:]
self.ovvv[:,:,:,:] = buf1[:,:,nocc:,nocc:]
for a in range(nvir):
orbva = orbv[:,a].reshape(-1,1)
buf = ao2mofn(cc._scf.mol, (orbva,orbv,orbv,orbv), compact=0)
if mo_coeff[0].dtype == np.double: buf = buf.real
buf = buf.reshape((1,nvir,nvir,nvir))
for b in range(nvir):
if spin[nocc+a] != spin[nocc+b]:
buf[0,b,:,:] = 0.
for c in range(nvir):
if spin[nocc+b] != spin[nocc+c]:
buf[:,:,b,c] = buf[:,:,c,b] = 0.
buf1 = buf - buf.transpose(0,3,2,1)
buf1 = buf1.transpose(0,2,1,3)
self.vvvv[a] = buf1[:]
cput1 = log.timer_debug1('transforming vvvv', *cput1)
log.timer('CCSD integral transformation', *cput0)
[docs]
def uspatial2spin(cc, moidx, mo_coeff):
'''Convert the results of an unrestricted mean-field calculation to spin-orbital form.
Spin-orbital ordering is determined by orbital energy without regard for spin.
Returns:
fock : (nso,nso) ndarray
The Fock matrix in the basis of spin-orbitals
so_coeff : (nao, nso) ndarray
The matrix of spin-orbital coefficients in the AO basis
spin : (nso,) ndarary
The spin (0 or 1) of each spin-orbital
'''
dm = cc._scf.make_rdm1(cc.mo_coeff, cc.mo_occ)
fockao = cc._scf.get_hcore() + cc._scf.get_veff(cc.mol, dm)
fockab = []
for a in range(2):
fockab.append( reduce(numpy.dot, (mo_coeff[a].T, fockao[a], mo_coeff[a])) )
nocc = cc.nocc
nao = cc.mo_coeff[0].shape[0]
nmo = cc.nmo
so_coeff = np.zeros((nao, nmo), dtype=mo_coeff[0].dtype)
nocc_a = int(sum(cc.mo_occ[0]*moidx[0]))
nocc_b = int(sum(cc.mo_occ[1]*moidx[1]))
nmo_a = fockab[0].shape[0]
nmo_b = fockab[1].shape[0]
nvir_a = nmo_a - nocc_a
#nvir_b = nmo_b - nocc_b
oa = range(0,nocc_a)
ob = range(nocc_a,nocc)
va = range(nocc,nocc+nvir_a)
vb = range(nocc+nvir_a,nmo)
spin = np.zeros(nmo, dtype=int)
spin[oa] = 0
spin[ob] = 1
spin[va] = 0
spin[vb] = 1
so_coeff[:,oa] = mo_coeff[0][:,:nocc_a]
so_coeff[:,ob] = mo_coeff[1][:,:nocc_b]
so_coeff[:,va] = mo_coeff[0][:,nocc_a:nmo_a]
so_coeff[:,vb] = mo_coeff[1][:,nocc_b:nmo_b]
fock = np.zeros((nmo, nmo), dtype=fockab[0].dtype)
fock[np.ix_(oa,oa)] = fockab[0][:nocc_a,:nocc_a]
fock[np.ix_(oa,va)] = fockab[0][:nocc_a,nocc_a:]
fock[np.ix_(va,oa)] = fockab[0][nocc_a:,:nocc_a]
fock[np.ix_(va,va)] = fockab[0][nocc_a:,nocc_a:]
fock[np.ix_(ob,ob)] = fockab[1][:nocc_b,:nocc_b]
fock[np.ix_(ob,vb)] = fockab[1][:nocc_b,nocc_b:]
fock[np.ix_(vb,ob)] = fockab[1][nocc_b:,:nocc_b]
fock[np.ix_(vb,vb)] = fockab[1][nocc_b:,nocc_b:]
# Do not sort because it's different to the orbital ordering generated by
# get_frozen_mask function in AO-direct vvvv contraction
# idxo = np.diagonal(fock[:nocc,:nocc]).argsort()
# idxv = nocc + np.diagonal(fock[nocc:,nocc:]).argsort()
# idx = np.concatenate((idxo,idxv))
# spin = spin[idx]
# so_coeff = so_coeff[:,idx]
# fock = fock[:, idx][idx]
return fock, so_coeff, spin
class _IMDS:
# Exactly the same as RCCSD IMDS except
# -- rintermediates --> uintermediates
# -- Loo, Lvv, cc_Fov --> Foo, Fvv, Fov
# -- One less 2-virtual intermediate
def __init__(self, cc):
self.cc = cc
self._made_shared = False
self.made_ip_imds = False
self.made_ea_imds = False
self.made_ee_imds = False
def _make_shared(self):
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.cc.stdout, self.cc.verbose)
t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris
self.Foo = imd.Foo(t1,t2,eris)
self.Fvv = imd.Fvv(t1,t2,eris)
self.Fov = imd.Fov(t1,t2,eris)
# 2 virtuals
self.Wovvo = imd.Wovvo(t1,t2,eris)
self.Woovv = eris.oovv
log.timer('EOM-CCSD shared intermediates', *cput0)
def make_ip(self):
if self._made_shared is False:
self._make_shared()
self._made_shared = True
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.cc.stdout, self.cc.verbose)
t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris
# 0 or 1 virtuals
self.Woooo = imd.Woooo(t1,t2,eris)
self.Wooov = imd.Wooov(t1,t2,eris)
self.Wovoo = imd.Wovoo(t1,t2,eris)
self.made_ip_imds = True
log.timer('EOM-CCSD IP intermediates', *cput0)
def make_ea(self):
if self._made_shared is False:
self._make_shared()
self._made_shared = True
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.cc.stdout, self.cc.verbose)
t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris
# 3 or 4 virtuals
self.Wvovv = imd.Wvovv(t1,t2,eris)
self.Wvvvv = imd.Wvvvv(t1,t2,eris)
self.Wvvvo = imd.Wvvvo(t1,t2,eris,self.Wvvvv)
self.made_ea_imds = True
log.timer('EOM-CCSD EA intermediates', *cput0)
def make_ee(self):
if self._made_shared is False:
self._make_shared()
self._made_shared = True
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.cc.stdout, self.cc.verbose)
t1,t2,eris = self.cc.t1, self.cc.t2, self.cc.eris
if self.made_ip_imds is False:
# 0 or 1 virtuals
self.Woooo = imd.Woooo(t1,t2,eris)
self.Wooov = imd.Wooov(t1,t2,eris)
self.Wovoo = imd.Wovoo(t1,t2,eris)
if self.made_ea_imds is False:
# 3 or 4 virtuals
self.Wvovv = imd.Wvovv(t1,t2,eris)
self.Wvvvv = imd.Wvvvv(t1,t2,eris)
self.Wvvvo = imd.Wvvvo(t1,t2,eris,self.Wvvvv)
self.made_ee_imds = True
log.timer('EOM-CCSD EE intermediates', *cput0)
if __name__ == '__main__':
from pyscf import scf
from pyscf import gto
mol = gto.Mole()
mol.atom = [['O', (0., 0., 0.)],
['O', (1.21, 0., 0.)]]
mol.basis = 'cc-pvdz'
mol.spin = 2
mol.build()
mf = scf.UHF(mol)
print(mf.scf())
# Freeze 1s electrons
frozen = [[0,1], [0,1]]
# also acceptable
#frozen = 4
ucc = UCCSD(mf, frozen=frozen)
ecc, t1, t2 = ucc.kernel()
print(ecc - -0.3486987472235819)
exit()
mol = gto.Mole()
mol.atom = [
[8 , (0. , 0. , 0.)],
[1 , (0. , -0.757 , 0.587)],
[1 , (0. , 0.757 , 0.587)]]
mol.basis = 'cc-pvdz'
mol.spin = 0
mol.build()
mf = scf.UHF(mol)
print(mf.scf())
mycc = UCCSD(mf)
ecc, t1, t2 = mycc.kernel()
print(ecc - -0.2133432712431435)
e,v = mycc.ipccsd(nroots=8)
print(e[0] - 0.4335604332073799)
print(e[2] - 0.5187659896045407)
print(e[4] - 0.6782876002229172)
mycc.verbose = 5
e,v = mycc.eaccsd(nroots=8)
print(e[0] - 0.16737886338859731)
print(e[2] - 0.24027613852009164)
print(e[4] - 0.51006797826488071)
print("e=", e)
e,v = mycc.eeccsd(nroots=4)
print(e[0] - 0.2757159395886167)
print(e[1] - 0.2757159395886167)
print(e[2] - 0.2757159395886167)
print(e[3] - 0.3005716731825082)
