#!/usr/bin/env python
# Copyright 2014-2019 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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
# James D. McClain
# Timothy Berkelbach <tim.berkelbach@gmail.com>
#
import numpy as np
from pyscf import lib
from pyscf import ao2mo
from pyscf.lib import logger, module_method
from pyscf.cc import ccsd
from pyscf.cc import rintermediates as imd
from pyscf import __config__
[docs]
def kernel(eom, nroots=1, koopmans=False, guess=None, left=False,
eris=None, imds=None, **kwargs):
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(eom.stdout, eom.verbose)
if eom.verbose >= logger.WARN:
eom.check_sanity()
eom.dump_flags()
if imds is None:
imds = eom.make_imds(eris)
matvec, diag = eom.gen_matvec(imds, left=left, **kwargs)
size = eom.vector_size()
nroots = min(nroots, size)
if guess is not None:
user_guess = True
for g in guess:
assert g.size == size
else:
user_guess = False
guess = eom.get_init_guess(nroots, koopmans, diag)
def precond(r, e0, x0):
return r/(e0-diag+1e-12)
# GHF or customized RHF/UHF may be of complex type
real_system = (eom._cc._scf.mo_coeff[0].dtype == np.double)
eig = lib.davidson_nosym1
if user_guess or koopmans:
assert len(guess) == nroots
def eig_close_to_init_guess(w, v, nroots, envs):
x0 = lib.linalg_helper._gen_x0(envs['v'], envs['xs'])
s = np.dot(np.asarray(guess).conj(), np.asarray(x0).T)
snorm = np.einsum('pi,pi->i', s.conj(), s)
idx = np.argsort(-snorm)[:nroots]
return lib.linalg_helper._eigs_cmplx2real(w, v, idx, real_system)
conv, es, vs = eig(matvec, guess, precond, pick=eig_close_to_init_guess,
tol=eom.conv_tol, max_cycle=eom.max_cycle,
max_space=eom.max_space, nroots=nroots, verbose=log)
else:
def pickeig(w, v, nroots, envs):
real_idx = np.where(abs(w.imag) < 1e-3)[0]
return lib.linalg_helper._eigs_cmplx2real(w, v, real_idx, real_system)
conv, es, vs = eig(matvec, guess, precond, pick=pickeig,
tol=eom.conv_tol, max_cycle=eom.max_cycle,
max_space=eom.max_space, nroots=nroots, verbose=log)
if eom.verbose >= logger.INFO:
for n, en, vn, convn in zip(range(nroots), es, vs, conv):
r1, r2 = eom.vector_to_amplitudes(vn)
if isinstance(r1, np.ndarray):
qp_weight = np.linalg.norm(r1)**2
else: # for EOM-UCCSD
r1 = np.hstack([x.ravel() for x in r1])
qp_weight = np.linalg.norm(r1)**2
logger.info(eom, 'EOM-CCSD root %d E = %.16g qpwt = %.6g conv = %s',
n, en, qp_weight, convn)
log.timer('EOM-CCSD', *cput0)
if nroots == 1:
return conv[0], es[0].real, vs[0]
else:
return conv, es.real, vs
[docs]
class EOM(lib.StreamObject):
_keys = {
'mol', 'max_space', 'max_cycle', 'conv_tol', 'partition',
'e', 'v', 'nocc', 'nmo',
}
def __init__(self, cc):
self.mol = cc.mol
self._cc = cc
self.verbose = cc.verbose
self.stdout = cc.stdout
self.max_memory = cc.max_memory
self.max_space = getattr(__config__, 'eom_rccsd_EOM_max_space', 20)
self.max_cycle = getattr(__config__, 'eom_rccsd_EOM_max_cycle', cc.max_cycle)
self.conv_tol = getattr(__config__, 'eom_rccsd_EOM_conv_tol', cc.conv_tol)
self.partition = getattr(__config__, 'eom_rccsd_EOM_partition', None)
##################################################
# don't modify the following attributes, they are not input options
self.e = None
self.v = None
self.nocc = cc.nocc
self.nmo = cc.nmo
[docs]
def dump_flags(self, verbose=None):
logger.info(self, '')
logger.info(self, '******** %s ********', self.__class__)
logger.info(self, 'max_space = %d', self.max_space)
logger.info(self, 'max_cycle = %d', self.max_cycle)
logger.info(self, 'conv_tol = %s', self.conv_tol)
logger.info(self, 'partition = %s', self.partition)
#logger.info(self, 'nocc = %d', self.nocc)
#logger.info(self, 'nmo = %d', self.nmo)
logger.info(self, 'max_memory %d MB (current use %d MB)',
self.max_memory, lib.current_memory()[0])
return self
[docs]
def reset(self, mol=None):
self._cc.reset(mol)
return self
def _sort_left_right_eigensystem(eom, right_converged, right_evals, right_evecs,
left_converged, left_evals, left_evecs, tol=1e-6):
'''Ensures the left and right eigenvectors correspond to the same eigenvalue.
Note:
Useful for perturbative methods that need both eigenstates. Right now, just
simply checks for equality between left and right eigenvalues, but can be
extended to make sure the overlap between states is sufficiently large.
Kwargs:
eom : :class:`EOM`
Class holding EOM results.
right_converged : array-like of bool
Whether the right eigenstates converged.
right_evals : array-like
Eigenvalues of right eigenstates.
right_evecs : array-like of ndarray
Eigenvectors of right eigenstates.
left_converged : array-like of bool
Whether the left eigenstates converged.
left_evals : array-like
Eigenvalues of left eigenstates.
left_evecs : array-like of ndarray
Eigenvectors of left eigenstates.
tol : float
Tolerance for determining whether a left and right eigenvalue
should be considered equal.
'''
log = logger.Logger(eom.stdout, eom.verbose)
right_evecs, left_evecs = [np.atleast_2d(x) for x in [right_evecs, left_evecs]]
right_evals, left_evals = [np.atleast_1d(x) for x in [right_evals, left_evals]]
right_converged, left_converged = [np.atleast_1d(x) for x in [right_converged, left_converged]]
srt_right_idx = []
srt_left_idx = []
left_idx = [idx for idx in range(len(left_evals)) if left_converged[idx]]
right_idx = [idx for idx in range(len(right_evals)) if right_converged[idx]]
if len(right_idx) != len(left_idx):
log.warn('Number of converged left and right eigenvalues are not equal.\n'
' No. Left = %3d, No. Right = %3d.' %
(len(left_idx), len(right_idx)))
for ir_idx, ir in enumerate(right_idx):
found = False
for il_idx, il in enumerate(left_idx):
if abs(right_evals[ir] - left_evals[il]) < tol:
found = True
srt_right_idx.append(ir)
srt_left_idx.append(il)
break
if found:
left_idx.pop(il_idx)
else:
log.warn('No converged left eigenvalue corresponding to right eigenvalue '
'%.6g (right idx=%3d).\nWill not perform perturbation on this state.'
% (right_evals[ir], ir))
log.info('Resulting left/right eigenstates:')
log.info('Left Eigen (idx) <-> Right Eigen (idx)')
for il, ir in zip(srt_left_idx, srt_right_idx):
log.info('%10.6g (%3d) %10.6g (%3d)',
left_evals[il], il, right_evals[ir], ir)
return (right_evals[srt_right_idx], right_evecs[srt_right_idx], left_evecs[srt_left_idx])
[docs]
def perturbed_ccsd_kernel(eom, nroots=1, koopmans=False, right_guess=None,
left_guess=None, eris=None, imds=None):
'''Wrapper for running perturbative excited-states that require both left
and right amplitudes.'''
if imds is None:
imds = eom.make_imds(eris=eris)
# Right eigenvectors
r_converged, r_e, r_v = \
kernel(eom, nroots, koopmans=koopmans, guess=right_guess, left=False,
eris=eris, imds=imds)
# Left eigenvectors
l_converged, l_e, l_v = \
kernel(eom, nroots, koopmans=koopmans, guess=right_guess, left=True,
eris=eris, imds=imds)
e, r_v, l_v = _sort_left_right_eigensystem(eom, r_converged, r_e, r_v, l_converged, l_e, l_v)
e_star = eom.ccsd_star_contract(e, r_v, l_v, imds=imds)
return e_star
########################################
# EOM-IP-CCSD
########################################
[docs]
def ipccsd(eom, nroots=1, left=False, koopmans=False, guess=None,
partition=None, eris=None, imds=None):
'''Calculate (N-1)-electron charged excitations via IP-EOM-CCSD.
Kwargs:
nroots : int
Number of roots (eigenvalues) requested
partition : bool or str
Use a matrix-partitioning for the doubles-doubles block.
Can be None, 'mp' (Moller-Plesset, i.e. orbital energies on the diagonal),
or 'full' (full diagonal elements).
koopmans : bool
Calculate Koopmans'-like (quasiparticle) excitations only, targeting via
overlap.
guess : list of ndarray
List of guess vectors to use for targeting via overlap.
'''
if partition is not None:
eom.partition = partition.lower()
assert eom.partition in ['mp','full']
eom.converged, eom.e, eom.v \
= kernel(eom, nroots, koopmans, guess, left, eris=eris, imds=imds)
return eom.e, eom.v
[docs]
def ipccsd_star(eom, nroots=1, koopmans=False, right_guess=None,
left_guess=None, eris=None, imds=None):
"""Calculates CCSD* perturbative correction.
Simply calls the relevant `kernel()` function and `perturb_star` of the
`eom` class.
Returns:
e_t_a_star (list of float):
The IP-CCSD* energy.
"""
return perturbed_ccsd_kernel(eom, nroots=nroots, koopmans=koopmans,
right_guess=right_guess, left_guess=left_guess, eris=eris,
imds=imds)
[docs]
def vector_to_amplitudes_ip(vector, nmo, nocc):
nvir = nmo - nocc
r1 = vector[:nocc].copy()
r2 = vector[nocc:].copy().reshape(nocc,nocc,nvir)
return r1, r2
[docs]
def amplitudes_to_vector_ip(r1, r2):
vector = np.hstack((r1, r2.ravel()))
return vector
[docs]
def spatial2spin_ip(rx, orbspin=None):
from pyscf.cc import eom_uccsd
if rx.ndim == 1:
r1 = (rx, np.zeros_like(rx))
return eom_uccsd.spatial2spin_ip(r1, orbspin)
else:
rx0 = np.zeros_like(rx)
r2 = (rx.transpose(1,0,2)-rx, rx0, -rx, rx0)
return eom_uccsd.spatial2spin_ip(r2, orbspin)
[docs]
def ipccsd_matvec(eom, vector, imds=None, diag=None):
# Ref: Nooijen and Snijders, J. Chem. Phys. 102, 1681 (1995) Eqs.(8)-(9)
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
# 1h-1h block
Hr1 = -np.einsum('ki,k->i', imds.Loo, r1)
#1h-2h1p block
Hr1 += 2*np.einsum('ld,ild->i', imds.Fov, r2)
Hr1 += -np.einsum('kd,kid->i', imds.Fov, r2)
Hr1 += -2*np.einsum('klid,kld->i', imds.Wooov, r2)
Hr1 += np.einsum('lkid,kld->i', imds.Wooov, r2)
# 2h1p-1h block
Hr2 = -np.einsum('kbij,k->ijb', imds.Wovoo, r1)
# 2h1p-2h1p block
if eom.partition == 'mp':
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr2 += lib.einsum('bd,ijd->ijb', fvv, r2)
Hr2 += -lib.einsum('ki,kjb->ijb', foo, r2)
Hr2 += -lib.einsum('lj,ilb->ijb', foo, r2)
elif eom.partition == 'full':
diag_matrix2 = eom.vector_to_amplitudes(diag, nmo, nocc)[1]
Hr2 += diag_matrix2 * r2
else:
Hr2 += lib.einsum('bd,ijd->ijb', imds.Lvv, r2)
Hr2 += -lib.einsum('ki,kjb->ijb', imds.Loo, r2)
Hr2 += -lib.einsum('lj,ilb->ijb', imds.Loo, r2)
Hr2 += lib.einsum('klij,klb->ijb', imds.Woooo, r2)
Hr2 += 2*lib.einsum('lbdj,ild->ijb', imds.Wovvo, r2)
Hr2 += -lib.einsum('kbdj,kid->ijb', imds.Wovvo, r2)
Hr2 += -lib.einsum('lbjd,ild->ijb', imds.Wovov, r2) #typo in Ref
Hr2 += -lib.einsum('kbid,kjd->ijb', imds.Wovov, r2)
tmp = 2*np.einsum('lkdc,kld->c', imds.Woovv, r2)
tmp += -np.einsum('kldc,kld->c', imds.Woovv, r2)
Hr2 += -np.einsum('c,ijcb->ijb', tmp, imds.t2)
vector = eom.amplitudes_to_vector(Hr1, Hr2)
return vector
[docs]
def lipccsd_matvec(eom, vector, imds=None, diag=None):
'''For left eigenvector'''
# Note this is not the same left EA equations used by Nooijen and Bartlett.
# Small changes were made so that the same type L2 basis was used for both the
# left EA and left IP equations. You will note more similarity for these
# equations to the left IP equations than for the left EA equations by Nooijen.
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
# 1h-1h block
Hr1 = -np.einsum('ki,i->k', imds.Loo, r1)
#1h-2h1p block
Hr1 += -np.einsum('kbij,ijb->k', imds.Wovoo, r2)
# 2h1p-1h block
Hr2 = -np.einsum('kd,l->kld', imds.Fov, r1)
Hr2 += 2.*np.einsum('ld,k->kld', imds.Fov, r1)
Hr2 += -np.einsum('klid,i->kld', 2.*imds.Wooov-imds.Wooov.transpose(1,0,2,3), r1)
# 2h1p-2h1p block
if eom.partition == 'mp':
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr2 += lib.einsum('bd,klb->kld', fvv, r2)
Hr2 += -lib.einsum('ki,ild->kld', foo, r2)
Hr2 += -lib.einsum('lj,kjd->kld', foo, r2)
elif eom.partition == 'full':
diag_matrix2 = eom.vector_to_amplitudes(diag, nmo, nocc)[1]
Hr2 += diag_matrix2 * r2
else:
Hr2 += lib.einsum('bd,klb->kld', imds.Lvv, r2)
Hr2 += -lib.einsum('ki,ild->kld', imds.Loo, r2)
Hr2 += -lib.einsum('lj,kjd->kld', imds.Loo, r2)
Hr2 += lib.einsum('lbdj,kjb->kld', 2.*imds.Wovvo-imds.Wovov.transpose(0,1,3,2), r2)
Hr2 += -lib.einsum('kbdj,ljb->kld', imds.Wovvo, r2)
Hr2 += lib.einsum('klij,ijd->kld', imds.Woooo, r2)
Hr2 += -lib.einsum('kbid,ilb->kld', imds.Wovov, r2)
tmp = np.einsum('ijcb,ijb->c', imds.t2, r2)
Hr2 += -np.einsum('lkdc,c->kld', 2.*imds.Woovv-imds.Woovv.transpose(1,0,2,3), tmp)
vector = eom.amplitudes_to_vector(Hr1, Hr2)
return vector
[docs]
def ipccsd_diag(eom, imds=None):
if imds is None: imds = eom.make_imds()
t1, t2 = imds.t1, imds.t2
dtype = np.result_type(t1, t2)
nocc, nvir = t1.shape
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr1 = -np.diag(imds.Loo)
Hr2 = np.zeros((nocc,nocc,nvir), dtype)
for i in range(nocc):
for j in range(nocc):
for b in range(nvir):
if eom.partition == 'mp':
Hr2[i,j,b] += fvv[b,b]
Hr2[i,j,b] += -foo[i,i]
Hr2[i,j,b] += -foo[j,j]
else:
Hr2[i,j,b] += imds.Lvv[b,b]
Hr2[i,j,b] += -imds.Loo[i,i]
Hr2[i,j,b] += -imds.Loo[j,j]
Hr2[i,j,b] += imds.Woooo[i,j,i,j]
Hr2[i,j,b] +=2*imds.Wovvo[j,b,b,j]
Hr2[i,j,b] += -imds.Wovvo[i,b,b,i]*(i==j)
Hr2[i,j,b] += -imds.Wovov[j,b,j,b]
Hr2[i,j,b] += -imds.Wovov[i,b,i,b]
Hr2[i,j,b] += -2*np.dot(imds.Woovv[j,i,b,:], t2[i,j,:,b])
Hr2[i,j,b] += np.dot(imds.Woovv[i,j,b,:], t2[i,j,:,b])
vector = eom.amplitudes_to_vector(Hr1, Hr2)
return vector
[docs]
def ipccsd_star_contract(eom, ipccsd_evals, ipccsd_evecs, lipccsd_evecs, imds=None):
from pyscf.cc.ccsd_t import _sort_eri, _sort_t2_vooo_
cpu1 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(eom.stdout, eom.verbose)
if imds is None:
imds = eom.make_imds()
t1, t2 = imds.t1, imds.t2
eris = imds.eris
nocc, nvir = t1.shape
nmo = nocc + nvir
dtype = np.result_type(t1, t2, eris.ovoo.dtype)
if eom._cc.incore_complete:
ftmp = None
eris_vvop = np.zeros((nvir,nvir,nocc,nmo), dtype)
else:
ftmp = lib.H5TmpFile()
eris_vvop = ftmp.create_dataset('vvop', (nvir,nvir,nocc,nmo), dtype)
orbsym = _sort_eri(eom._cc, eris, nocc, nvir, eris_vvop, log)
mo_energy, t1T, t2T, vooo, fvo, restore_t2_inplace = \
_sort_t2_vooo_(eom._cc, orbsym, t1, t2, eris)
cpu1 = log.timer_debug1('CCSD(T) sort_eri', *cpu1)
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now)
blksize = min(nvir, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nocc**3*6))))
mo_e_occ = np.asarray(mo_energy[:nocc])
mo_e_vir = np.asarray(mo_energy[nocc:])
def contract_l2p(l1, l2, a0, a1, b0, b1, cache_vvop, out=None):
'''Create perturbed l2.'''
if out is None:
out = np.zeros((nocc,)*3 + (a1-a0,b1-b0), dtype=dtype)
out += 0.5*np.einsum('abij,k->ijkab', cache_vvop[:,:,:,:nocc].conj(), l1)
out += lib.einsum('abie,jke->ijkab', cache_vvop[:,:,:,nocc:].conj(), l2)
out += -lib.einsum('bjkm,ima->ijkab', vooo[b0:b1], l2[:,:,a0:a1])
out += -lib.einsum('bjim,mka->ijkab', vooo[b0:b1], l2[:,:,a0:a1])
return out
def contract_pl2p(l1, l2, a0, a1, b0, b1, cache_vvop_a, cache_vvop_b):
'''Create P(ia|jb) of perturbed l2.'''
out = contract_l2p(l1, l2, a0, a1, b0, b1, cache_vvop_a)
out += contract_l2p(l1, l2, b0, b1, a0, a1, cache_vvop_b).transpose(1,0,2,4,3) # P(ia|jb)
return out
def contract_r2p(r1, r2, a0, a1, b0, b1, cache_vvop, out=None):
'''Create perturbed r2.'''
if out is None:
out = np.zeros((nocc,)*3 + (a1-a0,b1-b0), dtype=dtype)
tmp = np.einsum('mkbe,m->bke', eris.oovv[:,:,b0:b1,:], r1)
out += -lib.einsum('bke,aeji->ijkab', tmp, t2T[a0:a1])
tmp = np.einsum('mebj,m->bej', eris.ovvo[:,:,b0:b1,:], r1)
out += -lib.einsum('bej,aeki->ijkab', tmp, t2T[a0:a1])
tmp = np.einsum('mjnk,n->mjk', eris.oooo, r1)
out += lib.einsum('mjk,abmi->ijkab', tmp, t2T[a0:a1,b0:b1])
out += lib.einsum('abie,kje->ijkab', cache_vvop[:,:,:,nocc:], r2)
out += -lib.einsum('bjkm,mia->ijkab', vooo[b0:b1].conj(), r2[:,:,a0:a1])
out += -lib.einsum('bjim,kma->ijkab', vooo[b0:b1].conj(), r2[:,:,a0:a1])
return out
def contract_pr2p(r1, r2, a0, a1, b0, b1, cache_vvop_a, cache_vvop_b):
'''Create P(ia|jb) of perturbed r2.'''
out = contract_r2p(r1, r2, a0, a1, b0, b1, cache_vvop_a)
out += contract_r2p(r1, r2, b0, b1, a0, a1, cache_vvop_b).transpose(1,0,2,4,3) # P(ia|jb)
return out
ipccsd_evecs = np.array(ipccsd_evecs)
lipccsd_evecs = np.array(lipccsd_evecs)
e = []
ipccsd_evecs, lipccsd_evecs = [np.atleast_2d(x) for x in [ipccsd_evecs, lipccsd_evecs]]
ipccsd_evals = np.atleast_1d(ipccsd_evals)
for eval_, evec_, levec_ in zip(ipccsd_evals, ipccsd_evecs, lipccsd_evecs):
l1, l2 = eom.vector_to_amplitudes(levec_)
r1, r2 = eom.vector_to_amplitudes(evec_)
ldotr = np.dot(l1, r1) + np.dot(l2.ravel(), r2.ravel())
l1 /= ldotr
l2 /= ldotr
l2 = 1./3*(l2 + 2.*l2.transpose(1,0,2))
deltaE = 0.0
eijk = (mo_e_occ[:,None,None,None,None] +
mo_e_occ[None,:,None,None,None] +
mo_e_occ[None,None,:,None,None] + eval_)
for a0, a1 in lib.prange_tril(0, nvir, blksize):
b0, b1 = 0, a1
eijkab = (eijk - mo_e_vir[a0:a1][None,None,None,:,None] -
mo_e_vir[b0:b1][None,None,None,None,:])
eijkab = 1./eijkab
vvov_a = eris_vvop[a0:a1,b0:b1,:,:]
vvov_b = eris_vvop[b0:b1,a0:a1,:,:]
lijkab = contract_pl2p(l1, l2, a0, a1, b0, b1, vvov_a, vvov_b)
rijkab = contract_pr2p(r1, r2, a0, a1, b0, b1, vvov_a, vvov_b)
lijkab = 4.*lijkab \
- 2.*lijkab.transpose(1,0,2,3,4) \
- 2.*lijkab.transpose(2,1,0,3,4) \
- 2.*lijkab.transpose(0,2,1,3,4) \
+ 1.*lijkab.transpose(1,2,0,3,4) \
+ 1.*lijkab.transpose(2,0,1,3,4)
# Symmetry factors (1 for a == b, 2 for a < b)
fac = 2*np.ones_like(rijkab, dtype=int)
triu_idx = np.triu_indices(a1-a0, a0+1, m=b1-b0)
fac[:,:,:,triu_idx[0],triu_idx[1]] = 0
fac[:,:,:,np.arange(a1-a0),np.arange(a0,b1)] = 1
eijkab *= fac
deltaE += np.einsum('ijkab,ijkab,ijkab', lijkab, rijkab, eijkab)
deltaE = 0.5*deltaE.real
logger.info(eom, "ipccsd energy, star energy, delta energy = %16.12f, %16.12f, %16.12f",
eval_, eval_+deltaE, deltaE)
e.append(eval_+deltaE)
t2 = restore_t2_inplace(t2T)
return e
[docs]
class EOMIP(EOM):
[docs]
def get_init_guess(self, nroots=1, koopmans=True, diag=None):
size = self.vector_size()
nroots = min(nroots, size)
guess = []
if koopmans:
dtype = getattr(diag, 'dtype', np.double)
for n in range(nroots):
g = np.zeros(int(size), dtype)
g[self.nocc-n-1] = 1.0
guess.append(g)
else:
if diag is None:
diag = self.get_diag()
dtype = getattr(diag, 'dtype', np.double)
idx = diag.argsort()[:nroots]
for i in idx:
g = np.zeros(int(size), dtype)
g[i] = 1.0
guess.append(g)
return guess
kernel = ipccsd
ipccsd = ipccsd
ipccsd_star = ipccsd_star
matvec = ipccsd_matvec
l_matvec = lipccsd_matvec
get_diag = ipccsd_diag
ccsd_star_contract = ipccsd_star_contract
[docs]
def ipccsd_star_contract(self, ipccsd_evals, ipccsd_evecs, lipccsd_evecs, imds=None):
return self.ccsd_star_contract(ipccsd_evals, ipccsd_evecs, lipccsd_evecs, imds=imds)
[docs]
def gen_matvec(self, imds=None, left=False, **kwargs):
if imds is None: imds = self.make_imds()
diag = self.get_diag(imds)
if left:
matvec = lambda xs: [self.l_matvec(x, imds, diag) for x in xs]
else:
matvec = lambda xs: [self.matvec(x, imds, diag) for x in xs]
return matvec, diag
amplitudes_to_vector = staticmethod(amplitudes_to_vector_ip)
vector_to_amplitudes = module_method(vector_to_amplitudes_ip,
absences=['nmo', 'nocc'])
spatial2spin = staticmethod(spatial2spin_ip)
[docs]
def vector_size(self):
nocc = self.nocc
nvir = self.nmo - nocc
return nocc + nocc*nocc*nvir
[docs]
def make_imds(self, eris=None):
imds = _IMDS(self._cc, eris=eris)
imds.make_ip(self.partition)
return imds
@property
def eip(self):
return self.e
[docs]
class EOMIP_Ta(EOMIP):
'''Class for EOM IPCCSD(T)*(a) method by Matthews and Stanton.'''
[docs]
def make_imds(self, eris=None):
imds = _IMDS(self._cc, eris=eris)
imds.make_t3p2_ip(self._cc, self.partition)
return imds
########################################
# EOM-EA-CCSD
########################################
[docs]
def eaccsd(eom, nroots=1, left=False, koopmans=False, guess=None,
partition=None, eris=None, imds=None):
'''Calculate (N+1)-electron charged excitations via EA-EOM-CCSD.
Args:
See also ipccd()
'''
return ipccsd(eom, nroots, left, koopmans, guess, partition, eris, imds)
[docs]
def eaccsd_star(eom, nroots=1, koopmans=False, right_guess=None,
left_guess=None, eris=None, imds=None, **kwargs):
"""Calculates CCSD* perturbative correction.
Args:
See also ipccd_star()
"""
return perturbed_ccsd_kernel(eom, nroots=nroots, koopmans=koopmans,
right_guess=right_guess, left_guess=left_guess, eris=eris,
imds=imds)
[docs]
def vector_to_amplitudes_ea(vector, nmo, nocc):
nvir = nmo - nocc
r1 = vector[:nvir].copy()
r2 = vector[nvir:].copy().reshape(nocc,nvir,nvir)
return r1, r2
[docs]
def amplitudes_to_vector_ea(r1, r2):
vector = np.hstack((r1, r2.ravel()))
return vector
[docs]
def spatial2spin_ea(rx, orbspin=None):
from pyscf.cc import eom_uccsd
if rx.ndim == 1:
r1 = (rx, np.zeros_like(rx))
return eom_uccsd.spatial2spin_ea(r1, orbspin)
else:
rx0 = np.zeros_like(rx)
r2 = (rx-rx.transpose(0,2,1), rx0, rx, rx0)
return eom_uccsd.spatial2spin_ea(r2, orbspin)
[docs]
def eaccsd_matvec(eom, vector, imds=None, diag=None):
# Ref: Nooijen and Bartlett, J. Chem. Phys. 102, 3629 (1995) Eqs.(30)-(31)
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
nvir = nmo - nocc
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
# Eq. (37)
# 1p-1p block
Hr1 = np.einsum('ac,c->a', imds.Lvv, r1)
# 1p-2p1h block
Hr1 += np.einsum('ld,lad->a', 2.*imds.Fov, r2)
Hr1 += np.einsum('ld,lda->a', -imds.Fov, r2)
Hr1 += np.einsum('alcd,lcd->a', 2.*imds.Wvovv-imds.Wvovv.transpose(0,1,3,2), r2)
# Eq. (38)
# 2p1h-1p block
Hr2 = np.einsum('abcj,c->jab', imds.Wvvvo, r1)
# 2p1h-2p1h block
if eom.partition == 'mp':
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr2 += lib.einsum('ac,jcb->jab', fvv, r2)
Hr2 += lib.einsum('bd,jad->jab', fvv, r2)
Hr2 += -lib.einsum('lj,lab->jab', foo, r2)
elif eom.partition == 'full':
diag_matrix2 = eom.vector_to_amplitudes(diag, nmo, nocc)[1]
Hr2 += diag_matrix2 * r2
else:
Hr2 += lib.einsum('ac,jcb->jab', imds.Lvv, r2)
Hr2 += lib.einsum('bd,jad->jab', imds.Lvv, r2)
Hr2 += -lib.einsum('lj,lab->jab', imds.Loo, r2)
Hr2 += lib.einsum('lbdj,lad->jab', 2.*imds.Wovvo-imds.Wovov.transpose(0,1,3,2), r2)
Hr2 += -lib.einsum('lajc,lcb->jab', imds.Wovov, r2)
Hr2 += -lib.einsum('lbcj,lca->jab', imds.Wovvo, r2)
for a in range(nvir):
Hr2[:,a,:] += lib.einsum('bcd,jcd->jb', imds.Wvvvv[a], r2)
tmp = np.einsum('klcd,lcd->k', 2.*imds.Woovv-imds.Woovv.transpose(0,1,3,2), r2)
Hr2 += -np.einsum('k,kjab->jab', tmp, imds.t2)
vector = eom.amplitudes_to_vector(Hr1,Hr2)
return vector
[docs]
def leaccsd_matvec(eom, vector, imds=None, diag=None):
# Note this is not the same left EA equations used by Nooijen and Bartlett.
# Small changes were made so that the same type L2 basis was used for both the
# left EA and left IP equations. You will note more similarity for these
# equations to the left IP equations than for the left EA equations by Nooijen.
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
nvir = nmo - nocc
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
# Eq. (30)
# 1p-1p block
Hr1 = np.einsum('ac,a->c', imds.Lvv, r1)
# 1p-2p1h block
Hr1 += np.einsum('abcj,jab->c', imds.Wvvvo, r2)
# Eq. (31)
# 2p1h-1p block
Hr2 = 2.*np.einsum('c,ld->lcd', r1, imds.Fov)
Hr2 += -np.einsum('d,lc->lcd', r1, imds.Fov)
Hr2 += np.einsum('a,alcd->lcd', r1, 2.*imds.Wvovv-imds.Wvovv.transpose(0,1,3,2))
# 2p1h-2p1h block
if eom.partition == 'mp':
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr2 += lib.einsum('lad,ac->lcd', r2, fvv)
Hr2 += lib.einsum('lcb,bd->lcd', r2, fvv)
Hr2 += -lib.einsum('jcd,lj->lcd', r2, foo)
elif eom.partition == 'full':
diag_matrix2 = eom.vector_to_amplitudes(diag, nmo, nocc)[1]
Hr2 += diag_matrix2 * r2
else:
Hr2 += lib.einsum('lad,ac->lcd', r2, imds.Lvv)
Hr2 += lib.einsum('lcb,bd->lcd', r2, imds.Lvv)
Hr2 += -lib.einsum('jcd,lj->lcd', r2, imds.Loo)
Hr2 += lib.einsum('jcb,lbdj->lcd', r2, 2.*imds.Wovvo-imds.Wovov.transpose(0,1,3,2))
Hr2 += -lib.einsum('lajc,jab->lcb', imds.Wovov, r2)
Hr2 += -lib.einsum('lbcj,jab->lca', imds.Wovvo, r2)
for a in range(nvir):
Hr2 += lib.einsum('lb,bcd->lcd', r2[:,a,:], imds.Wvvvv[a])
tmp = np.einsum('ijcb,ibc->j', imds.t2, r2)
Hr2 += -np.einsum('kjfe,j->kef', 2.*imds.Woovv-imds.Woovv.transpose(0,1,3,2),tmp)
vector = eom.amplitudes_to_vector(Hr1,Hr2)
return vector
[docs]
def eaccsd_diag(eom, imds=None):
if imds is None: imds = eom.make_imds()
t1, t2 = imds.t1, imds.t2
dtype = np.result_type(t1, t2)
nocc, nvir = t1.shape
fock = imds.eris.fock
foo = fock[:nocc,:nocc]
fvv = fock[nocc:,nocc:]
Hr1 = np.diag(imds.Lvv)
Hr2 = np.zeros((nocc,nvir,nvir), dtype)
for a in range(nvir):
if eom.partition != 'mp':
_Wvvvva = np.array(imds.Wvvvv[a])
for b in range(nvir):
for j in range(nocc):
if eom.partition == 'mp':
Hr2[j,a,b] += fvv[a,a]
Hr2[j,a,b] += fvv[b,b]
Hr2[j,a,b] += -foo[j,j]
else:
Hr2[j,a,b] += imds.Lvv[a,a]
Hr2[j,a,b] += imds.Lvv[b,b]
Hr2[j,a,b] += -imds.Loo[j,j]
Hr2[j,a,b] += 2*imds.Wovvo[j,b,b,j]
Hr2[j,a,b] += -imds.Wovov[j,b,j,b]
Hr2[j,a,b] += -imds.Wovov[j,a,j,a]
Hr2[j,a,b] += -imds.Wovvo[j,b,b,j]*(a==b)
Hr2[j,a,b] += _Wvvvva[b,a,b]
Hr2[j,a,b] += -2*np.dot(imds.Woovv[:,j,a,b], t2[:,j,a,b])
Hr2[j,a,b] += np.dot(imds.Woovv[:,j,b,a], t2[:,j,a,b])
vector = eom.amplitudes_to_vector(Hr1,Hr2)
return vector
[docs]
def eaccsd_star_contract(eom, eaccsd_evals, eaccsd_evecs, leaccsd_evecs, imds=None):
cpu1 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(eom.stdout, eom.verbose)
if imds is None:
imds = eom.make_imds()
t1, t2 = imds.t1, imds.t2
eris = imds.eris
nocc, nvir = t1.shape
dtype = np.result_type(t1, t2, eris.ovoo.dtype)
# Notice we do not use `sort_eri` as compared to the eaccsd_star.
# The sort_eri does not produce eri's that are read-in quickly for the current contraction
# scheme. Here, we have that the block loop is over occupied indices whereas in the
# sort_eri it is done over virtual indices (due to the permutation over occupied indices
# in ipccsd_star versus virtual indices in eaccsd_star).
cpu1 = log.timer_debug1('CCSD(T) sort_eri', *cpu1) # Left if new sort_eri implemented
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*6))))
mo_energy = np.asarray(eris.mo_energy)
mo_e_occ = np.asarray(mo_energy[:nocc])
mo_e_vir = np.asarray(mo_energy[nocc:])
def contract_l2p(l1, l2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j, out=None):
'''Create perturbed l2.'''
if out is None:
out = np.zeros((i1-i0,j1-j0) + (nvir,)*3, dtype=dtype)
out += -0.5*np.einsum('iajb,c->ijabc', eris.ovov[i0:i1,:,j0:j1], l1)
out += lib.einsum('iajm,mbc->ijabc', eris.ovoo[i0:i1,:,j0:j1], l2)
out -= lib.einsum('iaeb,jec->ijabc', cache_ovvv_i, l2[j0:j1])
out -= lib.einsum('jbec,iae->ijabc', cache_ovvv_j, l2[i0:i1])
return out
def contract_pl2p(l1, l2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j):
'''Create P(ia|jb) of perturbed l2.'''
out = contract_l2p(l1, l2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j)
tmp = contract_l2p(l1, l2, j0, j1, i0, i1, cache_ovvv_j, cache_ovvv_i)
tmp = tmp.transpose(1,0,3,2,4) # P(ia|jb)
out = out + tmp
return out
def _get_vvvv(eris):
if eris.vvvv is None and getattr(eris, 'vvL', None) is not None: # DF eris
vvL = np.asarray(eris.vvL)
nvir = int(np.sqrt(eris.vvL.shape[0]*2))
return ao2mo.restore(1, lib.dot(vvL, vvL.T), nvir)
elif eris.vvvv.ndim == 2:
nvir = int(np.sqrt(eris.vvvv.shape[0]*2))
return ao2mo.restore(1, np.asarray(eris.vvvv), nvir)
else:
return eris.vvvv
def contract_r2p(r1, r2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j, out=None):
'''Create perturbed r2.'''
if out is None:
out = np.zeros((i1-i0,j1-j0) + (nvir,)*3, dtype=dtype)
tmp = np.einsum('becf,f->bce', _get_vvvv(eris), r1)
out += -lib.einsum('bce,ijae->ijabc', tmp, t2[i0:i1,j0:j1])
tmp = np.einsum('mjce,e->mcj', eris.oovv[:,j0:j1], r1)
out += lib.einsum('mcj,imab->ijabc', tmp, t2[i0:i1])
tmp = np.einsum('jbem,e->mbj', eris.ovvo[j0:j1], r1)
out += lib.einsum('mbj,imac->ijabc', tmp, t2[i0:i1])
out += lib.einsum('iajm,mbc->ijabc', eris.ovoo[i0:i1,:,j0:j1].conj(), r2)
out += -lib.einsum('iaeb,jec->ijabc', cache_ovvv_i.conj(), r2[j0:j1])
out += -lib.einsum('jbec,iae->ijabc', cache_ovvv_j.conj(), r2[i0:i1])
return out
def contract_pr2p(r1, r2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j):
'''Create P(ia|jb) of perturbed r2.'''
out = contract_r2p(r1, r2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j)
tmp = contract_r2p(r1, r2, j0, j1, i0, i1, cache_ovvv_j, cache_ovvv_i)
tmp = tmp.transpose(1,0,3,2,4) # P(ia|jb)
out = out + tmp
return out
eaccsd_evecs = np.array(eaccsd_evecs)
leaccsd_evecs = np.array(leaccsd_evecs)
e = []
eaccsd_evecs, leaccsd_evecs = [np.atleast_2d(x) for x in [eaccsd_evecs, leaccsd_evecs]]
eaccsd_evals = np.atleast_1d(eaccsd_evals)
for eval_, evec_, levec_ in zip(eaccsd_evals, eaccsd_evecs, leaccsd_evecs):
l1, l2 = eom.vector_to_amplitudes(levec_)
r1, r2 = eom.vector_to_amplitudes(evec_)
ldotr = np.dot(l1, r1) + np.dot(l2.ravel(),r2.ravel())
l1 /= ldotr
l2 /= ldotr
l2 = 1./3*(1.*l2 + 2.*l2.transpose(0,2,1))
r2 = r2.transpose(0,2,1)
deltaE = 0.0
eabc = (mo_e_vir[None,None,:,None,None] +
mo_e_vir[None,None,None,:,None] +
mo_e_vir[None,None,None,None,:] - eval_)
for i0, i1 in lib.prange_tril(0, nocc, blksize):
j0, j1 = 0, i1
eijabc = (mo_e_occ[i0:i1][:,None,None,None,None] +
mo_e_occ[j0:j1][None,:,None,None,None] - eabc)
eijabc = 1./eijabc
cache_ovvv_i = eris.get_ovvv(slice(i0,i1))
cache_ovvv_j = eris.get_ovvv(slice(j0,j1))
lijabc = contract_pl2p(l1, l2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j)
rijabc = contract_pr2p(r1, r2, i0, i1, j0, j1, cache_ovvv_i, cache_ovvv_j)
lijabc = 4.*lijabc \
- 2.*lijabc.transpose(0,1,3,2,4) \
- 2.*lijabc.transpose(0,1,4,3,2) \
- 2.*lijabc.transpose(0,1,2,4,3) \
+ 1.*lijabc.transpose(0,1,3,4,2) \
+ 1.*lijabc.transpose(0,1,4,2,3)
# Symmetry factors (1 for a == b, 2 for a < b)
fac = 2*np.ones_like(rijabc, dtype=int)
triu_idx = np.triu_indices(i1-i0,i0+1,m=j1-j0)
fac[triu_idx[0],triu_idx[1],:,:,:] = 0
fac[np.arange(i1-i0),np.arange(i0,j1)] = 1
eijabc *= fac
deltaE += np.einsum('ijabc,ijabc,ijabc',lijabc,rijabc,eijabc)
deltaE = 0.5*deltaE.real
logger.info(eom, "eaccsd energy, star energy, delta energy = %16.12f, %16.12f, %16.12f",
eval_, eval_+deltaE, deltaE)
e.append(eval_+deltaE)
return e
[docs]
class EOMEA(EOM):
[docs]
def get_init_guess(self, nroots=1, koopmans=True, diag=None):
size = self.vector_size()
nroots = min(nroots, size)
guess = []
if koopmans:
dtype = getattr(diag, 'dtype', np.double)
for n in range(nroots):
g = np.zeros(size, dtype)
g[n] = 1.0
guess.append(g)
else:
if diag is None:
diag = self.get_diag()
dtype = getattr(diag, 'dtype', np.double)
idx = diag.argsort()[:nroots]
for i in idx:
g = np.zeros(size, dtype)
g[i] = 1.0
guess.append(g)
return guess
kernel = eaccsd
eaccsd = eaccsd
eaccsd_star = eaccsd_star
matvec = eaccsd_matvec
l_matvec = leaccsd_matvec
get_diag = eaccsd_diag
ccsd_star_contract = eaccsd_star_contract
[docs]
def eaccsd_star_contract(self, eaccsd_evals, eaccsd_evecs, leaccsd_evecs, imds=None):
return self.ccsd_star_contract(eaccsd_evals, eaccsd_evecs, leaccsd_evecs, imds=imds)
[docs]
def gen_matvec(self, imds=None, left=False, **kwargs):
if imds is None: imds = self.make_imds()
diag = self.get_diag(imds)
if left:
matvec = lambda xs: [self.l_matvec(x, imds, diag) for x in xs]
else:
matvec = lambda xs: [self.matvec(x, imds, diag) for x in xs]
return matvec, diag
amplitudes_to_vector = staticmethod(amplitudes_to_vector_ea)
vector_to_amplitudes = module_method(vector_to_amplitudes_ea,
absences=['nmo', 'nocc'])
spatial2spin = staticmethod(spatial2spin_ea)
[docs]
def vector_size(self):
nocc = self.nocc
nvir = self.nmo - nocc
return nvir + nocc*nvir*nvir
[docs]
def make_imds(self, eris=None):
imds = _IMDS(self._cc, eris=eris)
imds.make_ea(self.partition)
return imds
@property
def eea(self):
return self.e
[docs]
class EOMEA_Ta(EOMEA):
'''Class for EOM EACCSD(T)*(a) method by Matthews and Stanton.'''
[docs]
def make_imds(self, eris=None):
imds = _IMDS(self._cc, eris=eris)
imds.make_t3p2_ea(self._cc, self.partition)
return imds
########################################
# EOM-EE-CCSD
########################################
#TODO: double spin-flip EOM-EE
[docs]
def eeccsd(eom, nroots=1, koopmans=False, guess=None, eris=None, imds=None):
'''Calculate N-electron neutral excitations via EOM-EE-CCSD.
Kwargs:
nroots : int
Number of roots (eigenvalues) requested
koopmans : bool
Calculate Koopmans'-like (1p1h) excitations only, targeting via
overlap.
guess : list of ndarray
List of guess vectors to use for targeting via overlap.
'''
if eris is None: eris = eom._cc.ao2mo()
if imds is None: imds = eom.make_imds(eris)
spinvec_size = eom.vector_size()
nroots = min(nroots, spinvec_size)
diag_eeS, diag_eeT, diag_sf = eom.get_diag(imds)
guess_eeS = []
guess_eeT = []
guess_sf = []
if guess:
for g in guess:
if g is None: # beta->alpha spin-flip excitation
pass
elif g.size == diag_eeS.size:
guess_eeS.append(g)
elif g.size == diag_eeT.size:
guess_eeT.append(g)
else:
guess_sf.append(g)
nroots_eeS = len(guess_eeS)
nroots_eeT = len(guess_eeT)
nroots_sf = len(guess_sf)
if len(guess) != nroots:
logger.warn(eom, 'Number of states in initial guess %d does not '
'equal to nroots %d.', len(guess), nroots)
else:
deeS = np.sort(diag_eeS)[:nroots]
deeT = np.sort(diag_eeT)[:nroots]
dsf = np.sort(diag_sf)[:nroots]
dmax = np.sort(np.hstack([deeS,deeT,dsf,dsf]))[nroots-1]
nroots_eeS = np.count_nonzero(deeS <= dmax)
nroots_eeT = np.count_nonzero(deeT <= dmax)
nroots_sf = np.count_nonzero(dsf <= dmax)
guess_eeS = guess_eeT = guess_sf = None
def eomee_sub(cls, nroots, guess, diag):
ee_sub = cls(eom._cc)
ee_sub.__dict__.update(eom.__dict__)
e, v = ee_sub.kernel(nroots, koopmans, guess, eris, imds, diag=diag)
if nroots == 1:
e, v = [e], [v]
ee_sub.converged = [ee_sub.converged]
return list(ee_sub.converged), list(e), list(v)
e0 = e1 = e2 = []
v0 = v1 = v2 = []
conv0 = conv1 = conv2 = []
if nroots_eeS > 0:
conv0, e0, v0 = eomee_sub(EOMEESinglet, nroots_eeS, guess_eeS, diag_eeS)
if nroots_eeT > 0:
conv2, e2, v2 = eomee_sub(EOMEETriplet, nroots_eeT, guess_eeT, diag_eeT)
if nroots_sf > 0:
conv1, e1, v1 = eomee_sub(EOMEESpinFlip, nroots_sf, guess_sf, diag_sf)
# The associated solutions of beta->alpha excitations
e1 = e1 + e1
conv1 = conv1 + conv1
v1 = v1 + [None] * len(v1)
# beta->alpha spin-flip excitations, the coefficients are (-r1, (-r2[0], r2[1]))
# as below. The EOMEESpinFlip class only handles alpha->beta excitations.
# Setting beta->alpha to None to bypass the vectors in initial guess
#for i in range(nroots_sf):
# r1, r2 = vector_to_amplitudes_eomsf(v1[i], eom.nmo, eom.nocc)
# v1.append(amplitudes_to_vector_eomsf(-r1, (-r2[0], r2[1])))
e = np.hstack([e0,e2,e1])
idx = e.argsort()
e = e[idx]
conv = conv0 + conv2 + conv1
conv = [conv[x] for x in idx]
v = v0 + v2 + v1
v = [v[x] for x in idx]
if nroots == 1:
conv = conv[0]
e = e[0]
v = v[0]
eom.converged = conv
eom.e = e
eom.v = v
return eom.e, eom.v
[docs]
def eomee_ccsd_singlet(eom, nroots=1, koopmans=False, guess=None,
eris=None, imds=None, diag=None):
'''EOM-EE-CCSD singlet
'''
eom.converged, eom.e, eom.v \
= kernel(eom, nroots, koopmans, guess, eris=eris, imds=imds, diag=diag)
return eom.e, eom.v
[docs]
def eomee_ccsd_triplet(eom, nroots=1, koopmans=False, guess=None,
eris=None, imds=None, diag=None):
'''EOM-EE-CCSD triplet
'''
return eomee_ccsd_singlet(eom, nroots, koopmans, guess, eris, imds, diag)
[docs]
def eomsf_ccsd(eom, nroots=1, koopmans=False, guess=None,
eris=None, imds=None, diag=None):
'''Spin flip EOM-EE-CCSD
'''
return eomee_ccsd_singlet(eom, nroots, koopmans, guess, eris, imds, diag)
vector_to_amplitudes_ee = vector_to_amplitudes_singlet = ccsd.vector_to_amplitudes
amplitudes_to_vector_ee = amplitudes_to_vector_singlet = ccsd.amplitudes_to_vector
[docs]
def spatial2spin_singlet(rx, orbspin=None):
from pyscf.cc import eom_uccsd
if rx.ndim == 2:
r1a = rx * .5**.5
r1 = (r1a, r1a)
return eom_uccsd.spatial2spin_eomee(r1, orbspin)
else:
r2ab = rx * .5**.5
r2aa = r2ab - r2ab.transpose(1,0,2,3)
r2 = (r2aa, r2ab, r2aa)
return eom_uccsd.spatial2spin_eomee(r2, orbspin)
[docs]
def amplitudes_to_vector_eomsf(t1, t2, out=None):
nocc, nvir = t1.shape
t2baaa, t2aaba = t2
otril = np.tril_indices(nocc, k=-1)
vtril = np.tril_indices(nvir, k=-1)
baaa = np.take(t2baaa.reshape(nocc*nocc,nvir*nvir),
vtril[0]*nvir+vtril[1], axis=1)
vector = np.hstack((t1.ravel(), baaa.ravel(), t2aaba[otril].ravel()))
return vector
[docs]
def vector_to_amplitudes_eomsf(vector, nmo, nocc):
nvir = nmo - nocc
t1 = vector[:nocc*nvir].reshape(nocc,nvir).copy()
pvec = vector[t1.size:]
nbaaa = nocc*nocc*nvir*(nvir-1)//2
naaba = nocc*(nocc-1)//2*nvir*nvir
t2baaa = np.zeros((nocc*nocc,nvir*nvir), dtype=vector.dtype)
t2aaba = np.zeros((nocc*nocc,nvir*nvir), dtype=vector.dtype)
otril = np.tril_indices(nocc, k=-1)
vtril = np.tril_indices(nvir, k=-1)
v = pvec[:nbaaa].reshape(nocc*nocc,nvir*(nvir-1)//2)
t2baaa[:,vtril[0]*nvir+vtril[1]] = v
t2baaa[:,vtril[1]*nvir+vtril[0]] = -v
v = pvec[nbaaa:nbaaa+naaba].reshape(-1,nvir*nvir)
t2aaba[otril[0]*nocc+otril[1]] = v
t2aaba[otril[1]*nocc+otril[0]] = -v
t2baaa = t2baaa.reshape(nocc,nocc,nvir,nvir)
t2aaba = t2aaba.reshape(nocc,nocc,nvir,nvir)
return t1, (t2baaa, t2aaba)
[docs]
def spatial2spin_eomsf(rx, orbspin=None):
from pyscf.cc import eom_uccsd
if isinstance(rx, np.ndarray) and rx.ndim == 2:
r1 = (rx, np.zeros_like(rx))
return eom_uccsd.spatial2spin_eomsf(r1, orbspin)
else:
rx0 = np.zeros_like(rx[0])
r2 = (rx0, rx[1], rx[0], rx0)
return eom_uccsd.spatial2spin_eomsf(r2, orbspin)
[docs]
def amplitudes_to_vector_triplet(t1, t2, out=None):
t2aa, t2ab = t2
dtype = np.result_type(t1, t2aa, t2ab)
nocc, nvir = t1.shape
nov = nocc * nvir
size1 = nov + nocc*(nocc-1)//2*nvir*(nvir-1)//2
size = size1 + nov*(nov+1)//2
vector = np.ndarray(size, dtype, buffer=out)
ccsd.amplitudes_to_vector_s4(t1, t2[0], out=vector)
t2ab = t2[1].transpose(0,2,1,3).reshape(nov,nov)
lib.pack_tril(t2ab, out=vector[size1:])
return vector
[docs]
def vector_to_amplitudes_triplet(vector, nmo, nocc):
nvir = nmo - nocc
nov = nocc * nvir
size1 = nov + nocc*(nocc-1)//2*nvir*(nvir-1)//2
size = size1 + nov*(nov+1)//2
t1, t2aa = ccsd.vector_to_amplitudes_s4(vector[:size1], nmo, nocc)
t2ab = lib.unpack_tril(vector[size1:size], filltriu=2)
t2ab = t2ab.reshape(nocc,nvir,nocc,nvir).transpose(0,2,1,3).copy()
return t1, (t2aa, t2ab)
[docs]
def spatial2spin_triplet(rx, orbspin=None):
from pyscf.cc import eom_uccsd
if isinstance(rx, np.ndarray) and rx.ndim == 2:
r1a = rx * .5**.5
r1 = (r1a, r1a)
return eom_uccsd.spatial2spin_eomee(r1, orbspin)
else:
r2ab = rx * .5**.5
r2aa, r2ab = rx
r2aa = r2ab - r2ab.transpose(1,0,2,3)
r2 = (r2aa, r2ab, -r2aa)
return eom_uccsd.spatial2spin_eomee(r2, orbspin)
[docs]
def eeccsd_matvec_singlet(eom, vector, imds=None):
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
nvir = nmo - nocc
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
t1, t2, eris = imds.t1, imds.t2, imds.eris
nocc, nvir = t1.shape
Hr1 = lib.einsum('ae,ie->ia', imds.Fvv, r1)
Hr1 -= lib.einsum('mi,ma->ia', imds.Foo, r1)
Hr1 += np.einsum('me,imae->ia',imds.Fov, r2) * 2
Hr1 -= np.einsum('me,imea->ia',imds.Fov, r2)
#:eris_vvvv = ao2mo.restore(1,np.asarray(eris.vvvv), t1.shape[1])
#:Hr2 += lib.einsum('ijef,aebf->ijab', tau2, eris_vvvv) * .5
tau2 = _make_tau(r2, r1, t1, fac=2)
Hr2 = eom._cc._add_vvvv(None, tau2, eris, with_ovvv=False, t2sym='jiba')
woOoO = np.asarray(imds.woOoO)
Hr2 += lib.einsum('mnij,mnab->ijab', woOoO, r2)
Hr2 *= .5
woOoO = None
Hr2 += lib.einsum('be,ijae->ijab', imds.Fvv , r2)
Hr2 -= lib.einsum('mj,imab->ijab', imds.Foo , r2)
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now - Hr2.size*8e-6)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*3))))
for p0,p1 in lib.prange(0, nocc, blksize):
ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1]
theta = r2[p0:p1] * 2 - r2[p0:p1].transpose(0,1,3,2)
Hr1 += lib.einsum('mfae,mife->ia', ovvv, theta)
theta = None
tmp = lib.einsum('meaf,ijef->maij', ovvv, tau2)
Hr2 -= lib.einsum('ma,mbij->ijab', t1[p0:p1], tmp)
tmp = lib.einsum('meaf,me->af', ovvv, r1[p0:p1]) * 2
tmp -= lib.einsum('mfae,me->af', ovvv, r1[p0:p1])
Hr2 += lib.einsum('af,ijfb->ijab', tmp, t2)
ovvv = tmp = None
tau2 = None
Hr2 -= lib.einsum('mbij,ma->ijab', imds.woVoO, r1)
blksize = min(nvir, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nocc*nvir**2*2))))
for p0, p1 in lib.prange(0, nvir, nocc):
Hr2 += lib.einsum('ejab,ie->ijab', np.asarray(imds.wvOvV[p0:p1]), r1[:,p0:p1])
woVVo = np.asarray(imds.woVVo)
tmp = lib.einsum('mbej,imea->jiab', woVVo, r2)
Hr2 += tmp
tmp *= .5
Hr2 += tmp.transpose(0,1,3,2)
tmp = None
woVvO = woVVo * .5
woVVo = None
woVvO += np.asarray(imds.woVvO)
theta = r2*2 - r2.transpose(0,1,3,2)
Hr1 += np.einsum('maei,me->ia', woVvO, r1) * 2
Hr2 += lib.einsum('mbej,imae->ijab', woVvO, theta)
woVvO = None
woOoV = np.asarray(imds.woOoV)
Hr1-= lib.einsum('mnie,mnae->ia', woOoV, theta)
tmp = lib.einsum('nmie,me->ni', woOoV, r1) * 2
tmp-= lib.einsum('mnie,me->ni', woOoV, r1)
Hr2 -= lib.einsum('ni,njab->ijab', tmp, t2)
tmp = woOoV = None
eris_ovov = np.asarray(eris.ovov)
tmp = np.einsum('mfne,mf->en', eris_ovov, r1) * 2
tmp -= np.einsum('menf,mf->en', eris_ovov, r1)
tmp = np.einsum('en,nb->eb', tmp, t1)
tmp += lib.einsum('menf,mnbf->eb', eris_ovov, theta)
Hr2 -= lib.einsum('eb,ijea->jiab', tmp, t2)
tmp = None
tmp = lib.einsum('nemf,imef->ni', eris_ovov, theta)
Hr1 -= lib.einsum('na,ni->ia', t1, tmp)
Hr2 -= lib.einsum('mj,miab->ijba', tmp, t2)
tmp = theta = None
tau2 = _make_tau(r2, r1, t1, fac=2)
tmp = lib.einsum('menf,ijef->mnij', eris_ovov, tau2)
tau2 = None
tau = _make_tau(t2, t1, t1)
tau *= .5
Hr2 += lib.einsum('mnij,mnab->ijab', tmp, tau)
tau = tmp = eris_ovov = None
Hr2 = Hr2 + Hr2.transpose(1,0,3,2)
vector = eom.amplitudes_to_vector(Hr1, Hr2)
return vector
[docs]
def eeccsd_matvec_triplet(eom, vector, imds=None):
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
nvir = nmo - nocc
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
r2aa, r2ab = r2
t1, t2, eris = imds.t1, imds.t2, imds.eris
nocc, nvir = t1.shape
Hr1 = lib.einsum('ae,ie->ia', imds.Fvv, r1)
Hr1 -= lib.einsum('mi,ma->ia', imds.Foo, r1)
Hr1 += np.einsum('me,imae->ia',imds.Fov, r2aa)
Hr1 += np.einsum('ME,iMaE->ia',imds.Fov, r2ab)
tau2ab = np.einsum('ia,jb->ijab', r1, t1)
tau2ab-= np.einsum('ia,jb->ijab', t1, r1)
tau2ab+= r2ab
tau2aa = np.einsum('ia,jb->ijab', r1, t1)
tau2aa-= np.einsum('ia,jb->jiab', r1, t1)
tau2aa = tau2aa - tau2aa.transpose(0,1,3,2)
tau2aa+= r2aa
#:eris_vvvv = ao2mo.restore(1,np.asarray(eris.vvvv), t1.shape[1])
#:Hr2aa += lib.einsum('ijef,aebf->ijab', tau2aa, eris_vvvv) * .25
#:Hr2ab += lib.einsum('ijef,aebf->ijab', tau2ab, eris_vvvv) * .5
Hr2aa = eom._cc._add_vvvv(None, tau2aa, eris, with_ovvv=False, t2sym='jiba')
Hr2ab = eom._cc._add_vvvv(None, tau2ab, eris, with_ovvv=False, t2sym='-jiba')
woOoO = np.asarray(imds.woOoO)
Hr2aa += lib.einsum('mnij,mnab->ijab', woOoO, r2aa)
Hr2ab += lib.einsum('mNiJ,mNaB->iJaB', woOoO, r2ab)
Hr2aa *= .25
Hr2ab *= .5
woOoO = None
Hr2aa += lib.einsum('be,ijae->ijab', imds.Fvv*.5, r2aa)
Hr2aa -= lib.einsum('mj,imab->ijab', imds.Foo*.5, r2aa)
Hr2ab += lib.einsum('BE,iJaE->iJaB', imds.Fvv, r2ab)
Hr2ab -= lib.einsum('MJ,iMaB->iJaB', imds.Foo, r2ab)
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now - Hr2aa.size*8e-6)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*3))))
tmp1 = np.zeros((nvir,nvir), dtype=r1.dtype)
for p0,p1 in lib.prange(0, nocc, blksize):
ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1]
theta = r2aa[:,p0:p1] + r2ab[:,p0:p1]
Hr1 += lib.einsum('mfae,imef->ia', ovvv, theta)
theta = None
tmpaa = lib.einsum('meaf,ijef->maij', ovvv, tau2aa)
tmpab = lib.einsum('meAF,iJeF->mAiJ', ovvv, tau2ab)
tmp1 += lib.einsum('mfae,me->af', ovvv, r1[p0:p1])
Hr2aa+= lib.einsum('mb,maij->ijab', t1[p0:p1]*.5, tmpaa)
Hr2ab-= lib.einsum('mb,mAiJ->iJbA', t1[p0:p1], tmpab)
ovvv = tmpaa = tmpab = None
tau2aa = tau2ab = None
woVVo = np.asarray(imds.woVVo)
Hr1 += np.einsum('maei,me->ia', woVVo, r1)
Hr2aa += lib.einsum('mbej,imae->ijba', woVVo, r2ab)
Hr2ab += lib.einsum('MBEJ,iMEa->iJaB', woVVo, r2aa)
Hr2ab += lib.einsum('MbeJ,iMeA->iJbA', woVVo, r2ab)
woVvO = np.asarray(imds.woVvO)
wovvo = woVVo + woVvO
theta = r2aa + r2ab
tmp = lib.einsum('mbej,imae->ijab', wovvo, theta)
woVVo = woVvO = wovvo = None
woOoV = np.asarray(imds.woOoV)
Hr1 -= lib.einsum('mnie,mnae->ia', woOoV, theta)
tmpa = lib.einsum('mnie,me->ni', woOoV, r1)
tmp += lib.einsum('ni,njab->ijab', tmpa, t2)
tmp -= lib.einsum('af,ijfb->ijab', tmp1, t2)
tmp -= lib.einsum('mbij,ma->ijab', imds.woVoO, r1)
blksize = min(nvir, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nocc*nvir**2*2))))
for p0,p1 in lib.prange(0, nvir, blksize):
tmp += lib.einsum('ejab,ie->ijab', np.asarray(imds.wvOvV[p0:p1]), r1[:,p0:p1])
Hr2aa += tmp
Hr2ab += tmp
tmp = woOoV = None
eris_ovov = np.asarray(eris.ovov)
tmpa = -lib.einsum('menf,imfe->ni', eris_ovov, theta)
Hr1 += lib.einsum('na,ni->ia', t1, tmpa)
tmp = lib.einsum('mj,imab->ijab', tmpa, t2)
tmp1 = np.einsum('menf,mf->en', eris_ovov, r1)
tmpa = lib.einsum('en,nb->eb', tmp1, t1)
tmpa-= lib.einsum('menf,mnbf->eb', eris_ovov, theta)
tmp += lib.einsum('eb,ijae->ijab', tmpa, t2)
Hr2aa += tmp
Hr2ab -= tmp
tmp = theta = tmp1 = tmpa = None
tau2aa = np.einsum('ia,jb->ijab', r1, t1)
tau2aa-= np.einsum('ia,jb->jiab', r1, t1)
tau2aa = tau2aa - tau2aa.transpose(0,1,3,2)
tau2aa+= r2aa
tmpaa = lib.einsum('menf,ijef->mnij', eris_ovov, tau2aa)
tau2aa = None
tmpaa *= .25
tau = _make_tau(t2, t1, t1)
Hr2aa += lib.einsum('mnij,mnab->ijab', tmpaa, tau)
tmpaa = tau = None
tau2ab = np.einsum('ia,jb->ijab', r1, t1)
tau2ab-= np.einsum('ia,jb->ijab', t1, r1)
tau2ab+= r2ab
tmpab = lib.einsum('meNF,iJeF->mNiJ', eris_ovov, tau2ab)
tau2ab = None
tmpab *= .5
tau = _make_tau(t2, t1, t1)
Hr2ab += lib.einsum('mNiJ,mNaB->iJaB', tmpab, tau)
tmpab = tau = None
eris_ovov = None
Hr2aa = Hr2aa - Hr2aa.transpose(0,1,3,2)
Hr2aa = Hr2aa - Hr2aa.transpose(1,0,2,3)
Hr2ab = Hr2ab - Hr2ab.transpose(1,0,3,2)
vector = eom.amplitudes_to_vector(Hr1, (Hr2aa,Hr2ab))
return vector
[docs]
def eeccsd_matvec_sf(eom, vector, imds=None):
'''Spin flip EOM-CCSD'''
if imds is None: imds = eom.make_imds()
nocc = eom.nocc
nmo = eom.nmo
nvir = nmo - nocc
t1, t2, eris = imds.t1, imds.t2, imds.eris
r1, r2 = eom.vector_to_amplitudes(vector, nmo, nocc)
r2baaa, r2aaba = r2
nocc, nvir = t1.shape
Hr1 = np.einsum('ae,ie->ia', imds.Fvv, r1)
Hr1 -= np.einsum('mi,ma->ia', imds.Foo, r1)
Hr1 += np.einsum('me,imae->ia', imds.Fov, r2baaa)
Hr1 += np.einsum('me,imae->ia', imds.Fov, r2aaba)
tau2baaa = np.einsum('ia,jb->ijab', r1, t1)
tau2baaa += r2baaa * .5
tau2baaa = tau2baaa - tau2baaa.transpose(0,1,3,2)
tau2aaba = np.einsum('ia,jb->ijab', r1, t1)
tau2aaba += r2aaba * .5
tau2aaba = tau2aaba - tau2aaba.transpose(1,0,2,3)
#:eris_vvvv = ao2mo.restore(1,np.asarray(eris.vvvv), t1.shape[1])
#:Hr2baaa += .5*lib.einsum('ijef,aebf->ijab', tau2baaa, eris_vvvv)
#:Hr2aaba += .5*lib.einsum('ijef,aebf->ijab', tau2aaba, eris_vvvv)
Hr2aaba = eom._cc._add_vvvv(None, tau2aaba, eris, with_ovvv=False, t2sym='-jiab')
Hr2baaa = eom._cc._add_vvvv(None, tau2baaa, eris, with_ovvv=False, t2sym=False)
woOoO = np.asarray(imds.woOoO)
Hr2baaa += lib.einsum('mnij,mnab->ijab', woOoO, r2baaa)
Hr2aaba += lib.einsum('mnij,mnab->ijab', woOoO, r2aaba)
Hr2aaba *= .5
Hr2baaa *= .5
woOoO = None
Hr2baaa -= lib.einsum('mj,imab->ijab', imds.Foo*.5, r2baaa)
Hr2aaba -= lib.einsum('mj,imab->ijab', imds.Foo*.5, r2aaba)
Hr2baaa -= lib.einsum('mj,miab->jiab', imds.Foo*.5, r2baaa)
Hr2aaba -= lib.einsum('mj,miab->jiab', imds.Foo*.5, r2aaba)
Hr2baaa += lib.einsum('be,ijae->ijab', imds.Fvv*.5, r2baaa)
Hr2aaba += lib.einsum('be,ijae->ijab', imds.Fvv*.5, r2aaba)
Hr2baaa += lib.einsum('be,ijea->ijba', imds.Fvv*.5, r2baaa)
Hr2aaba += lib.einsum('be,ijea->ijba', imds.Fvv*.5, r2aaba)
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now - Hr2aaba.size*8e-6)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*3))))
for p0,p1 in lib.prange(0, nocc, blksize):
ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1]
theta = r2baaa[:,p0:p1] + r2aaba[:,p0:p1]
Hr1 += lib.einsum('mfae,imef->ia', ovvv, theta)
theta = None
tmp1aaba = lib.einsum('meaf,ijef->maij', ovvv, tau2baaa)
Hr2baaa -= lib.einsum('mb,maij->ijba', t1[p0:p1]*.5, tmp1aaba)
tmp1aaba = tmp1baaa = tmp1abaa = tmp2aaba = None
tmp2aaba = lib.einsum('meaf,ijfe->maij', ovvv, tau2baaa)
Hr2baaa -= lib.einsum('mb,maij->ijab', t1[p0:p1]*.5, tmp2aaba)
tmp1aaba = tmp1baaa = tmp1abaa = tmp2aaba = None
tmp1baaa = lib.einsum('meaf,ijef->maij', ovvv, tau2aaba)
Hr2aaba -= lib.einsum('mb,maij->ijba', t1[p0:p1]*.5, tmp1baaa)
tmp1aaba = tmp1baaa = tmp1abaa = tmp2aaba = None
tmp1abaa = lib.einsum('meaf,ijfe->maij', ovvv, tau2aaba)
Hr2aaba -= lib.einsum('mb,maij->ijab', t1[p0:p1]*.5, tmp1abaa)
tmp1aaba = tmp1baaa = tmp1abaa = tmp2aaba = None
tmp = lib.einsum('mfae,me->af', ovvv, r1[p0:p1])
tmp = lib.einsum('af,jibf->ijab', tmp, t2)
Hr2baaa -= tmp
Hr2aaba -= tmp
tmp = ovvv = None
tau2aaba = tau2baaa = None
tmp = lib.einsum('mbij,ma->ijab', imds.woVoO, r1)
Hr2baaa -= tmp
Hr2aaba -= tmp
tmp = None
woOoV = np.asarray(imds.woOoV)
Hr1 -= lib.einsum('mnie,mnae->ia', woOoV, r2aaba)
Hr1 -= lib.einsum('mnie,mnae->ia', woOoV, r2baaa)
tmp = lib.einsum('mnie,me->ni', woOoV, r1)
tmp = lib.einsum('ni,njab->ijab', tmp, t2)
Hr2baaa += tmp
Hr2aaba += tmp
tmp = woOoV = None
for p0,p1 in lib.prange(0, nvir, nocc):
tmp = lib.einsum('ejab,ie->ijab', np.asarray(imds.wvOvV[p0:p1]), r1[:,p0:p1])
Hr2baaa += tmp
Hr2aaba += tmp
tmp = None
woVVo = np.asarray(imds.woVVo)
Hr1 += np.einsum('maei,me->ia', woVVo, r1)
Hr2baaa += lib.einsum('mbej,miea->jiba', woVVo, r2baaa)
Hr2aaba += lib.einsum('mbej,miea->jiba', woVVo, r2aaba)
woVVo = None
woVvO = np.asarray(imds.woVvO)
Hr2baaa += lib.einsum('mbej,imae->ijab', woVvO, r2aaba)
Hr2aaba += lib.einsum('mbej,imae->ijab', woVvO, r2baaa)
woVvO = woVvO + np.asarray(imds.woVVo)
Hr2baaa += lib.einsum('mbej,imae->ijab', woVvO, r2baaa)
Hr2aaba += lib.einsum('mbej,imae->ijab', woVvO, r2aaba)
woVvO = None
eris_ovov = np.asarray(eris.ovov)
theta = r2aaba + r2baaa
tmp = lib.einsum('nfme,imfe->ni', eris_ovov, theta)
Hr1 -= np.einsum('na,ni->ia', t1, tmp)
Hr2baaa -= lib.einsum('mj,imba->jiab', tmp, t2)
Hr2aaba -= lib.einsum('mj,imba->jiab', tmp, t2)
tmp = np.einsum('menf,mf->en', eris_ovov, r1)
tmp = np.einsum('en,nb->eb', tmp, t1)
tmp-= lib.einsum('menf,mnbf->eb', eris_ovov, theta)
Hr2baaa += lib.einsum('ea,ijbe->jiab', tmp, t2)
Hr2aaba += lib.einsum('ea,ijbe->jiab', tmp, t2)
theta = tmp = None
tau2baaa = np.einsum('ia,jb->ijab', r1, t1)
tau2baaa += r2baaa * .5
tau2baaa = tau2baaa - tau2baaa.transpose(0,1,3,2)
tau = _make_tau(t2, t1, t1)
tmp1aaba = lib.einsum('menf,ijef->mnij', eris_ovov, tau2baaa)
tau2baaa = None
Hr2baaa += .5*lib.einsum('mnij,mnab->ijab', tmp1aaba, tau)
tau = tmp1aaba = None
tau2aaba = np.einsum('ia,jb->ijab', r1, t1)
tau2aaba += r2aaba * .5
tau2aaba = tau2aaba - tau2aaba.transpose(1,0,2,3)
tau = _make_tau(t2, t1, t1)
tmp1baaa = lib.einsum('menf,ijef->mnij', eris_ovov, tau2aaba)
tau2aaba = None
Hr2aaba += .5*lib.einsum('mnij,mnab->ijab', tmp1baaa, tau)
tau = tmp1baaa = None
eris_ovov = None
Hr2baaa = Hr2baaa - Hr2baaa.transpose(0,1,3,2)
Hr2aaba = Hr2aaba - Hr2aaba.transpose(1,0,2,3)
vector = eom.amplitudes_to_vector(Hr1, (Hr2baaa,Hr2aaba))
return vector
[docs]
def eeccsd_diag(eom, imds=None):
if imds is None: imds = eom.make_imds()
eris = imds.eris
t1, t2 = imds.t1, imds.t2
dtype = np.result_type(t1, t2)
tau = _make_tau(t2, t1, t1)
nocc, nvir = t1.shape
Fo = imds.Foo.diagonal()
Fv = imds.Fvv.diagonal()
Wovab = np.einsum('iaai->ia', imds.woVVo)
Wovaa = Wovab + np.einsum('iaai->ia', imds.woVvO)
eia = lib.direct_sum('-i+a->ia', Fo, Fv)
Hr1aa = eia + Wovaa
Hr1ab = eia + Wovab
eris_ovov = np.asarray(eris.ovov)
Wvvab = np.einsum('mnab,manb->ab', tau, eris_ovov)
Wvvaa = .5*Wvvab - .5*np.einsum('mnba,manb->ab', tau, eris_ovov)
ijb = np.einsum('iejb,ijeb->ijb', eris_ovov, t2)
Hr2ab = lib.direct_sum('iJB+a->iJaB',-ijb, Fv)
jab = np.einsum('kajb,kjab->jab', eris_ovov, t2)
Hr2ab+= lib.direct_sum('-i-jab->ijab', Fo, jab)
jib = np.einsum('iejb,ijbe->jib', eris_ovov, t2)
jib = jib + jib.transpose(1,0,2)
jib-= ijb + ijb.transpose(1,0,2)
jba = np.einsum('kajb,jkab->jba', eris_ovov, t2)
jba = jba + jba.transpose(0,2,1)
jba-= jab + jab.transpose(0,2,1)
Hr2aa = lib.direct_sum('jib+a->jiba', jib, Fv)
Hr2aa+= lib.direct_sum('-i+jba->ijba', Fo, jba)
eris_ovov = None
Hr2baaa = lib.direct_sum('ijb+a->ijba',-ijb, Fv)
Hr2baaa += Wovaa.reshape(1,nocc,1,nvir)
Hr2baaa += Wovab.reshape(nocc,1,1,nvir)
Hr2baaa = Hr2baaa + Hr2baaa.transpose(0,1,3,2)
Hr2baaa+= lib.direct_sum('-i+jab->ijab', Fo, jba)
Hr2baaa-= Fo.reshape(1,-1,1,1)
Hr2aaba = lib.direct_sum('-i-jab->ijab', Fo, jab)
Hr2aaba += Wovaa.reshape(1,nocc,1,nvir)
Hr2aaba += Wovab.reshape(1,nocc,nvir,1)
Hr2aaba = Hr2aaba + Hr2aaba.transpose(1,0,2,3)
Hr2aaba+= lib.direct_sum('ijb+a->ijab', jib, Fv)
Hr2aaba+= Fv.reshape(1,1,1,-1)
Hr2ab += Wovaa.reshape(1,nocc,1,nvir)
Hr2ab += Wovab.reshape(nocc,1,1,nvir)
Hr2ab = Hr2ab + Hr2ab.transpose(1,0,3,2)
Hr2aa += Wovaa.reshape(1,nocc,1,nvir) * 2
Hr2aa = Hr2aa + Hr2aa.transpose(0,1,3,2)
Hr2aa = Hr2aa + Hr2aa.transpose(1,0,2,3)
Hr2aa *= .5
Wooab = np.einsum('ijij->ij', imds.woOoO)
Wooaa = Wooab - np.einsum('ijji->ij', imds.woOoO)
Hr2aa += Wooaa.reshape(nocc,nocc,1,1)
Hr2ab += Wooab.reshape(nocc,nocc,1,1)
Hr2baaa += Wooab.reshape(nocc,nocc,1,1)
Hr2aaba += Wooaa.reshape(nocc,nocc,1,1)
#:eris_ovvv = lib.unpack_tril(np.asarray(eris.ovvv).reshape(nocc*nvir,nvir**2)).reshape(nocc,nvir,nvir,nvir)
#:tmp = np.einsum('mb,mbaa->ab', t1, eris_ovvv)
#:Wvvaa += np.einsum('mb,maab->ab', t1, eris_ovvv)
mem_now = lib.current_memory()[0]
max_memory = max(0, eom.max_memory - mem_now)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*3))))
tmp = np.zeros((nvir,nvir), dtype=dtype)
for p0,p1 in lib.prange(0, nocc, blksize):
ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1]
tmp += np.einsum('mb,mbaa->ab', t1[p0:p1], ovvv)
Wvvaa += np.einsum('mb,maab->ab', t1[p0:p1], ovvv)
ovvv = None
Wvvaa -= tmp
Wvvab -= tmp
Wvvab -= tmp.T
Wvvaa = Wvvaa + Wvvaa.T
if eris.vvvv is None: # AO-direct CCSD, vvvv is not generated.
pass
elif eris.vvvv.ndim == 4:
eris_vvvv = ao2mo.restore(1,np.asarray(eris.vvvv), t1.shape[1])
tmp = np.einsum('aabb->ab', eris_vvvv)
Wvvaa += tmp
Wvvaa -= np.einsum('abba->ab', eris_vvvv)
Wvvab += tmp
else:
for i in range(nvir):
i0 = i*(i+1)//2
vvv = lib.unpack_tril(np.asarray(eris.vvvv[i0:i0+i+1]))
tmp = np.einsum('bb->b', vvv[i])
Wvvaa[i] += tmp
Wvvab[i] += tmp
tmp = np.einsum('bb->b', vvv[:,:i+1,i])
Wvvaa[i,:i+1] -= tmp
Wvvaa[:i ,i] -= tmp[:i]
vvv = None
Hr2aa += Wvvaa.reshape(1,1,nvir,nvir)
Hr2ab += Wvvab.reshape(1,1,nvir,nvir)
Hr2baaa += Wvvaa.reshape(1,1,nvir,nvir)
Hr2aaba += Wvvab.reshape(1,1,nvir,nvir)
vec_eeS = amplitudes_to_vector_singlet(Hr1aa, Hr2ab)
vec_eeT = amplitudes_to_vector_triplet(Hr1aa, (Hr2aa,Hr2ab))
vec_sf = amplitudes_to_vector_eomsf(Hr1ab, (Hr2baaa,Hr2aaba))
return vec_eeS, vec_eeT, vec_sf
[docs]
class EOMEE(EOM):
[docs]
def get_init_guess(self, nroots=1, koopmans=True, diag=None):
if diag is None:
diag = self.get_diag()
if koopmans:
nocc = self.nocc
nvir = self.nmo - nocc
idx = diag[:nocc*nvir].argsort()
else:
idx = diag.argsort()
size = self.vector_size()
dtype = getattr(diag, 'dtype', np.double)
nroots = min(nroots, size)
guess = []
for i in idx[:nroots]:
g = np.zeros(size, dtype)
g[i] = 1.0
guess.append(g)
return guess
kernel = eeccsd
eeccsd = eeccsd
get_diag = eeccsd_diag
[docs]
def vector_size(self):
'''size of the vector based on spin-orbital basis'''
nocc = self.nocc
nvir = self.nmo - nocc
return nocc*nvir + nocc*nocc*nvir*nvir
[docs]
def make_imds(self, eris=None):
imds = _IMDS(self._cc, eris=eris)
imds.make_ee()
return imds
@property
def eee(self):
return self.e
[docs]
class EOMEESinglet(EOMEE):
kernel = eomee_ccsd_singlet
eomee_ccsd_singlet = eomee_ccsd_singlet
matvec = eeccsd_matvec_singlet
[docs]
def get_diag(self, imds=None):
return eeccsd_diag(self, imds=None)[0]
[docs]
def gen_matvec(self, imds=None, diag=None, **kwargs):
if imds is None: imds = self.make_imds()
if diag is None: diag = self.get_diag(imds)
matvec = lambda xs: [self.matvec(x, imds) for x in xs]
return matvec, diag
amplitudes_to_vector = staticmethod(amplitudes_to_vector_singlet)
vector_to_amplitudes = module_method(vector_to_amplitudes_singlet,
absences=['nmo', 'nocc'])
spatial2spin = staticmethod(spatial2spin_singlet)
[docs]
def vector_size(self):
nocc = self.nocc
nvir = self.nmo - nocc
nov = nocc * nvir
return nov + nov*(nov+1)//2
[docs]
class EOMEETriplet(EOMEE):
kernel = eomee_ccsd_triplet
eomee_ccsd_triplet = eomee_ccsd_triplet
matvec = eeccsd_matvec_triplet
[docs]
def get_diag(self, imds=None):
return eeccsd_diag(self, imds=None)[1]
[docs]
def gen_matvec(self, imds=None, diag=None, **kwargs):
if imds is None: imds = self.make_imds()
if diag is None: diag = self.get_diag(imds)
matvec = lambda xs: [self.matvec(x, imds) for x in xs]
return matvec, diag
amplitudes_to_vector = staticmethod(amplitudes_to_vector_triplet)
vector_to_amplitudes = module_method(vector_to_amplitudes_triplet,
absences=['nmo', 'nocc'])
spatial2spin = staticmethod(spatial2spin_triplet)
[docs]
def vector_size(self):
nocc = self.nocc
nvir = self.nmo - nocc
nov = nocc * nvir
return nov + nocc*(nocc-1)//2*nvir*(nvir-1)//2 + nov*(nov+1)//2
[docs]
class EOMEESpinFlip(EOMEE):
kernel = eomsf_ccsd
eomsf_ccsd = eomsf_ccsd
matvec = eeccsd_matvec_sf
[docs]
def get_diag(self, imds=None):
return eeccsd_diag(self, imds=None)[2]
[docs]
def gen_matvec(self, imds=None, diag=None, **kwargs):
if imds is None: imds = self.make_imds()
if diag is None: diag = self.get_diag(imds)
matvec = lambda xs: [self.matvec(x, imds) for x in xs]
return matvec, diag
amplitudes_to_vector = staticmethod(amplitudes_to_vector_eomsf)
vector_to_amplitudes = module_method(vector_to_amplitudes_eomsf,
absences=['nmo', 'nocc'])
spatial2spin = staticmethod(spatial2spin_eomsf)
[docs]
def vector_size(self):
nocc = self.nocc
nvir = self.nmo - nocc
nbaaa = nocc*nocc*nvir*(nvir-1)//2
naaba = nocc*(nocc-1)//2*nvir*nvir
return nocc*nvir + nbaaa + naaba
#TODO: Check whether EOM methods works with rccsd.RCCSD when orbitals are complex
ccsd.CCSD.EOMIP = lib.class_as_method(EOMIP)
ccsd.CCSD.EOMIP_Ta = lib.class_as_method(EOMIP_Ta)
ccsd.CCSD.EOMEA = lib.class_as_method(EOMEA)
ccsd.CCSD.EOMEA_Ta = lib.class_as_method(EOMEA_Ta)
ccsd.CCSD.EOMEE = lib.class_as_method(EOMEE)
ccsd.CCSD.EOMEESinglet = lib.class_as_method(EOMEESinglet)
ccsd.CCSD.EOMEETriplet = lib.class_as_method(EOMEETriplet)
ccsd.CCSD.EOMEESpinFlip = lib.class_as_method(EOMEESpinFlip)
class _IMDS:
def __init__(self, cc, eris=None):
self.verbose = cc.verbose
self.stdout = cc.stdout
self.max_memory = cc.max_memory
self.t1 = cc.t1
self.t2 = cc.t2
if eris is None:
eris = cc.ao2mo()
self.eris = eris
self._made_shared_2e = False
def _make_shared_1e(self):
cput0 = (logger.process_clock(), logger.perf_counter())
t1, t2, eris = self.t1, self.t2, self.eris
self.Loo = imd.Loo(t1, t2, eris)
self.Lvv = imd.Lvv(t1, t2, eris)
self.Fov = imd.cc_Fov(t1, t2, eris)
logger.timer_debug1(self, 'EOM-CCSD shared one-electron '
'intermediates', *cput0)
return self
def _make_shared_2e(self):
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.stdout, self.verbose)
t1, t2, eris = self.t1, self.t2, self.eris
# 2 virtuals
self.Wovov = imd.Wovov(t1, t2, eris)
self.Wovvo = imd.Wovvo(t1, t2, eris)
self.Woovv = np.asarray(eris.ovov).transpose(0,2,1,3)
self._made_shared_2e = True
log.timer_debug1('EOM-CCSD shared two-electron intermediates', *cput0)
return self
def make_ip(self, ip_partition=None):
self._make_shared_1e()
if not self._made_shared_2e and ip_partition != 'mp':
self._make_shared_2e()
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.stdout, self.verbose)
t1, t2, eris = self.t1, self.t2, self.eris
# 0 or 1 virtuals
if ip_partition != 'mp':
self.Woooo = imd.Woooo(t1, t2, eris)
self.Wooov = imd.Wooov(t1, t2, eris)
self.Wovoo = imd.Wovoo(t1, t2, eris)
log.timer_debug1('EOM-CCSD IP intermediates', *cput0)
return self
def make_t3p2_ip(self, cc, ip_partition=None):
assert (ip_partition is None)
cput0 = (logger.process_clock(), logger.perf_counter())
t1, t2, eris = cc.t1, cc.t2, self.eris
delta_E_corr, pt1, pt2, Wovoo, Wvvvo = \
imd.get_t3p2_imds_slow(cc, t1, t2, eris)
self.t1 = pt1
self.t2 = pt2
self._made_shared_2e = False # Force update
self.make_ip() # Make after t1/t2 updated
self.Wovoo = self.Wovoo + Wovoo
logger.timer_debug1(self, 'EOM-CCSD(T)a IP intermediates', *cput0)
return self
def make_ea(self, ea_partition=None):
self._make_shared_1e()
if not self._made_shared_2e and ea_partition != 'mp':
self._make_shared_2e()
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.stdout, self.verbose)
t1, t2, eris = self.t1, self.t2, self.eris
# 3 or 4 virtuals
self.Wvovv = imd.Wvovv(t1, t2, eris)
if ea_partition == 'mp':
self.Wvvvo = imd.Wvvvo(t1, t2, eris)
else:
self.Wvvvv = imd.Wvvvv(t1, t2, eris)
self.Wvvvo = imd.Wvvvo(t1, t2, eris, self.Wvvvv)
log.timer_debug1('EOM-CCSD EA intermediates', *cput0)
return self
def make_t3p2_ea(self, cc, ea_partition=None):
assert (ea_partition is None)
cput0 = (logger.process_clock(), logger.perf_counter())
t1, t2, eris = cc.t1, cc.t2, self.eris
delta_E_corr, pt1, pt2, Wovoo, Wvvvo = \
imd.get_t3p2_imds_slow(cc, t1, t2, eris)
self.t1 = pt1
self.t2 = pt2
self._made_shared_2e = False # Force update
self.make_ea() # Make after t1/t2 updated
self.Wvvvo = self.Wvvvo + Wvvvo
logger.timer_debug1(self, 'EOM-CCSD(T)a EA intermediates', *cput0)
return self
def make_ee(self):
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(self.stdout, self.verbose)
t1, t2, eris = self.t1, self.t2, self.eris
dtype = np.result_type(t1, t2)
if np.iscomplexobj(t2):
raise NotImplementedError('Complex integrals are not supported in EOM-EE-CCSD')
nocc, nvir = t1.shape
fswap = lib.H5TmpFile()
self.saved = lib.H5TmpFile()
self.wvOvV = self.saved.create_dataset('wvOvV', (nvir,nocc,nvir,nvir), dtype.char)
self.woVvO = self.saved.create_dataset('woVvO', (nocc,nvir,nvir,nocc), dtype.char)
self.woVVo = self.saved.create_dataset('woVVo', (nocc,nvir,nvir,nocc), dtype.char)
self.woOoV = self.saved.create_dataset('woOoV', (nocc,nocc,nocc,nvir), dtype.char)
foo = eris.fock[:nocc,:nocc]
fov = eris.fock[:nocc,nocc:]
fvv = eris.fock[nocc:,nocc:]
self.Fov = np.zeros((nocc,nvir), dtype=dtype)
self.Foo = np.zeros((nocc,nocc), dtype=dtype)
self.Fvv = np.zeros((nvir,nvir), dtype=dtype)
#:eris_ovvv = lib.unpack_tril(np.asarray(eris.ovvv).reshape(nocc*nvir,nvir**2)).reshape(nocc,nvir,nvir,nvir)
#:self.Fvv = np.einsum('mf,mfae->ae', t1, eris_ovvv) * 2
#:self.Fvv -= np.einsum('mf,meaf->ae', t1, eris_ovvv)
#:self.woVvO = lib.einsum('jf,mebf->mbej', t1, eris_ovvv)
#:self.woVVo = lib.einsum('jf,mfbe->mbej',-t1, eris_ovvv)
#:tau = _make_tau(t2, t1, t1)
#:self.woVoO = 0.5 * lib.einsum('mebf,ijef->mbij', eris_ovvv, tau)
#:self.woVoO += 0.5 * lib.einsum('mfbe,ijfe->mbij', eris_ovvv, tau)
eris_ovoo = np.asarray(eris.ovoo)
woVoO = np.empty((nocc,nvir,nocc,nocc), dtype=dtype)
tau = _make_tau(t2, t1, t1)
theta = t2*2 - t2.transpose(0,1,3,2)
mem_now = lib.current_memory()[0]
max_memory = max(0, self.max_memory - mem_now)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nvir**3*3))))
for seg, (p0,p1) in enumerate(lib.prange(0, nocc, blksize)):
ovvv = eris.get_ovvv(slice(p0,p1)) # ovvv = eris.ovvv[p0:p1]
# transform integrals (ia|bc) -> (ac|ib)
fswap['ebmf/%d'%seg] = np.einsum('mebf->ebmf', ovvv)
self.Fvv += np.einsum('mf,mfae->ae', t1[p0:p1], ovvv) * 2
self.Fvv -= np.einsum('mf,meaf->ae', t1[p0:p1], ovvv)
woVoO[p0:p1] = lib.einsum('mebf,ijef->mbij', ovvv, tau)
woVvO = lib.einsum('jf,mebf->mbej', t1, ovvv)
woVVo = lib.einsum('jf,mfbe->mbej',-t1, ovvv)
ovvv = None
eris_ovov = np.asarray(eris.ovov[p0:p1])
woOoV = lib.einsum('if,mfne->mnie', t1, eris_ovov)
woOoV+= eris_ovoo[:,:,p0:p1].transpose(2,0,3,1)
self.woOoV[p0:p1] = woOoV
woOoV = None
tmp = lib.einsum('njbf,mfne->mbej', t2, eris_ovov)
woVvO -= tmp * .5
woVVo += tmp
ovoo = lib.einsum('menf,jf->menj', eris_ovov, t1)
woVvO -= lib.einsum('nb,menj->mbej', t1, ovoo)
ovoo = lib.einsum('mfne,jf->menj', eris_ovov, t1)
woVVo += lib.einsum('nb,menj->mbej', t1, ovoo)
ovoo = None
ovov = eris_ovov * 2 - eris_ovov.transpose(0,3,2,1)
woVvO += lib.einsum('njfb,menf->mbej', theta, ovov) * .5
self.Fov[p0:p1] = np.einsum('nf,menf->me', t1, ovov)
tilab = np.einsum('ia,jb->ijab', t1[p0:p1], t1) * .5
tilab += t2[p0:p1]
self.Foo += lib.einsum('mief,menf->ni', tilab, ovov)
self.Fvv -= lib.einsum('mnaf,menf->ae', tilab, ovov)
eris_ovov = ovov = tilab = None
woVvO -= lib.einsum('nb,menj->mbej', t1, eris_ovoo[p0:p1,:,:])
woVVo += lib.einsum('nb,nemj->mbej', t1, eris_ovoo[:,:,p0:p1])
woVvO += np.asarray(eris.ovvo[p0:p1]).transpose(0,2,1,3)
woVVo -= np.asarray(eris.oovv[p0:p1]).transpose(0,2,3,1)
self.woVvO[p0:p1] = woVvO
self.woVVo[p0:p1] = woVVo
self.Fov += fov
self.Foo += foo + 0.5*np.einsum('me,ie->mi', self.Fov+fov, t1)
self.Fvv += fvv - 0.5*np.einsum('me,ma->ae', self.Fov+fov, t1)
# 0 or 1 virtuals
woOoO = lib.einsum('je,nemi->mnij', t1, eris_ovoo)
woOoO = woOoO + woOoO.transpose(1,0,3,2)
woOoO += np.asarray(eris.oooo).transpose(0,2,1,3)
tmp = lib.einsum('meni,jneb->mbji', eris_ovoo, t2)
woVoO -= tmp.transpose(0,1,3,2) * .5
woVoO -= tmp
tmp = None
ovoo = eris_ovoo*2 - eris_ovoo.transpose(2,1,0,3)
woVoO += lib.einsum('nemi,njeb->mbij', ovoo, theta) * .5
self.Foo += np.einsum('ne,nemi->mi', t1, ovoo)
ovoo = None
eris_ovov = np.asarray(eris.ovov)
woOoO += lib.einsum('ijef,menf->mnij', tau, eris_ovov)
self.woOoO = self.saved['woOoO'] = woOoO
woVoO -= lib.einsum('nb,mnij->mbij', t1, woOoO)
woOoO = None
tmpoovv = lib.einsum('njbf,nemf->ejmb', t2, eris_ovov)
ovov = eris_ovov*2 - eris_ovov.transpose(0,3,2,1)
eris_ovov = None
tmpovvo = lib.einsum('nifb,menf->eimb', theta, ovov)
ovov = None
tmpovvo *= -.5
tmpovvo += tmpoovv * .5
woVoO -= lib.einsum('ie,ejmb->mbij', t1, tmpovvo)
woVoO -= lib.einsum('ie,ejmb->mbji', t1, tmpoovv)
woVoO += eris_ovoo.transpose(3,1,2,0)
# 3 or 4 virtuals
eris_ovvo = np.asarray(eris.ovvo)
tmpovvo -= eris_ovvo.transpose(1,3,0,2)
fswap['ovvo'] = tmpovvo
tmpovvo = None
eris_oovv = np.asarray(eris.oovv)
tmpoovv -= eris_oovv.transpose(3,1,0,2)
fswap['oovv'] = tmpoovv
tmpoovv = None
woVoO += lib.einsum('mebj,ie->mbij', eris_ovvo, t1)
woVoO += lib.einsum('mjbe,ie->mbji', eris_oovv, t1)
woVoO += lib.einsum('me,ijeb->mbij', self.Fov, t2)
self.woVoO = self.saved['woVoO'] = woVoO
woVoO = eris_ovvo = eris_oovv = None
#:theta = t2*2 - t2.transpose(0,1,3,2)
#:eris_ovvv = lib.unpack_tril(np.asarray(eris.ovvv).reshape(nocc*nvir,nvir**2)).reshape(nocc,nvir,nvir,nvir)
#:ovvv = eris_ovvv*2 - eris_ovvv.transpose(0,3,2,1)
#:tmpab = lib.einsum('mebf,miaf->eiab', eris_ovvv, t2)
#:tmpab = tmpab + tmpab.transpose(0,1,3,2) * .5
#:tmpab-= lib.einsum('mfbe,mifa->eiba', ovvv, theta) * .5
#:self.wvOvV += eris_ovvv.transpose(2,0,3,1).conj()
#:self.wvOvV -= tmpab
nsegs = len(fswap['ebmf'])
def load_ebmf(slice):
dat = [fswap['ebmf/%d'%i][slice] for i in range(nsegs)]
return np.concatenate(dat, axis=2)
mem_now = lib.current_memory()[0]
max_memory = max(0, self.max_memory - mem_now)
blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*1e6/8/(nocc*nvir**2*4))))
for p0, p1 in lib.prange(0, nvir, blksize):
#:wvOvV = lib.einsum('mebf,miaf->eiab', ovvv, t2)
#:wvOvV += lib.einsum('mfbe,miaf->eiba', ovvv, t2)
#:wvOvV -= lib.einsum('mfbe,mifa->eiba', ovvv, t2)*2
#:wvOvV += lib.einsum('mebf,mifa->eiba', ovvv, t2)
ebmf = load_ebmf(slice(p0, p1))
wvOvV = lib.einsum('ebmf,miaf->eiab', ebmf, t2)
wvOvV = -.5 * wvOvV.transpose(0,1,3,2) - wvOvV
# Using the permutation symmetry (em|fb) = (em|bf)
efmb = load_ebmf((slice(None), slice(p0, p1)))
wvOvV += np.einsum('ebmf->bmfe', efmb.conj())
# tmp = (mf|be) - (me|bf)*.5
tmp = -.5 * ebmf
tmp += efmb.transpose(1,0,2,3)
ebmf = None
wvOvV += lib.einsum('efmb,mifa->eiba', tmp, theta)
tmp = None
wvOvV += lib.einsum('meni,mnab->eiab', eris_ovoo[:,p0:p1], tau)
wvOvV -= lib.einsum('me,miab->eiab', self.Fov[:,p0:p1], t2)
wvOvV += lib.einsum('ma,eimb->eiab', t1, fswap['ovvo'][p0:p1])
wvOvV += lib.einsum('ma,eimb->eiba', t1, fswap['oovv'][p0:p1])
self.wvOvV[p0:p1] = wvOvV
self.made_ee_imds = True
log.timer('EOM-CCSD EE intermediates', *cput0)
return self
def _make_tau(t2, t1, r1, fac=1, out=None):
tau = np.einsum('ia,jb->ijab', t1, r1)
tau = tau + tau.transpose(1,0,3,2)
tau *= fac * .5
tau += t2
return tau
