# Copyright 2014-2022 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: Samragni Banerjee <samragnibanerjee4@gmail.com>
# Alexander Sokolov <alexander.y.sokolov@gmail.com>
#
import time
import numpy as np
import pyscf.ao2mo as ao2mo
import pyscf.adc
import pyscf.adc.radc
from pyscf.adc import radc_ao2mo
import itertools
from itertools import product
from pyscf import lib
from pyscf.pbc import scf
from pyscf.pbc import df
from pyscf.pbc import mp
from pyscf.lib import logger
from pyscf.pbc.adc import kadc_rhf
from pyscf.pbc.adc import kadc_ao2mo
from pyscf.pbc.adc import dfadc
from pyscf import __config__
from pyscf.pbc.mp.kmp2 import (get_nocc, get_nmo, padding_k_idx,_padding_k_idx,
padded_mo_coeff, get_frozen_mask, _add_padding)
from pyscf.pbc.cc.kccsd_rhf import _get_epq
from pyscf.pbc.cc.kccsd_t_rhf import _get_epqr
from pyscf.pbc.lib import kpts_helper
from pyscf.lib.parameters import LOOSE_ZERO_TOL, LARGE_DENOM # noqa
from pyscf.pbc import tools
import h5py
import tempfile
[docs]
def vector_size(adc):
nkpts = adc.nkpts
nocc = adc.nocc
nvir = adc.nmo - adc.nocc
n_singles = nvir
n_doubles = nkpts * nkpts * nocc * nvir * nvir
size = n_singles + n_doubles
return size
[docs]
def get_imds(adc, eris=None):
cput0 = (time.process_time(), time.time())
log = logger.Logger(adc.stdout, adc.verbose)
if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"):
raise NotImplementedError(adc.method)
method = adc.method
nkpts = adc.nkpts
mo_energy = adc.mo_energy
mo_coeff = adc.mo_coeff
nocc = adc.nocc
nvir = adc.nmo - adc.nocc
kconserv = adc.khelper.kconserv
mo_coeff, mo_energy = _add_padding(adc, mo_coeff, mo_energy)
e_occ = [mo_energy[k][:nocc] for k in range(nkpts)]
e_vir = [mo_energy[k][nocc:] for k in range(nkpts)]
e_occ = np.array(e_occ)
e_vir = np.array(e_vir)
idn_vir = np.identity(nvir)
if eris is None:
eris = adc.transform_integrals()
eris_ovov = eris.ovov
# a-b block
# Zeroth-order terms
M_ab = np.empty((nkpts,nvir,nvir),dtype=mo_coeff.dtype)
for ka in range(nkpts):
kb = ka
M_ab[ka] = lib.einsum('ab,a->ab', idn_vir, e_vir[ka])
for kl in range(nkpts):
for km in range(nkpts):
kd = kconserv[kl,ka,km]
# Second-order terms
t2_1_mla = adc.t2[0][km,kl,ka]
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('mlad,lbmd->ab',t2_1_mla, eris_ovov[kl,kb,km],optimize=True)
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('mlad,ldmb->ab',t2_1_mla, eris_ovov[kl,kd,km],optimize=True)
t2_1_lma = adc.t2[0][kl,km,ka]
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('lmad,lbmd->ab',t2_1_lma, eris_ovov[kl,kb,km],optimize=True)
M_ab[ka] -= 0.5 * \
lib.einsum('lmad,lbmd->ab',t2_1_lma, eris_ovov[kl,kb,km],optimize=True)
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('lmad,ldmb->ab',t2_1_lma, eris_ovov[kl,kd,km],optimize=True)
del t2_1_lma
t2_1_lmb = adc.t2[0][kl,km,kb]
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('lmbd,lamd->ab',t2_1_lmb.conj(),
eris_ovov[kl,ka,km].conj(),optimize=True)
M_ab[ka] -= 0.5 * \
lib.einsum('lmbd,lamd->ab',t2_1_lmb.conj(),
eris_ovov[kl,ka,km].conj(),optimize=True)
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('lmbd,ldma->ab',t2_1_lmb.conj(),
eris_ovov[kl,kd,km].conj(),optimize=True)
del t2_1_lmb
t2_1_mlb = adc.t2[0][km,kl,kb]
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('mlbd,lamd->ab',t2_1_mlb.conj(),
eris_ovov[kl,ka,km].conj(),optimize=True)
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('mlbd,ldma->ab',t2_1_mlb.conj(),
eris_ovov[kl,kd,km].conj(),optimize=True)
del t2_1_mlb
cput0 = log.timer_debug1("Completed M_ab second-order terms ADC(2) calculation", *cput0)
if(method =='adc(3)'):
eris_oovv = eris.oovv
eris_ovvo = eris.ovvo
t1_2 = adc.t1[0]
for ka in range(nkpts):
kb = ka
for kl in range(nkpts):
if isinstance(eris.ovvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nocc:
chnk_size = nocc
a = 0
for p in range(0,nocc,chnk_size):
kd = kconserv[ka,kb,kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kd], eris.Lvv[ka,kb], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
M_ab[ka] += 2. * lib.einsum('ld,ldab->ab',t1_2[kl]
[a:a+k], eris_ovvv,optimize=True)
del eris_ovvv
kd = kconserv[kb,ka,kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kd], eris.Lvv[kb,ka], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
M_ab[ka] += 2. * lib.einsum('ld,ldba->ab',t1_2[kl]
[a:a+k].conj(), eris_ovvv.conj(),optimize=True)
del eris_ovvv
kd = kconserv[ka,kb,kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kb], eris.Lvv[ka,kd], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
M_ab[ka] -= lib.einsum('ld,lbad->ab',t1_2[kl][a:a+k],
eris_ovvv,optimize=True)
del eris_ovvv
kd = kconserv[kb,ka,kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,ka], eris.Lvv[kb,kd], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
M_ab[ka] -= lib.einsum('ld,labd->ab',t1_2[kl,a:a+k].conj(),
eris_ovvv.conj(),optimize=True)
del eris_ovvv
a += k
else :
eris_ovvv = eris.ovvv
kd = kconserv[ka,kb,kl]
M_ab[ka] += 2. * lib.einsum('ld,ldab->ab',t1_2[kl],
eris_ovvv[kl,kd,ka],optimize=True)
M_ab[ka] -= lib.einsum('ld,lbad->ab',t1_2[kl],
eris_ovvv[kl,kb,ka],optimize=True)
kd = kconserv[kb,ka,kl]
M_ab[ka] += 2. * lib.einsum('ld,ldba->ab',t1_2[kl].conj(),
eris_ovvv[kl,kd,kb].conj(),optimize=True)
M_ab[ka] -= lib.einsum('ld,labd->ab',t1_2[kl].conj(),
eris_ovvv[kl,ka,kb].conj(),optimize=True)
cput0 = log.timer_debug1("Completed M_ab ovvv ADC(3) calculation", *cput0)
for ka in range(nkpts):
kb = ka
for kl in range(nkpts):
for km in range(nkpts):
kd = kconserv[km,ka,kl]
t2_1 = adc.t2[0]
t2_2_lma = adc.t2[1][kl,km,ka]
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('lmad,lbmd->ab',t2_2_lma, eris_ovov[kl,kb,km],optimize=True)
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('lmad,ldmb->ab',t2_2_lma, eris_ovov[kl,kd,km],optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('lmad,lbmd->ab',
t2_2_lma, eris_ovov[kl,kb,km],optimize=True)
t2_2_mla = adc.t2[1][km,kl,ka]
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('mlad,lbmd->ab',t2_2_mla, eris_ovov[kl,kb,km],optimize=True)
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('mlad,ldmb->ab',t2_2_mla, eris_ovov[kl,kd,km],optimize=True)
t2_2_lmb = adc.t2[1][kl,km,kb]
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('lmbd,lamd->ab',t2_2_lmb.conj(),
eris_ovov[kl,ka,km].conj(),optimize=True)
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('lmbd,ldma->ab',t2_2_lmb.conj(),
eris_ovov[kl,kd,km].conj(),optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('lmbd,lamd->ab',t2_2_lmb.conj(),
eris_ovov[kl,ka,km].conj(),optimize=True)
t2_2_mlb = adc.t2[1][km,kl,kb]
M_ab[ka] += 0.5 * 0.5 * \
lib.einsum('mlbd,mdla->ab',t2_2_mlb.conj(),
eris_ovov[km,kd,kl].conj(),optimize=True)
M_ab[ka] -= 0.5 * 0.5 * \
lib.einsum('mlbd,ldma->ab',t2_2_mlb.conj(),
eris_ovov[kl,kd,km].conj(),optimize=True)
del t2_2_mlb
log.timer_debug1("Starting the small integrals calculation")
for kn, ke, kd in kpts_helper.loop_kkk(nkpts):
kl = kconserv[ke,kn,kd]
km = kconserv[kb,kl,kd]
temp_t2_v_1 = lib.einsum(
'lned,mlbd->nemb',t2_1[kl,kn,ke], t2_1[km,kl,kb].conj(),optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('nemb,nmae->ab',
temp_t2_v_1, eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] += 0.5 *2. * lib.einsum('nemb,neam->ab',
temp_t2_v_1, eris_ovvo[kn,ke,ka], optimize=True)
M_ab[ka] += 0.5 *2. * lib.einsum('nema,nebm->ab',
temp_t2_v_1, eris_ovvo[kn,ke,kb], optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('nema,nmbe->ab',
temp_t2_v_1, eris_oovv[kn,km,kb], optimize=True)
del temp_t2_v_1
temp_t2_v_1 = lib.einsum(
'nled,lmbd->mbne',t2_1[kn,kl,ke], t2_1[kl,km,kb].conj(),optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('mbne,nmae->ab',
temp_t2_v_1, eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] += 0.5 *2. * lib.einsum('mbne,neam->ab',
temp_t2_v_1, eris_ovvo[kn,ke,ka], optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('mane,nmbe->ab',
temp_t2_v_1, eris_oovv[kn,km,kb], optimize=True)
M_ab[ka] += 0.5 *2. * lib.einsum('mane,nebm->ab',
temp_t2_v_1, eris_ovvo[kn,ke,kb], optimize=True)
del temp_t2_v_1
temp_t2_v_2 = lib.einsum(
'nled,mlbd->nemb',t2_1[kn,kl,ke], t2_1[km,kl,kb].conj(),optimize=True)
M_ab[ka] += 0.5 * 2. * lib.einsum('nemb,nmae->ab',
temp_t2_v_2, eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] -= 0.5 * 4. * lib.einsum('nemb,neam->ab',
temp_t2_v_2, eris_ovvo[kn,ke,ka], optimize=True)
M_ab[ka] += 0.5 * 2. * lib.einsum('nema,nmbe->ab',
temp_t2_v_2, eris_oovv[kn,km,kb], optimize=True)
M_ab[ka] -= 0.5 * 4. * lib.einsum('nema,nebm->ab',
temp_t2_v_2, eris_ovvo[kn,ke,kb], optimize=True)
del temp_t2_v_2
temp_t2_v_3 = lib.einsum(
'lned,lmbd->nemb',t2_1[kl,kn,ke], t2_1[kl,km,kb].conj(),optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('nemb,neam->ab',
temp_t2_v_3, eris_ovvo[kn,ke,ka], optimize=True)
M_ab[ka] += 0.5 * 2. * lib.einsum('nemb,nmae->ab',
temp_t2_v_3, eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] += 0.5 * 2. * lib.einsum('nema,nmbe->ab',
temp_t2_v_3, eris_oovv[kn,km,kb], optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('nema,nebm->ab',
temp_t2_v_3, eris_ovvo[kn,ke,kb], optimize=True)
del temp_t2_v_3
kl = kconserv[ke,kn,kd]
km = kconserv[ke,kl,kd]
temp_t2_v_8 = lib.einsum(
'lned,mled->mn',t2_1[kl,kn,ke], t2_1[km,kl,ke].conj(),optimize=True)
M_ab[ka] += 2.* lib.einsum('mn,nmab->ab',temp_t2_v_8,
eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] -= lib.einsum('mn,nbam->ab', temp_t2_v_8,
eris_ovvo[kn,kb,ka], optimize=True)
del temp_t2_v_8
temp_t2_v_9 = lib.einsum(
'nled,mled->mn',t2_1[kn,kl,ke], t2_1[km,kl,ke].conj(),optimize=True)
M_ab[ka] -= 4. * lib.einsum('mn,nmab->ab',temp_t2_v_9,
eris_oovv[kn,km,ka], optimize=True)
M_ab[ka] += 2. * lib.einsum('mn,nbam->ab',temp_t2_v_9,
eris_ovvo[kn,kb,ka], optimize=True)
del temp_t2_v_9
for km in range(nkpts):
for kl in range(nkpts):
kf = kconserv[km,kl,ka]
temp_t2 = adc.imds.t2_1_vvvv[:]
t2_1_mla = adc.t2[0][km,kl,ka]
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('mlaf,mlbf->ab',t2_1_mla,
temp_t2[km,kl,kb].conj(), optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('mlaf,lmbf->ab',t2_1_mla,
temp_t2[kl,km,kb].conj(), optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('mlaf,mlfb->ab',t2_1_mla,
temp_t2[km,kl,kf].conj(), optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('mlaf,lmfb->ab',t2_1_mla,
temp_t2[kl,km,kf].conj(), optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('mlaf,mlbf->ab',
t2_1_mla, temp_t2[km,kl,kb].conj(), optimize=True)
del t2_1_mla
t2_1_lma = adc.t2[0][kl,km,ka]
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('lmaf,mlbf->ab',t2_1_lma,
temp_t2[km,kl,kb].conj(), optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('lmaf,lmbf->ab',t2_1_lma,
temp_t2[kl,km,kb].conj(), optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('lmaf,mlfb->ab',t2_1_lma,
temp_t2[km,kl,kf].conj(), optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('lmaf,lmfb->ab',t2_1_lma,
temp_t2[kl,km,kf].conj(), optimize=True)
del t2_1_lma
kd = kconserv[km,ka,kl]
t2_1_mlb = adc.t2[0][km,kl,kb]
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('mlad,mlbd->ab', temp_t2[km,
kl,ka].conj(), t2_1_mlb, optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('lmad,mlbd->ab', temp_t2[kl,
km,ka].conj(), t2_1_mlb, optimize=True)
M_ab[ka] -= 0.5 * lib.einsum('mlad,mlbd->ab',
temp_t2[km,kl,ka].conj(), t2_1_mlb, optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('lmad,mlbd->ab', temp_t2[kl,
km,ka].conj(), t2_1_mlb, optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('mlad,mlbd->ab', temp_t2[km,
kl,ka].conj(), t2_1_mlb, optimize=True)
del t2_1_mlb
t2_1_lmb = adc.t2[0][kl,km,kb]
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('mlad,lmbd->ab', temp_t2[km,
kl,ka].conj(), t2_1_lmb, optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('lmad,lmbd->ab', temp_t2[kl,
km,ka].conj(), t2_1_lmb, optimize=True)
M_ab[ka] -= 0.5 * 0.25* \
lib.einsum('lmad,lmbd->ab', temp_t2[kl,
km,ka].conj(), t2_1_lmb, optimize=True)
M_ab[ka] += 0.5 * 0.25* \
lib.einsum('mlad,lmbd->ab', temp_t2[km,
kl,ka].conj(), t2_1_lmb, optimize=True)
del temp_t2
del t2_1_lmb
for kd in range(nkpts):
kf = kconserv[km,kd,kl]
ke = kconserv[kb,ka,kf]
if isinstance(eris.vvvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nvir:
chnk_size = nvir
a = 0
for p in range(0,nvir,chnk_size):
eris_vvvv = dfadc.get_vvvv_df(
adc, eris.Lvv[kb,ka], eris.Lvv[ke,kf], p, chnk_size)/nkpts
k = eris_vvvv.shape[0]
M_ab[ka,a:a+k] -= lib.einsum('mldf,mled,aebf->ab',t2_1[km,kl,kd],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += lib.einsum('mldf,lmed,aebf->ab',t2_1[km,kl,kd],
t2_1[kl,km,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += lib.einsum('lmdf,mled,aebf->ab',t2_1[kl,km,kd],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] -= lib.einsum('lmdf,lmed,aebf->ab',t2_1[kl,km,kd],
t2_1[kl,km,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += 2.*lib.einsum(
'mlfd,mled,aebf->ab',t2_1[km,kl,kf], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
del eris_vvvv
eris_vvvv = dfadc.get_vvvv_df(
adc, eris.Lvv[kb,kf], eris.Lvv[ke,ka], p, chnk_size)/nkpts
M_ab[ka,a:a+k] += 0.5*lib.einsum(
'mldf,mled,aefb->ab',t2_1[km,kl,kd], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= 0.5*lib.einsum(
'mldf,lmed,aefb->ab',t2_1[km,kl,kd], t2_1[kl,km,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= 0.5*lib.einsum(
'lmdf,mled,aefb->ab',t2_1[kl,km,kd], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] += 0.5*lib.einsum(
'lmdf,lmed,aefb->ab',t2_1[kl,km,kd], t2_1[kl,km,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= lib.einsum('mlfd,mled,aefb->ab',t2_1[km,kl,kf],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
del eris_vvvv
a += k
elif isinstance(eris.vvvv, np.ndarray):
eris_vvvv = eris.vvvv
M_ab[ka] -= lib.einsum('mldf,mled,aebf->ab',t2_1[km,kl,kd],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kb], optimize=True)
M_ab[ka] += lib.einsum('mldf,lmed,aebf->ab',t2_1[km,kl,kd],
t2_1[kl,km,ke].conj(), eris_vvvv[ka,ke,kb], optimize=True)
M_ab[ka] += lib.einsum('lmdf,mled,aebf->ab',t2_1[kl,km,kd],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kb], optimize=True)
M_ab[ka] -= lib.einsum('lmdf,lmed,aebf->ab',t2_1[kl,km,kd],
t2_1[kl,km,ke].conj(), eris_vvvv[ka,ke,kb], optimize=True)
M_ab[ka] += 2.*lib.einsum('mlfd,mled,aebf->ab',t2_1[km,kl,kf],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kb], optimize=True)
M_ab[ka] += 0.5*lib.einsum('mldf,mled,aefb->ab',t2_1[km,kl,kd],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kf], optimize=True)
M_ab[ka] -= 0.5*lib.einsum('mldf,lmed,aefb->ab',t2_1[km,kl,kd],
t2_1[kl,km,ke].conj(), eris_vvvv[ka,ke,kf], optimize=True)
M_ab[ka] -= 0.5*lib.einsum('lmdf,mled,aefb->ab',t2_1[kl,km,kd],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kf], optimize=True)
M_ab[ka] += 0.5*lib.einsum('lmdf,lmed,aefb->ab',t2_1[kl,km,kd],
t2_1[kl,km,ke].conj(), eris_vvvv[ka,ke,kf], optimize=True)
M_ab[ka] -= lib.einsum('mlfd,mled,aefb->ab',t2_1[km,kl,kf],
t2_1[km,kl,ke].conj(), eris_vvvv[ka,ke,kf], optimize=True)
else :
chnk_size = adc.chnk_size
if chnk_size > nvir:
chnk_size = nvir
a = 0
for p in range(0,nvir,chnk_size):
eris_vvvv = eris.vvvv[kb,ke,ka,p:p+chnk_size]
k = eris_vvvv.shape[0]
M_ab[ka,a:a+k] -= lib.einsum('mldf,mled,aebf->ab',t2_1[km,kl,kd],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += lib.einsum('mldf,lmed,aebf->ab',t2_1[km,kl,kd],
t2_1[kl,km,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += lib.einsum('lmdf,mled,aebf->ab',t2_1[kl,km,kd],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] -= lib.einsum('lmdf,lmed,aebf->ab',t2_1[kl,km,kd],
t2_1[kl,km,ke].conj(), eris_vvvv, optimize=True)
M_ab[ka,a:a+k] += 2.*lib.einsum(
'mlfd,mled,aebf->ab',t2_1[km,kl,kf], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
del eris_vvvv
eris_vvvv = eris.vvvv[ka,ke,kf,p:p+chnk_size]
M_ab[ka,a:a+k] += 0.5*lib.einsum(
'mldf,mled,aefb->ab',t2_1[km,kl,kd], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= 0.5*lib.einsum(
'mldf,lmed,aefb->ab',t2_1[km,kl,kd], t2_1[kl,km,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= 0.5*lib.einsum(
'lmdf,mled,aefb->ab',t2_1[kl,km,kd], t2_1[km,kl,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] += 0.5*lib.einsum(
'lmdf,lmed,aefb->ab',t2_1[kl,km,kd], t2_1[kl,km,ke].conj(), eris_vvvv,
optimize=True)
M_ab[ka,a:a+k] -= lib.einsum('mlfd,mled,aefb->ab',t2_1[km,kl,kf],
t2_1[km,kl,ke].conj(), eris_vvvv, optimize=True)
del eris_vvvv
a += k
cput0 = log.timer_debug1("Completed M_ab ADC(3) calculation", *cput0)
return M_ab
[docs]
def get_diag(adc,kshift,M_ab=None,eris=None):
log = logger.Logger(adc.stdout, adc.verbose)
if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"):
raise NotImplementedError(adc.method)
if M_ab is None:
M_ab = adc.get_imds()
nkpts = adc.nkpts
kconserv = adc.khelper.kconserv
nocc = adc.nocc
nvir = adc.nmo - adc.nocc
n_singles = nvir
n_doubles = nkpts * nkpts * nocc * nvir * nvir
dim = n_singles + n_doubles
s1 = 0
f1 = n_singles
s2 = f1
f2 = s2 + n_doubles
mo_energy = adc.mo_energy
mo_coeff = adc.mo_coeff
nocc = adc.nocc
nmo = adc.nmo
nvir = nmo - nocc
mo_coeff, mo_energy = _add_padding(adc, mo_coeff, mo_energy)
e_occ = [mo_energy[k][:nocc] for k in range(nkpts)]
e_vir = [mo_energy[k][nocc:] for k in range(nkpts)]
e_vir = np.array(e_vir)
e_occ = np.array(e_occ)
diag = np.zeros((dim), dtype=np.complex128)
doubles = np.zeros((nkpts,nkpts,nocc*nvir*nvir),dtype=np.complex128)
# Compute precond in h1-h1 block
M_ab_diag = np.diagonal(M_ab[kshift])
diag[s1:f1] = M_ab_diag.copy()
# Compute precond in 2p1h-2p1h block
for kj in range(nkpts):
for ka in range(nkpts):
kb = kconserv[kshift,ka,kj]
d_ab = e_vir[ka][:,None] + e_vir[kb]
d_i = e_occ[kj][:,None]
D_n = -d_i + d_ab.reshape(-1)
doubles[kj,ka] += D_n.reshape(-1)
diag[s2:f2] = doubles.reshape(-1)
log.timer_debug1("Completed ea_diag calculation")
return diag
[docs]
def matvec(adc, kshift, M_ab=None, eris=None):
if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"):
raise NotImplementedError(adc.method)
method = adc.method
nkpts = adc.nkpts
nocc = adc.nocc
kconserv = adc.khelper.kconserv
nvir = adc.nmo - adc.nocc
n_singles = nvir
n_doubles = nkpts * nkpts * nocc * nvir * nvir
s_singles = 0
f_singles = n_singles
s_doubles = f_singles
f_doubles = s_doubles + n_doubles
mo_energy = adc.mo_energy
mo_coeff = adc.mo_coeff
mo_coeff, mo_energy = _add_padding(adc, mo_coeff, mo_energy)
e_occ = [mo_energy[k][:nocc] for k in range(nkpts)]
e_vir = [mo_energy[k][nocc:] for k in range(nkpts)]
e_vir = np.array(e_vir)
e_occ = np.array(e_occ)
if M_ab is None:
M_ab = adc.get_imds()
#Calculate sigma vector
def sigma_(r):
cput0 = (time.process_time(), time.time())
log = logger.Logger(adc.stdout, adc.verbose)
r1 = r[s_singles:f_singles]
r2 = r[s_doubles:f_doubles]
r2 = r2.reshape(nkpts,nkpts,nocc,nvir,nvir)
s2 = np.zeros((nkpts,nkpts,nocc,nvir,nvir), dtype=np.complex128)
cell = adc.cell
kpts = adc.kpts
madelung = tools.madelung(cell, kpts)
############ ADC(2) ij block ############################
s1 = lib.einsum('ab,b->a',M_ab[kshift],r1)
########### ADC(2) coupling blocks #########################
for kb in range(nkpts):
for kc in range(nkpts):
ki = kconserv[kb,kshift, kc]
if isinstance(eris.ovvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nocc:
chnk_size = nocc
a = 0
for p in range(0,nocc,chnk_size):
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[ki,kc], eris.Lvv[kshift,kb], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
s1 += 2. * lib.einsum('icab,ibc->a', eris_ovvv.conj(),
r2[ki,kb,a:a+k], optimize=True)
s2[ki,kb,a:a+k] += lib.einsum('icab,a->ibc', eris_ovvv, r1, optimize=True)
del eris_ovvv
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[ki,kb], eris.Lvv[kshift,kc], p, chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
s1 -= lib.einsum('ibac,ibc->a', eris_ovvv.conj(),
r2[ki,kb,a:a+k], optimize=True)
del eris_ovvv
a += k
else :
eris_ovvv = eris.ovvv[:]
s1 += 2. * lib.einsum('icab,ibc->a',
eris_ovvv[ki,kc,kshift].conj(), r2[ki,kb], optimize=True)
s2[ki,kb] += lib.einsum('icab,a->ibc', eris_ovvv[ki,
kc,kshift], r1, optimize=True)
del eris_ovvv
######### ########## ADC(2) ajk - bil block ############################
s2[ki, kb] -= lib.einsum('i,ibc->ibc', e_occ[ki], r2[ki, kb])
s2[ki, kb] += lib.einsum('b,ibc->ibc', e_vir[kb], r2[ki, kb])
s2[ki, kb] += lib.einsum('c,ibc->ibc', e_vir[kc], r2[ki, kb])
################ ADC(3) ajk - bil block ############################
if (method == "adc(2)-x" or method == "adc(3)"):
eris_oovv = eris.oovv
eris_ovvo = eris.ovvo
for kx in range(nkpts):
for ky in range(nkpts):
ki = kconserv[ky,kshift, kx]
for kz in range(nkpts):
kw = kconserv[kx, kz, ky]
if isinstance(eris.vvvv, np.ndarray):
eris_vvvv = eris.vvvv.reshape(nkpts,nkpts,nkpts,nvir*nvir,nvir*nvir)
r2_1 = r2.reshape(nkpts,nkpts,nocc,nvir*nvir)
s2[ki, kx] += np.dot(r2_1[ki,kw],eris_vvvv[kx,ky,
kw].T.conj()).reshape(nocc,nvir,nvir)
elif isinstance(eris.vvvv, type(None)):
s2[ki,kx] += ea_contract_r_vvvv(adc,r2[ki,kw],eris.Lvv,kx,ky,kw)
else :
s2[ki,kx] += ea_contract_r_vvvv(adc,r2[ki,kw],eris.vvvv,kx,ky,kw)
for kj in range(nkpts):
kz = kconserv[ky,ki,kj]
s2[ki,kx] -= 0.5*lib.einsum('iyzj,jzx->ixy',
eris_ovvo[ki,ky,kz],r2[kj,kz],optimize=True)
s2[ki,kx] += lib.einsum('iyzj,jxz->ixy',eris_ovvo[ki,
ky,kz],r2[kj,kx],optimize=True)
s2[ki,kx] -= 0.5*lib.einsum('ijzy,jxz->ixy',
eris_oovv[ki,kj,kz],r2[kj,kx],optimize=True)
kz = kconserv[kj,ki,kx]
s2[ki,kx] -= 0.5*lib.einsum('ijzx,jzy->ixy',
eris_oovv[ki,kj,kz],r2[kj,kz],optimize=True)
kw = kconserv[kj,ki,kx]
s2[ki,kx] -= 0.5*lib.einsum('ijwx,jwy->ixy',
eris_oovv[ki,kj,kw],r2[kj,kw],optimize=True)
kw = kconserv[kj,ki,ky]
s2[ki,kx] -= 0.5*lib.einsum('ijwy,jxw->ixy',
eris_oovv[ki,kj,kw],r2[kj,kx],optimize=True)
s2[ki,kx] += lib.einsum('iywj,jxw->ixy',eris_ovvo[ki,
ky,kw],r2[kj,kx],optimize=True)
s2[ki,kx] -= 0.5*lib.einsum('iywj,jwx->ixy',
eris_ovvo[ki,ky,kw],r2[kj,kw],optimize=True)
if adc.exxdiv is not None:
s2 += -madelung * r2
if (method == "adc(3)"):
eris_ovoo = eris.ovoo
################ ADC(3) a - ibc block and ibc-a coupling blocks ########################
t2_1 = adc.t2[0]
for kz in range(nkpts):
for kw in range(nkpts):
kj = kconserv[kw, kshift, kz]
for kl in range(nkpts):
km = kconserv[kw, kl, kz]
t2_1_lmz = adc.t2[0][kl,km,kz]
ka = kconserv[km, kj, kl]
temp_1 = lib.einsum('lmzw,jzw->jlm',t2_1_lmz,r2[kj,kz])
temp = 0.25 * lib.einsum('lmzw,jzw->jlm',t2_1_lmz,r2[kj,kz])
temp -= 0.25 * lib.einsum('lmzw,jwz->jlm',t2_1_lmz,r2[kj,kw])
del t2_1_lmz
t2_1_mlz = adc.t2[0][km,kl,kz]
temp -= 0.25 * lib.einsum('mlzw,jzw->jlm',t2_1_mlz,r2[kj,kz])
temp += 0.25 * lib.einsum('mlzw,jwz->jlm',t2_1_mlz,r2[kj,kw])
del t2_1_mlz
s1 += lib.einsum('jlm,lamj->a',temp, eris_ovoo[kl,ka,km], optimize=True)
s1 -= lib.einsum('jlm,malj->a',temp, eris_ovoo[km,ka,kl], optimize=True)
s1 += lib.einsum('jlm,lamj->a',temp_1, eris_ovoo[kl,ka,km], optimize=True)
del temp
del temp_1
for kx in range(nkpts):
for ky in range(nkpts):
ki = kconserv[ky, kshift, kx]
for kl in range(nkpts):
km = kconserv[kx, kl, ky]
kb = kconserv[km,ki,kl]
temp = lib.einsum(
'b,lbmi->lmi',r1,eris_ovoo[kl,kb,km].conj(), optimize=True)
s2[ki,kx] += lib.einsum('lmi,lmxy->ixy',temp,
t2_1[kl,km,kx].conj(), optimize=True)
del temp
for kz in range(nkpts):
for kw in range(nkpts):
kj = kconserv[kz, kshift, kw]
for kl in range(nkpts):
kd = kconserv[kj, kw, kl]
t2_1_jlw = adc.t2[0][kj,kl,kw]
temp_s_a = lib.einsum('jlwd,jzw->lzd',t2_1_jlw,r2[kj,kz],optimize=True)
temp_s_a -= lib.einsum('jlwd,jwz->lzd',t2_1_jlw,r2[kj,kw],optimize=True)
temp_t2_r2_1 = lib.einsum('jlwd,jzw->lzd',t2_1_jlw,r2[kj,kz],optimize=True)
temp_t2_r2_1 -= lib.einsum('jlwd,jwz->lzd',t2_1_jlw,r2[kj,kw],optimize=True)
temp_t2_r2_1 += lib.einsum('jlwd,jzw->lzd',t2_1_jlw,r2[kj,kz],optimize=True)
del t2_1_jlw
t2_1_ljw = adc.t2[0][kl,kj,kw]
temp_s_a -= lib.einsum('ljwd,jzw->lzd',t2_1_ljw,r2[kj,kz],optimize=True)
temp_s_a += lib.einsum('ljwd,jwz->lzd',t2_1_ljw,r2[kj,kw],optimize=True)
temp_s_a += lib.einsum('ljdw,jzw->lzd',
t2_1[kl,kj,kd],r2[kj,kz],optimize=True)
temp_t2_r2_1 -= lib.einsum('ljwd,jzw->lzd',t2_1_ljw,r2[kj,kz],optimize=True)
del t2_1_ljw
temp_t2_r2_4 = lib.einsum(
'jldw,jwz->lzd',t2_1[kj,kl,kd],r2[kj,kw], optimize=True)
if isinstance(eris.ovvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nocc:
chnk_size = nocc
a = 0
for p in range(0,nocc,chnk_size):
ka = kconserv[kz, kd, kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kd], eris.Lvv[kz,ka], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
s1 += 0.5*lib.einsum('lzd,ldza->a',
temp_s_a[a:a+k],eris_ovvv,optimize=True)
s1 += 0.5*lib.einsum('lzd,ldza->a',
temp_t2_r2_1[a:a+k],eris_ovvv,optimize=True)
del eris_ovvv
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,ka], eris.Lvv[kz,kd], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
s1 -= 0.5*lib.einsum('lzd,lazd->a',
temp_s_a[a:a+k],eris_ovvv,optimize=True)
s1 -= 0.5*lib.einsum('lzd,lazd->a',
temp_t2_r2_4[a:a+k],eris_ovvv,optimize=True)
del eris_ovvv
a += k
else :
ka = kconserv[kz, kd, kl]
eris_ovvv = eris.ovvv[:]
s1 += 0.5*lib.einsum('lzd,ldza->a',temp_s_a,
eris_ovvv[kl,kd,kz],optimize=True)
s1 += 0.5*lib.einsum('lzd,ldza->a',temp_t2_r2_1,
eris_ovvv[kl,kd,kz],optimize=True)
s1 -= 0.5*lib.einsum('lzd,lazd->a',temp_s_a,
eris_ovvv[kl,ka,kz],optimize=True)
s1 -= 0.5*lib.einsum('lzd,lazd->a',temp_t2_r2_4,
eris_ovvv[kl,ka,kz],optimize=True)
del eris_ovvv
del temp_s_a
del temp_t2_r2_1
del temp_t2_r2_4
for kz in range(nkpts):
for kw in range(nkpts):
kj = kconserv[kz, kshift, kw]
for kl in range(nkpts):
kd = kconserv[kj, kz, kl]
t2_1_jlz = adc.t2[0][kj,kl,kz]
temp_s_a_1 = -lib.einsum('jlzd,jwz->lwd',
t2_1_jlz,r2[kj,kw],optimize=True)
temp_s_a_1 += lib.einsum('jlzd,jzw->lwd',
t2_1_jlz,r2[kj,kz],optimize=True)
temp_t2_r2_2 = -lib.einsum('jlzd,jwz->lwd',
t2_1_jlz,r2[kj,kw],optimize=True)
temp_t2_r2_2 += lib.einsum('jlzd,jzw->lwd',
t2_1_jlz,r2[kj,kz],optimize=True)
temp_t2_r2_2 -= lib.einsum('jlzd,jwz->lwd',
t2_1_jlz,r2[kj,kw],optimize=True)
del t2_1_jlz
t2_1_ljz = adc.t2[0][kl,kj,kz]
temp_s_a_1 += lib.einsum('ljzd,jwz->lwd',t2_1_ljz,r2[kj,kw],optimize=True)
temp_s_a_1 -= lib.einsum('ljzd,jzw->lwd',t2_1_ljz,r2[kj,kz],optimize=True)
temp_s_a_1 -= lib.einsum('ljdz,jwz->lwd',
t2_1[kl,kj,kd],r2[kj,kw],optimize=True)
temp_t2_r2_2 += lib.einsum('ljzd,jwz->lwd',t2_1_ljz,r2[kj,kw],optimize=True)
temp_t2_r2_3 = -lib.einsum('ljzd,jzw->lwd',t2_1_ljz,r2[kj,kz],optimize=True)
del t2_1_ljz
if isinstance(eris.ovvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nocc:
chnk_size = nocc
a = 0
for p in range(0,nocc,chnk_size):
ka = kconserv[kw, kd, kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kd], eris.Lvv[kw,ka], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
s1 -= 0.5*lib.einsum('lwd,ldwa->a',
temp_s_a_1[a:a+k],eris_ovvv,optimize=True)
s1 -= 0.5*lib.einsum('lwd,ldwa->a',
temp_t2_r2_2[a:a+k],eris_ovvv,optimize=True)
del eris_ovvv
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,ka], eris.Lvv[kw,kd], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
s1 += 0.5*lib.einsum('lwd,lawd->a',
temp_s_a_1[a:a+k],eris_ovvv,optimize=True)
s1 += 0.5*lib.einsum('lwd,lawd->a',
temp_t2_r2_3[a:a+k],eris_ovvv,optimize=True)
del eris_ovvv
a += k
else :
ka = kconserv[kw, kd, kl]
eris_ovvv = eris.ovvv[:]
s1 -= 0.5*lib.einsum('lwd,ldwa->a',temp_s_a_1,
eris_ovvv[kl,kd,kw],optimize=True)
s1 += 0.5*lib.einsum('lwd,lawd->a',temp_s_a_1,
eris_ovvv[kl,ka,kw],optimize=True)
s1 -= 0.5*lib.einsum('lwd,ldwa->a',temp_t2_r2_2,
eris_ovvv[kl,kd,kw],optimize=True)
s1 += 0.5*lib.einsum('lwd,lawd->a',temp_t2_r2_3,
eris_ovvv[kl,ka,kw],optimize=True)
del eris_ovvv
del temp_s_a_1
del temp_t2_r2_2
del temp_t2_r2_3
for kx in range(nkpts):
for ky in range(nkpts):
ki = kconserv[ky,kshift,kx]
for kl in range(nkpts):
temp = np.zeros((nocc,nvir,nvir),dtype=eris.oooo.dtype)
temp_1_1 = np.zeros((nocc,nvir,nvir),dtype=eris.oooo.dtype)
temp_2_1 = np.zeros((nocc,nvir,nvir),dtype=eris.oooo.dtype)
if isinstance(eris.ovvv, type(None)):
chnk_size = adc.chnk_size
if chnk_size > nocc:
chnk_size = nocc
a = 0
for p in range(0,nocc,chnk_size):
kd = kconserv[kl, kshift, kx]
kb = kconserv[kl,kd,kx]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kd], eris.Lvv[kx,kb], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
k = eris_ovvv.shape[0]
temp_1_1[a:a+k] += lib.einsum('ldxb,b->lxd',
eris_ovvv,r1,optimize=True)
temp_2_1[a:a+k] += lib.einsum('ldxb,b->lxd',
eris_ovvv,r1,optimize=True)
del eris_ovvv
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kb], eris.Lvv[kx,kd], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
temp_1_1[a:a+k] -= lib.einsum('lbxd,b->lxd',
eris_ovvv,r1,optimize=True)
del eris_ovvv
kd = kconserv[kl, kshift, ky]
kb = kconserv[ky,kd,kl]
eris_ovvv = dfadc.get_ovvv_df(
adc, eris.Lov[kl,kb], eris.Lvv[ky,kd], p,
chnk_size).reshape(-1,nvir,nvir,nvir)/nkpts
temp[a:a+k] -= lib.einsum('lbyd,b->lyd', eris_ovvv,r1,optimize=True)
del eris_ovvv
a += k
else :
kd = kconserv[ki, ky, kl]
kb = kconserv[kl,kd,kx]
eris_ovvv = eris.ovvv[:]
temp_1_1 += lib.einsum('ldxb,b->lxd',
eris_ovvv[kl,kd,kx],r1,optimize=True)
temp_1_1 -= lib.einsum('lbxd,b->lxd',
eris_ovvv[kl,kb,kx],r1,optimize=True)
temp_2_1 += lib.einsum('ldxb,b->lxd',
eris_ovvv[kl,kd,kx],r1,optimize=True)
kd = kconserv[kl, kshift, ky]
kb = kconserv[ky,kd,kl]
temp -= lib.einsum('lbyd,b->lyd', eris_ovvv[kl,kb,ky],r1,optimize=True)
del eris_ovvv
kd = kconserv[kl,ky,ki]
s2[ki,kx] += lib.einsum('lxd,lidy->ixy',temp_1_1.conj(),
t2_1[kl,ki,kd].conj(),optimize=True)
s2[ki,kx] += lib.einsum('lxd,ilyd->ixy',temp_2_1.conj(),
t2_1[ki,kl,ky].conj(),optimize=True)
s2[ki,kx] -= lib.einsum('lxd,ildy->ixy',temp_2_1.conj(),
t2_1[ki,kl,kd].conj(),optimize=True)
kd = kconserv[kl,kx,ki]
s2[ki,kx] += lib.einsum('lyd,lixd->ixy',temp.conj(),
t2_1[kl,ki,kx].conj(),optimize=True)
del temp
del temp_1_1
del temp_2_1
s2 = s2.reshape(-1)
s = np.hstack((s1,s2))
del s1
del s2
cput0 = log.timer_debug1("completed sigma vector calculation", *cput0)
return s
return sigma_
[docs]
def get_trans_moments(adc,kshift):
nmo = adc.nmo
T = []
for orb in range(nmo):
T_a = get_trans_moments_orbital(adc,orb,kshift)
T.append(T_a)
T = np.array(T)
return T
[docs]
def get_trans_moments_orbital(adc, orb, kshift):
if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"):
raise NotImplementedError(adc.method)
method = adc.method
nkpts = adc.nkpts
nocc = adc.nocc
nvir = adc.nmo - adc.nocc
t2_1 = adc.t2[0]
idn_vir = np.identity(nvir)
T1 = np.zeros((nvir),dtype=np.complex128)
T2 = np.zeros((nkpts,nkpts,nocc,nvir,nvir), dtype=np.complex128)
######## ADC(2) 1h part ############################################
if orb < nocc:
if (adc.approx_trans_moments is False or adc.method == "adc(3)"):
t1_2 = adc.t1[0]
T1 -= t1_2[kshift][orb,:]
for kj in range(nkpts):
for ka in range(nkpts):
kb = adc.khelper.kconserv[kj, ka, kshift]
ki = adc.khelper.kconserv[ka, kj, kb]
t2_1_t= t2_1[ki,kj,ka].transpose(1,0,2,3)
T2[kj,ka] -= t2_1_t[:,orb,:,:].conj()
else:
T1 += idn_vir[(orb-nocc), :]
for kk in range(nkpts):
for kc in range(nkpts):
kl = adc.khelper.kconserv[kc, kk, kshift]
ka = adc.khelper.kconserv[kc, kl, kk]
T1 -= 0.25* \
lib.einsum('klc,klac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_1[kk,kl,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('lkc,lkac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_1[kl,kk,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('klc,klac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_1[kk,kl,ka].conj(), optimize=True)
T1 += 0.25* \
lib.einsum('lkc,klac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_1[kk,kl,ka].conj(), optimize=True)
T1 += 0.25* \
lib.einsum('klc,lkac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_1[kl,kk,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('lkc,lkac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_1[kl,kk,ka].conj(), optimize=True)
######### ADC(3) 2p-1h part ############################################
if (adc.method == "adc(2)-x" and adc.approx_trans_moments is False) or (adc.method == "adc(3)"):
t2_2 = adc.t2[1]
if orb < nocc:
for kj in range(nkpts):
for ka in range(nkpts):
kb = adc.khelper.kconserv[kj, ka, kshift]
ki = adc.khelper.kconserv[ka, kj, kb]
t2_2_t = t2_2[ki,kj,ka].conj().transpose(1,0,2,3)
T2[kj,ka] -= t2_2_t[:,orb,:,:].conj()
########### ADC(3) 1p part ############################################
if(method=='adc(3)'):
if orb < nocc:
for kk in range(nkpts):
for kc in range(nkpts):
ka = adc.khelper.kconserv[kk, kc, kshift]
T1 += 0.5*lib.einsum('kac,ck->a',
t2_1[kk,kshift,kc][:,orb,:,:], t1_2[kc].T,optimize=True)
T1 -= 0.5*lib.einsum('kac,ck->a',
t2_1[kshift,kk,ka][orb,:,:,:], t1_2[kc].T,optimize=True)
T1 -= 0.5*lib.einsum('kac,ck->a',
t2_1[kshift,kk,ka][orb,:,:,:], t1_2[kc].T,optimize=True)
else:
for kk in range(nkpts):
for kc in range(nkpts):
kl = adc.khelper.kconserv[kk, kc, kshift]
ka = adc.khelper.kconserv[kl, kc, kk]
T1 -= 0.25* \
lib.einsum('klc,klac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_2[kk,kl,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('lkc,lkac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_2[kl,kk,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('klc,klac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_2[kk,kl,ka].conj(), optimize=True)
T1 += 0.25* \
lib.einsum('klc,lkac->a',t2_1[kk,kl,kshift][:,:,
(orb-nocc),:], t2_2[kl,kk,ka].conj(), optimize=True)
T1 += 0.25* \
lib.einsum('lkc,klac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_2[kk,kl,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('lkc,lkac->a',t2_1[kl,kk,kshift][:,:,
(orb-nocc),:], t2_2[kl,kk,ka].conj(), optimize=True)
T1 -= 0.25* \
lib.einsum('klac,klc->a',t2_1[kk,kl,ka].conj(),
t2_2[kk,kl,kshift][:,:,(orb-nocc),:],optimize=True)
T1 -= 0.25* \
lib.einsum('lkac,lkc->a',t2_1[kl,kk,ka].conj(),
t2_2[kl,kk,kshift][:,:,(orb-nocc),:],optimize=True)
T1 -= 0.25* \
lib.einsum('klac,klc->a',t2_1[kk,kl,ka].conj(),
t2_2[kk,kl,kshift][:,:,(orb-nocc),:],optimize=True)
T1 += 0.25* \
lib.einsum('klac,lkc->a',t2_1[kk,kl,ka].conj(),
t2_2[kl,kk,kshift][:,:,(orb-nocc),:],optimize=True)
T1 += 0.25* \
lib.einsum('lkac,klc->a',t2_1[kl,kk,ka].conj(),
t2_2[kk,kl,kshift][:,:,(orb-nocc),:],optimize=True)
T1 -= 0.25* \
lib.einsum('lkac,lkc->a',t2_1[kl,kk,ka].conj(),
t2_2[kl,kk,kshift][:,:,(orb-nocc),:],optimize=True)
del t2_2
del t2_1
for ka in range(nkpts):
for kb in range(nkpts):
ki = adc.khelper.kconserv[kb,kshift, ka]
T2[ki,ka] += T2[ki,ka] - T2[ki,kb].transpose(0,2,1)
T2 = T2.reshape(-1)
T = np.hstack((T1,T2))
return T
[docs]
def renormalize_eigenvectors(adc, kshift, U, nroots=1):
nkpts = adc.nkpts
nocc = adc.t2[0].shape[3]
nvir = adc.nmo - adc.nocc
n_singles = nvir
for I in range(U.shape[1]):
U1 = U[:n_singles,I]
U2 = U[n_singles:,I].reshape(nkpts,nkpts,nocc,nvir,nvir)
UdotU = np.dot(U1.conj().ravel(),U1.ravel())
for ka in range(nkpts):
for kb in range(nkpts):
ki = adc.khelper.kconserv[kb,kshift, ka]
UdotU += 2.*np.dot(U2[ki,ka].conj().ravel(), U2[ki,ka].ravel()) - \
np.dot(U2[ki,ka].conj().ravel(),
U2[ki,kb].transpose(0,2,1).ravel())
U[:,I] /= np.sqrt(UdotU)
U = U.reshape(-1,nroots)
return U
[docs]
def get_properties(adc, kshift, U, nroots=1):
#Transition moments
T = adc.get_trans_moments(kshift)
#Spectroscopic amplitudes
U = adc.renormalize_eigenvectors(kshift,U,nroots)
X = np.dot(T, U).reshape(-1, nroots)
#Spectroscopic factors
P = 2.0*lib.einsum("pi,pi->i", X.conj(), X)
P = P.real
return P,X
[docs]
class RADCEA(kadc_rhf.RADC):
'''restricted ADC for EA energies and spectroscopic amplitudes
Attributes:
verbose : int
Print level. Default value equals to :class:`Mole.verbose`
max_memory : float or int
Allowed memory in MB. Default value equals to :class:`Mole.max_memory`
incore_complete : bool
Avoid all I/O. Default is False.
method : string
nth-order ADC method. Options are : ADC(2), ADC(2)-X, ADC(3). Default is ADC(2).
conv_tol : float
Convergence threshold for Davidson iterations. Default is 1e-12.
max_cycle : int
Number of Davidson iterations. Default is 50.
max_space : int
Space size to hold trial vectors for Davidson iterative diagonalization. Default is 12.
Kwargs:
nroots : int
Number of roots (eigenvalues) requested. Default value is 1.
>>> myadc = adc.RADC(mf).run()
>>> myadcip = adc.RADC(myadc).run()
Saved results
e_ea : float or list of floats
EA energy (eigenvalue). For nroots = 1, it is a single float number. If nroots > 1,
it is a list of floats for the lowest nroots eigenvalues.
v_ea : array
Eigenvectors for each EA transition.
p_ea : float
Spectroscopic amplitudes for each EA transition.
'''
_keys = set((
'tol_residual','conv_tol', 'e_corr', 'method', 'mo_coeff', 'mo_energy_b',
't1', 'mo_energy_a', 'max_space', 't2', 'max_cycle',
'kpts', 'exxdiv', 'khelper', 'cell', 'nkop_chk', 'kop_npick', 'chnk_size',
))
def __init__(self, adc):
self.verbose = adc.verbose
self.stdout = adc.stdout
self.max_memory = adc.max_memory
self.max_space = adc.max_space
self.max_cycle = adc.max_cycle
self.conv_tol = adc.conv_tol
self.tol_residual = adc.tol_residual
self.t1 = adc.t1
self.t2 = adc.t2
self.method = adc.method
self.method_type = adc.method_type
self._scf = adc._scf
self.compute_properties = adc.compute_properties
self.approx_trans_moments = adc.approx_trans_moments
self.kpts = adc._scf.kpts
self.exxdiv = adc.exxdiv
self.verbose = adc.verbose
self.max_memory = adc.max_memory
self.method = adc.method
self.khelper = adc.khelper
self.cell = adc.cell
self.mo_coeff = adc.mo_coeff
self.mo_occ = adc.mo_occ
self.frozen = adc.frozen
self._nocc = adc._nocc
self._nmo = adc._nmo
self._nvir = adc._nvir
self.nkop_chk = adc.nkop_chk
self.kop_npick = adc.kop_npick
self.t2 = adc.t2
self.e_corr = adc.e_corr
self.mo_energy = adc.mo_energy
self.imds = adc.imds
self.chnk_size = adc.chnk_size
kernel = kadc_rhf.kernel
get_imds = get_imds
get_diag = get_diag
matvec = matvec
vector_size = vector_size
get_trans_moments = get_trans_moments
renormalize_eigenvectors = renormalize_eigenvectors
get_properties = get_properties
[docs]
def get_init_guess(self, nroots=1, diag=None, ascending=True):
if diag is None :
diag = self.get_diag()
idx = None
dtype = getattr(diag, 'dtype', np.complex128)
if ascending:
idx = np.argsort(diag)
else:
idx = np.argsort(diag)[::-1]
guess = np.zeros((diag.shape[0], nroots), dtype=dtype)
min_shape = min(diag.shape[0], nroots)
guess[:min_shape,:min_shape] = np.identity(min_shape)
g = np.zeros((diag.shape[0], nroots), dtype=dtype)
g[idx] = guess.copy()
guess = []
for p in range(g.shape[1]):
guess.append(g[:,p])
return guess
[docs]
def gen_matvec(self,kshift,imds=None, eris=None):
if imds is None:
imds = self.get_imds(eris)
diag = self.get_diag(kshift,imds,eris)
matvec = self.matvec(kshift, imds, eris)
return matvec, diag
[docs]
def ea_contract_r_vvvv(adc,r2,vvvv,ka,kb,kc):
nocc = r2.shape[0]
nvir = r2.shape[1]
nkpts = adc.nkpts
kconserv = adc.khelper.kconserv
kd = kconserv[ka, kc, kb]
r2 = np.ascontiguousarray(r2.reshape(nocc,-1))
r2_vvvv = np.zeros((nvir,nvir,nocc),dtype=r2.dtype)
chnk_size = adc.chnk_size
if chnk_size > nvir:
chnk_size = nvir
a = 0
if isinstance(vvvv, np.ndarray):
vv1 = vvvv[ka,kc]
vv2 = vvvv[kb,kd]
for p in range(0,nvir,chnk_size):
vvvv_p = dfadc.get_vvvv_df(adc, vv1, vv2, p, chnk_size)/nkpts
k = vvvv_p.shape[0]
vvvv_p = vvvv_p.reshape(-1,nvir*nvir)
r2_vvvv[a:a+k] += np.dot(vvvv_p.conj(),r2.T).reshape(-1,nvir,nocc)
del vvvv_p
a += k
else :
for p in range(0,nvir,chnk_size):
vvvv_p = vvvv[ka,kb,kc][p:p+chnk_size].reshape(-1,nvir*nvir)
k = vvvv_p.shape[0]
r2_vvvv[a:a+k] += np.dot(vvvv_p.conj(),r2.T).reshape(-1,nvir,nocc)
del vvvv_p
a += k
r2_vvvv = np.ascontiguousarray(r2_vvvv.transpose(2,0,1))
return r2_vvvv