Source code for pyscf.pbc.cc.eom_kccsd_ghf

#!/usr/bin/env python
# Copyright 2017-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authors: James D. McClain
#          Mario Motta
#          Yang Gao
#          Qiming Sun <osirpt.sun@gmail.com>
#          Jason Yu
#

import itertools

import numpy as np

from pyscf import lib
from pyscf.lib import logger
from pyscf.cc import eom_rccsd
from pyscf.pbc.lib import kpts_helper
from pyscf.lib.parameters import LOOSE_ZERO_TOL, LARGE_DENOM  # noqa
from pyscf.pbc.cc import kintermediates as imd
from pyscf.pbc.cc.kccsd_rhf import _get_epq
from pyscf.pbc.cc.kccsd_t_rhf import _get_epqr
from pyscf.pbc.mp.kmp2 import (get_frozen_mask, get_nocc, get_nmo,
                               padded_mo_coeff, padding_k_idx)  # noqa

einsum = lib.einsum

[docs] def kernel(eom, nroots=1, koopmans=False, guess=None, left=False, eris=None, imds=None, partition=None, kptlist=None, dtype=None, **kwargs): '''Calculate excitation energy via eigenvalue solver Kwargs: nroots : int Number of roots (eigenvalues) requested per k-point 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. left : bool If True, calculates left eigenvectors rather than right eigenvectors. eris : `object(uccsd._ChemistsERIs)` Holds uccsd electron repulsion integrals in chemist notation. imds : `object(_IMDS)` Holds eom intermediates in chemist notation. 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). kptlist : list List of k-point indices for which eigenvalues are requested. dtype : type Type for eigenvectors. ''' 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=eris) size = eom.vector_size() nroots = min(nroots,size) nkpts = eom.nkpts if kptlist is None: kptlist = range(nkpts) # Make the max number of roots the maximum number of occupied orbitals at any given # kpoint in the list for k, kshift in enumerate(kptlist): frozen_orbs = eom.mask_frozen(np.zeros(size, dtype=int), kshift, const=1) if isinstance(frozen_orbs, tuple): nfrozen = (np.sum(frozen_orbs[0]), np.sum(frozen_orbs[1])) nroots = min(nroots, size - nfrozen[0]) nroots = min(nroots, size - nfrozen[1]) else: nfrozen = np.sum(frozen_orbs) nroots = min(nroots, size - nfrozen) if dtype is None: dtype = np.result_type(*imds.t1) evals = np.zeros((len(kptlist),nroots)) evecs = np.zeros((len(kptlist),nroots,size), dtype) convs = np.zeros((len(kptlist),nroots), dtype) for k, kshift in enumerate(kptlist): matvec, diag = eom.gen_matvec(kshift, imds, left=left, **kwargs) diag = eom.mask_frozen(diag, kshift, const=LARGE_DENOM) user_guess = False if guess is not None: user_guess = True assert len(guess) == nroots for g in guess: assert g.size == size else: user_guess = False guess = eom.get_init_guess(kshift, nroots, koopmans, diag) for ig, g in enumerate(guess): guess_norm = np.linalg.norm(g) guess_norm_tol = LOOSE_ZERO_TOL if guess_norm < guess_norm_tol: raise ValueError('Guess vector (id=%d) with norm %.4g is below threshold %.4g.\n' 'This could possibly be due to masking/freezing orbitals.\n' 'Check your guess vector to make sure it has sufficiently large norm.' % (ig, guess_norm, guess_norm_tol)) def precond(r, e0, x0): return r/(e0-diag+1e-12) eig = lib.davidson_nosym1 if user_guess or koopmans: def pickeig(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_eigenvectors=False) conv_k, evals_k, evecs_k = eig(matvec, guess, precond, pick=pickeig, tol=eom.conv_tol, max_cycle=eom.max_cycle, max_space=eom.max_space, nroots=nroots, verbose=log) else: conv_k, evals_k, evecs_k = eig(matvec, guess, precond, tol=eom.conv_tol, max_cycle=eom.max_cycle, max_space=eom.max_space, nroots=nroots, verbose=log) evals_k = evals_k.real evals[k] = evals_k evecs[k] = evecs_k convs[k] = conv_k for n, en, vn in zip(range(nroots), evals_k, evecs_k): r1, r2 = eom.vector_to_amplitudes(vn, kshift=kshift) 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 = %0.6g', n, en, qp_weight) log.timer('EOM-CCSD', *cput0) return convs, evals, evecs
[docs] def enforce_2p_spin_doublet(r2, kconserv, kshift, orbspin, excitation): '''Enforces condition that net spin can only change by +/- 1/2''' assert (excitation in ['ip', 'ea']) if excitation == 'ip': nkpts, nocc, nvir = np.array(r2.shape)[[1, 3, 4]] elif excitation == 'ea': nkpts, nocc, nvir = np.array(r2.shape)[[1, 2, 3]] else: raise NotImplementedError idxoa = [np.where(orbspin[k][:nocc] == 0)[0] for k in range(nkpts)] idxob = [np.where(orbspin[k][:nocc] == 1)[0] for k in range(nkpts)] idxva = [np.where(orbspin[k][nocc:] == 0)[0] for k in range(nkpts)] idxvb = [np.where(orbspin[k][nocc:] == 1)[0] for k in range(nkpts)] if excitation == 'ip': for ki, kj in itertools.product(range(nkpts), repeat=2): if ki > kj: # Avoid double-counting of anti-symmetrization continue ka = kconserv[ki, kshift, kj] idxoaa = idxoa[ki][:,None] * nocc + idxoa[kj] #idxoab = idxoa[ki][:,None] * nocc + idxob[kj] #idxoba = idxob[ki][:,None] * nocc + idxoa[kj] idxobb = idxob[ki][:,None] * nocc + idxob[kj] r2_tmp = 0.5 * (r2[ki, kj] - r2[kj, ki].transpose(1, 0, 2)) r2_tmp = r2_tmp.reshape(nocc**2, nvir) # Zero out states with +/- 3 unpaired spins r2_tmp[idxobb.ravel()[:, None], idxva[ka]] = 0.0 r2_tmp[idxoaa.ravel()[:, None], idxvb[ka]] = 0.0 r2[ki, kj] = r2_tmp.reshape(nocc, nocc, nvir) r2[kj, ki] = -r2[ki, kj].transpose(1, 0, 2) # Enforce antisymmetry else: for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift, ka, kj] if ka > kb: # Avoid double-counting of anti-symmetrization continue idxvaa = idxva[ka][:,None] * nvir + idxva[kb] #idxvab = idxva[ka][:,None] * nvir + idxvb[kb] #idxvba = idxvb[ka][:,None] * nvir + idxva[kb] idxvbb = idxvb[ka][:,None] * nvir + idxvb[kb] r2_tmp = 0.5 * (r2[kj, ka] - r2[kj, kb].transpose(0, 2, 1)) r2_tmp = r2_tmp.reshape(nocc, nvir**2) # Zero out states with +/- 3 unpaired spins r2_tmp[idxoa[kshift], idxvbb.ravel()[:, None]] = 0.0 r2_tmp[idxob[kshift], idxvaa.ravel()[:, None]] = 0.0 r2[kj, ka] = r2_tmp.reshape(nocc, nvir, nvir) r2[kj, kb] = -r2[kj, ka].transpose(0, 2, 1) # Enforce antisymmetry return r2
[docs] def get_padding_k_idx(eom, cc): return padding_k_idx(cc, kind="split")
######################################## # EOM-IP-CCSD ########################################
[docs] def enforce_2p_spin_ip_doublet(r2, kconserv, kshift, orbspin): return enforce_2p_spin_doublet(r2, kconserv, kshift, orbspin, 'ip')
[docs] def spin2spatial_ip_doublet(r1, r2, kconserv, kshift, orbspin): '''Convert R1/R2 of spin orbital representation to R1/R2 of spatial orbital representation ''' nkpts, nocc, nvir = np.array(r2.shape)[[1, 3, 4]] idxoa = [np.where(orbspin[k][:nocc] == 0)[0] for k in range(nkpts)] idxob = [np.where(orbspin[k][:nocc] == 1)[0] for k in range(nkpts)] idxva = [np.where(orbspin[k][nocc:] == 0)[0] for k in range(nkpts)] idxvb = [np.where(orbspin[k][nocc:] == 1)[0] for k in range(nkpts)] nocc_a = len(idxoa[0]) # Assume nocc/nvir same for each k-point nocc_b = len(idxob[0]) nvir_a = len(idxva[0]) nvir_b = len(idxvb[0]) r1a = r1[idxoa[kshift]] r1b = r1[idxob[kshift]] r2aaa = np.zeros((nkpts,nkpts,nocc_a,nocc_a,nvir_a), dtype=r2.dtype) r2baa = np.zeros((nkpts,nkpts,nocc_b,nocc_a,nvir_a), dtype=r2.dtype) r2abb = np.zeros((nkpts,nkpts,nocc_a,nocc_b,nvir_b), dtype=r2.dtype) r2bbb = np.zeros((nkpts,nkpts,nocc_b,nocc_b,nvir_b), dtype=r2.dtype) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki, kshift, kj] idxoaa = idxoa[ki][:,None] * nocc + idxoa[kj] idxoab = idxoa[ki][:,None] * nocc + idxob[kj] idxoba = idxob[ki][:,None] * nocc + idxoa[kj] idxobb = idxob[ki][:,None] * nocc + idxob[kj] r2_tmp = r2[ki, kj].reshape(nocc**2, nvir) r2aaa_tmp = lib.take_2d(r2_tmp, idxoaa.ravel(), idxva[ka]) r2baa_tmp = lib.take_2d(r2_tmp, idxoba.ravel(), idxva[ka]) r2abb_tmp = lib.take_2d(r2_tmp, idxoab.ravel(), idxvb[ka]) r2bbb_tmp = lib.take_2d(r2_tmp, idxobb.ravel(), idxvb[ka]) r2aaa[ki, kj] = r2aaa_tmp.reshape(nocc_a, nocc_a, nvir_a) r2baa[ki, kj] = r2baa_tmp.reshape(nocc_b, nocc_a, nvir_a) r2abb[ki, kj] = r2abb_tmp.reshape(nocc_a, nocc_b, nvir_b) r2bbb[ki, kj] = r2bbb_tmp.reshape(nocc_b, nocc_b, nvir_b) return [r1a, r1b], [r2aaa, r2baa, r2abb, r2bbb]
[docs] def spatial2spin_ip_doublet(r1, r2, kconserv, kshift, orbspin=None): '''Convert R1/R2 of spatial orbital representation to R1/R2 of spin orbital representation ''' r1a, r1b = r1 r2aaa, r2baa, r2abb, r2bbb = r2 nkpts, nocc_a, nvir_a = np.array(r2aaa.shape)[[1, 3, 4]] nkpts, nocc_b, nvir_b = np.array(r2bbb.shape)[[1, 3, 4]] if orbspin is None: orbspin = np.zeros((nkpts, nocc_a+nocc_b+nvir_a+nvir_b), dtype=int) orbspin[:,1::2] = 1 nocc = nocc_a + nocc_b nvir = nvir_a + nvir_b idxoa = [np.where(orbspin[k][:nocc] == 0)[0] for k in range(nkpts)] idxob = [np.where(orbspin[k][:nocc] == 1)[0] for k in range(nkpts)] idxva = [np.where(orbspin[k][nocc:] == 0)[0] for k in range(nkpts)] idxvb = [np.where(orbspin[k][nocc:] == 1)[0] for k in range(nkpts)] r1 = np.zeros(nocc, dtype = r1a.dtype) r1[idxoa[kshift]] = r1a r1[idxob[kshift]] = r1b r2 = np.zeros((nkpts, nkpts, nocc**2, nvir), dtype = r2aaa.dtype) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki, kshift, kj] idxoaa = idxoa[ki][:,None] * nocc + idxoa[kj] idxoab = idxoa[ki][:,None] * nocc + idxob[kj] idxoba = idxob[ki][:,None] * nocc + idxoa[kj] idxobb = idxob[ki][:,None] * nocc + idxob[kj] r2aaa_tmp = r2aaa[ki,kj].reshape(nocc_a * nocc_a, nvir_a) r2baa_tmp = r2baa[ki,kj].reshape(nocc_b * nocc_a, nvir_a) r2abb_tmp = r2abb[ki,kj].reshape(nocc_a * nocc_b, nvir_b) r2bbb_tmp = r2bbb[ki,kj].reshape(nocc_b * nocc_b, nvir_b) lib.takebak_2d(r2[ki, kj], r2aaa_tmp, idxoaa.ravel(), idxva[ka]) lib.takebak_2d(r2[ki, kj], r2baa_tmp, idxoba.ravel(), idxva[ka]) lib.takebak_2d(r2[ki, kj], r2abb_tmp, idxoab.ravel(), idxvb[ka]) lib.takebak_2d(r2[ki, kj], r2bbb_tmp, idxobb.ravel(), idxvb[ka]) r2aba_tmp = - r2baa[kj,ki].reshape(nocc_a * nocc_b, nvir_a) r2bab_tmp = - r2abb[kj,ki].reshape(nocc_a * nocc_b, nvir_b) lib.takebak_2d(r2[ki, kj], r2aba_tmp, idxoab.T.ravel(), idxva[ka]) lib.takebak_2d(r2[ki, kj], r2bab_tmp, idxoba.T.ravel(), idxvb[ka]) r2 = r2.reshape(nkpts, nkpts, nocc, nocc, nvir) return r1, r2
[docs] def vector_to_amplitudes_ip(vector, kshift, nkpts, nmo, nocc, kconserv): nvir = nmo - nocc r1 = vector[:nocc].copy() r2_tril = vector[nocc:].copy().reshape(nkpts*nocc*(nkpts*nocc-1)//2,nvir) idx, idy = np.tril_indices(nkpts*nocc, -1) r2 = np.zeros((nkpts*nocc,nkpts*nocc,nvir), dtype=vector.dtype) r2[idx, idy] = r2_tril r2[idy, idx] = -r2_tril r2 = r2.reshape(nkpts,nocc,nkpts,nocc,nvir).transpose(0,2,1,3,4) return [r1,r2]
[docs] def amplitudes_to_vector_ip(r1, r2, kshift, kconserv): nkpts, nocc, nvir = np.asarray(r2.shape)[[0,2,4]] # From symmetry for aaa and bbb terms, only store lower # triangular part (ki,i) < (kj,j) idx, idy = np.tril_indices(nkpts*nocc, -1) r2 = r2.transpose(0,2,1,3,4).reshape(nkpts*nocc,nkpts*nocc,nvir) return np.hstack((r1, r2[idx,idy].ravel()))
[docs] def ipccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''2ph operators are of the form s_{ij}^{a }, i.e. 'ia' indices are coupled. This differs from the restricted case that uses s_{ij}^{ b}.''' if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nkpts = eom.nkpts kconserv = imds.kconserv r1, r2 = vector_to_amplitudes_ip(vector, kshift, nkpts, nmo, nocc, kconserv) Hr1 = -np.einsum('mi,m->i', imds.Foo[kshift], r1) for km in range(nkpts): Hr1 += np.einsum('me,mie->i', imds.Fov[km], r2[km, kshift]) for kn in range(nkpts): Hr1 += - 0.5 * np.einsum('nmie,mne->i', imds.Wooov[kn, km, kshift], r2[km, kn]) Hr2 = np.zeros_like(r2) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki, kshift, kj] Hr2[ki, kj] += lib.einsum('ae,ije->ija', imds.Fvv[ka], r2[ki, kj]) Hr2[ki, kj] -= lib.einsum('mi,mja->ija', imds.Foo[ki], r2[ki, kj]) Hr2[ki, kj] += lib.einsum('mj,mia->ija', imds.Foo[kj], r2[kj, ki]) Hr2[ki, kj] -= np.einsum('maji,m->ija', imds.Wovoo[kshift, ka, kj], r1) for km in range(nkpts): kn = kconserv[ki, km, kj] Hr2[ki, kj] += 0.5 * lib.einsum('mnij,mna->ija', imds.Woooo[km, kn, ki], r2[km, kn]) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki, kshift, kj] for km in range(nkpts): ke = kconserv[km, kshift, kj] Hr2[ki, kj] += lib.einsum('maei,mje->ija', imds.Wovvo[km, ka, ke], r2[km, kj]) ke = kconserv[km, kshift, ki] Hr2[ki, kj] -= lib.einsum('maej,mie->ija', imds.Wovvo[km, ka, ke], r2[km, ki]) tmp = lib.einsum('xymnef,xymnf->e', imds.Woovv[:, :, kshift], r2[:, :]) # contract_{km, kn} Hr2[:, :] += 0.5 * lib.einsum('e,yxjiea->xyija', tmp, imds.t2[:, :, kshift]) # sum_{ki, kj} vector = amplitudes_to_vector_ip(Hr1, Hr2, kshift, kconserv) return vector
[docs] def lipccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''2ph operators are of the form s_{ij}^{ b}, i.e. 'jb' indices are coupled. See also `ipccsd_matvec`''' if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nvir = nmo - nocc nkpts = eom.nkpts kconserv = imds.kconserv r1, r2 = vector_to_amplitudes_ip(vector, kshift, nkpts, nmo, nocc, kconserv) dtype = np.result_type(r1, r2) Hr1 = -lib.einsum('mi,i->m', imds.Foo[kshift], r1) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki, kshift, kj] Hr1 += -0.5 * lib.einsum('maji,ija->m', imds.Wovoo[kshift,ka,kj], r2[ki,kj]) Hr2 = np.zeros_like(r2) for km, kn in itertools.product(range(nkpts), repeat=2): ke = kconserv[km, kshift, kn] Hr2[km,kn] += -lib.einsum('nmie,i->mne', imds.Wooov[kn,km,kshift], r1) Hr2[km,kshift] += (km==ke)*lib.einsum('me,n->mne', imds.Fov[km], r1) Hr2[kshift,kn] -= (kn==ke)*lib.einsum('ne,m->mne', imds.Fov[kn], r1) for km, kn in itertools.product(range(nkpts), repeat=2): ke = kconserv[km, kshift, kn] Hr2[km,kn] += lib.einsum('ae,mna->mne', imds.Fvv[ke], r2[km,kn]) tmp1 = lib.einsum('mi,ine->mne', imds.Foo[km], r2[km,kn]) tmp1T = lib.einsum('ni,ime->mne', imds.Foo[kn], r2[kn,km]) Hr2[km,kn] += (-tmp1 + tmp1T) for ki in range(nkpts): kj = kconserv[km,ki,kn] Hr2[km,kn] += 0.5 * lib.einsum('mnij,ije->mne', imds.Woooo[km,kn,ki], r2[ki,kj]) ka = kconserv[ke,km,ki] tmp2 = lib.einsum('maei,ina->mne', imds.Wovvo[km,ka,ke], r2[ki,kn]) ka = kconserv[ke,kn,ki] tmp2T = lib.einsum('naei,ima->mne', imds.Wovvo[kn,ka,ke], r2[ki,km]) Hr2[km,kn] += (tmp2 - tmp2T) tmp = np.zeros(nvir, dtype=dtype) for ki, kj in itertools.product(range(nkpts), repeat=2): ka = kconserv[ki,kshift,kj] kf = kshift tmp += lib.einsum('ija,ijaf->f',r2[ki,kj],imds.t2[ki,kj,ka]) for km, kn in itertools.product(range(nkpts), repeat=2): ke = kconserv[km, kshift, kn] Hr2[km,kn] += 0.5 * lib.einsum('mnfe,f->mne', imds.Woovv[km,kn,kf], tmp) vector = amplitudes_to_vector_ip(Hr1, Hr2, kshift, kconserv) return vector
[docs] def ipccsd_diag(eom, kshift, imds=None): if imds is None: imds = eom.make_imds() t1 = imds.t1 nkpts, nocc, nvir = t1.shape kconserv = imds.kconserv Hr1 = -np.diag(imds.Foo[kshift]) Hr2 = np.zeros((nkpts,nkpts,nocc,nocc,nvir), dtype=t1.dtype) if eom.partition == 'mp': foo = eom.eris.fock[:,:nocc,:nocc] fvv = eom.eris.fock[:,nocc:,nocc:] for ki in range(nkpts): for kj in range(nkpts): ka = kconserv[ki,kshift,kj] Hr2[ki,kj] -= foo[ki].diagonal()[:,None,None] Hr2[ki,kj] -= foo[kj].diagonal()[None,:,None] Hr2[ki,kj] += fvv[ka].diagonal()[None,None,:] else: for ki in range(nkpts): for kj in range(nkpts): ka = kconserv[ki,kshift,kj] Hr2[ki,kj] -= imds.Foo[ki].diagonal()[:,None,None] Hr2[ki,kj] -= imds.Foo[kj].diagonal()[None,:,None] Hr2[ki,kj] += imds.Fvv[ka].diagonal()[None,None,:] if ki == kconserv[ki,kj,kj]: Hr2[ki,kj] += np.einsum('ijij->ij', imds.Woooo[ki, kj, ki])[:,:,None] Hr2[ki, kj] += lib.einsum('iaai->ia', imds.Wovvo[ki, ka, ka])[:,None,:] Hr2[ki, kj] += lib.einsum('jaaj->ja', imds.Wovvo[kj, ka, ka])[None,:,:] Hr2[ki, kj] += lib.einsum('ijea,jiea->ija',imds.Woovv[ki,kj,kshift], imds.t2[kj,ki,kshift]) vector = amplitudes_to_vector_ip(Hr1, Hr2, kshift, kconserv) return vector
[docs] def ipccsd_star_contract(eom, ipccsd_evals, ipccsd_evecs, lipccsd_evecs, kshift, imds=None): """ Returns: e_star (list of float): The IP-CCSD* energy. Notes: The user should check to make sure the right and left eigenvalues before running the perturbative correction. The 2hp right amplitudes are assumed to be of the form s^{a }_{ij}, i.e. the (ia) indices are coupled. Reference: Saeh, Stanton "...energy surfaces of radicals" JCP 111, 8275 (1999); DOI:10.1063/1.480171 """ assert (eom.partition is None) if imds is None: imds = eom.make_imds() t1, t2 = imds.t1, imds.t2 eris = imds.eris nkpts, nocc, nvir = t1.shape nmo = nocc + nvir dtype = np.result_type(t1, t2) kconserv = eom.kconserv mo_energy_occ = np.array([eris.mo_energy[ki][:nocc] for ki in range(nkpts)]) mo_energy_vir = np.array([eris.mo_energy[ki][nocc:] for ki in range(nkpts)]) mo_e_o = mo_energy_occ mo_e_v = mo_energy_vir ipccsd_evecs = np.array(ipccsd_evecs) lipccsd_evecs = np.array(lipccsd_evecs) e_star = [] ipccsd_evecs, lipccsd_evecs = [np.atleast_2d(x) for x in [ipccsd_evecs, lipccsd_evecs]] ipccsd_evals = np.atleast_1d(ipccsd_evals) for ip_eval, ip_evec, ip_levec in zip(ipccsd_evals, ipccsd_evecs, lipccsd_evecs): # Enforcing <L|R> = 1 l1, l2 = vector_to_amplitudes_ip(ip_levec, kshift, nkpts, nmo, nocc, kconserv) r1, r2 = vector_to_amplitudes_ip(ip_evec, kshift, nkpts, nmo, nocc, kconserv) ldotr = np.dot(l1, r1) + 0.5 * np.dot(l2.ravel(), r2.ravel()) logger.info(eom, 'Left-right amplitude overlap : %14.8e + 1j %14.8e', ldotr.real, ldotr.imag) if abs(ldotr) < 1e-7: logger.warn(eom, 'Small %s left-right amplitude overlap. Results ' 'may be inaccurate.', ldotr) l1 /= ldotr l2 /= ldotr deltaE = 0.0 + 1j*0.0 for ka, kb in itertools.product(range(nkpts), repeat=2): lijkab = np.zeros((nkpts,nkpts,nocc,nocc,nocc,nvir,nvir),dtype=dtype) rijkab = np.zeros((nkpts,nkpts,nocc,nocc,nocc,nvir,nvir),dtype=dtype) kklist = kpts_helper.get_kconserv3(eom._cc._scf.cell, eom._cc.kpts, [ka,kb,kshift,range(nkpts),range(nkpts)]) for ki, kj in itertools.product(range(nkpts), repeat=2): kk = kklist[ki,kj] #TODO: can reduce size of ijkab arrays since `kk` fixed from other k-points # lijkab update if kk == kshift and kb == kconserv[ki,ka,kj]: lijkab[ki,kj] += lib.einsum('ijab,k->ijkab', eris.oovv[ki,kj,ka], l1) km = kconserv[kj,ka,ki] tmp = lib.einsum('jima,mkb->ijkab', eris.ooov[kj,ki,km], l2[km,kk]) km = kconserv[kj,kb,ki] tmpT = lib.einsum('jimb,mka->ijkab', eris.ooov[kj,ki,km], l2[km,kk]) lijkab[ki,kj] += (-tmp + tmpT) ke = kconserv[ka,ki,kb] lijkab[ki,kj] += lib.einsum('ieab,jke->ijkab', eris.ovvv[ki,ke,ka], l2[kj,kk]) # rijkab update tmp = lib.einsum('mbke,m->bke', eris.ovov[kshift,kb,kk], r1) tmp = lib.einsum('bke,ijae->ijkab', tmp, t2[ki,kj,ka]) tmpT = lib.einsum('make,m->ake', eris.ovov[kshift,ka,kk], r1) tmpT = lib.einsum('ake,ijbe->ijkab', tmpT, t2[ki,kj,kb]) rijkab[ki,kj] -= (tmp - tmpT) km = kconserv[kj,kshift,kk] tmp = lib.einsum('mnjk,n->mjk', eris.oooo[km,kshift,kj], r1) tmp = lib.einsum('mjk,imab->ijkab', tmp, t2[ki,km,ka]) rijkab[ki,kj] += tmp km = kconserv[kj,ka,ki] tmp = lib.einsum('jima,mkb->ijkab', eris.ooov[kj,ki,km].conj(), r2[km,kk]) km = kconserv[kj,kb,ki] tmpT = lib.einsum('jimb,mka->ijkab', eris.ooov[kj,ki,km].conj(), r2[km,kk]) rijkab[ki,kj] -= (tmp - tmpT) ke = kconserv[ka,ki,kb] rijkab[ki,kj] += lib.einsum('ieab,jke->ijkab', eris.ovvv[ki,ke,ka].conj(), r2[kj,kk]) eijk = np.zeros((nkpts,nkpts,nocc,nocc,nocc), dtype=dtype) Plijkab = np.zeros_like(lijkab) Prijkab = np.zeros_like(rijkab) for ki, kj in itertools.product(range(nkpts), repeat=2): kk = kklist[ki,kj] # P(ijk) Plijkab[ki,kj] = (lijkab[ki,kj] + lijkab[kj,kk].transpose(2,0,1,3,4) + lijkab[kk,ki].transpose(1,2,0,3,4)) Prijkab[ki,kj] = (rijkab[ki,kj] + rijkab[kj,kk].transpose(2,0,1,3,4) + rijkab[kk,ki].transpose(1,2,0,3,4)) eijk[ki,kj] = _get_epqr([0,nocc,ki,mo_e_o,eom.nonzero_opadding], [0,nocc,kj,mo_e_o,eom.nonzero_opadding], [0,nocc,kk,mo_e_o,eom.nonzero_opadding]) # Creating denominator eab = _get_epq([0,nvir,ka,mo_e_v,eom.nonzero_vpadding], [0,nvir,kb,mo_e_v,eom.nonzero_vpadding], fac=[-1.,-1.]) eijkab = (eijk[:, :, :, :, :, None, None] + eab[None, None, None, None, None, :, :]) denom = eijkab + ip_eval denom = 1. / denom deltaE += lib.einsum('xyijkab,xyijkab,xyijkab', Plijkab, Prijkab, denom) deltaE *= 1./12 deltaE = deltaE.real logger.info(eom, "Exc. energy, delta energy = %16.12f, %16.12f", ip_eval + deltaE, deltaE) e_star.append(ip_eval + deltaE) return e_star
[docs] def ipccsd(eom, nroots=1, koopmans=False, guess=None, left=False, eris=None, imds=None, partition=None, kptlist=None, dtype=None, **kwargs): '''See `kernel()` for a description of arguments.''' if partition: eom.partition = partition.lower() assert eom.partition in ['mp','full'] if eom.partition in ['mp', 'full']: raise NotImplementedError eom.converged, eom.e, eom.v \ = kernel(eom, nroots, koopmans, guess, left, eris=eris, imds=imds, partition=partition, kptlist=kptlist, dtype=dtype) return eom.e, eom.v
[docs] def perturbed_ccsd_kernel(eom, nroots=1, koopmans=False, right_guess=None, left_guess=None, eris=None, imds=None, partition=None, kptlist=None, dtype=None): '''Wrapper for running perturbative excited-states that require both left and right amplitudes.''' from pyscf.cc.eom_rccsd import _sort_left_right_eigensystem if imds is None: imds = eom.make_imds(eris=eris) e_star = [] for k, kshift in enumerate(kptlist): # Right eigenvectors r_converged, r_e, r_v = \ kernel(eom, nroots, koopmans=koopmans, guess=right_guess, left=False, eris=eris, imds=imds, partition=partition, kptlist=[kshift,], dtype=dtype) # Left eigenvectors l_converged, l_e, l_v = \ kernel(eom, nroots, koopmans=koopmans, guess=right_guess, left=True, eris=eris, imds=imds, partition=partition, kptlist=[kshift,], dtype=dtype) ek, r_vk, l_vk = _sort_left_right_eigensystem(eom, r_converged[0], r_e[0], r_v[0], l_converged[0], l_e[0], l_v[0]) e_star.append(eom.ccsd_star_contract(ek, r_vk, l_vk, kshift, imds=imds)) return e_star
[docs] def ipccsd_star(eom, nroots=1, koopmans=False, right_guess=None, left_guess=None, eris=None, imds=None, partition=None, kptlist=None, dtype=None, **kwargs): '''See `kernel()` for a description of arguments.''' if partition: raise NotImplementedError return perturbed_ccsd_kernel(eom, nroots=nroots, koopmans=koopmans, right_guess=right_guess, left_guess=left_guess, eris=eris, imds=imds, partition=partition, kptlist=kptlist, dtype=dtype)
[docs] def mask_frozen_ip(eom, vector, kshift, const=LARGE_DENOM): '''Replaces all frozen orbital indices of `vector` with the value `const`.''' r1, r2 = eom.vector_to_amplitudes(vector, kshift=kshift) nkpts = eom.nkpts kconserv = eom.kconserv # Get location of padded elements in occupied and virtual space nonzero_opadding, nonzero_vpadding = eom.nonzero_opadding, eom.nonzero_vpadding new_r1 = const * np.ones_like(r1) new_r2 = const * np.ones_like(r2) new_r1[nonzero_opadding[kshift]] = r1[nonzero_opadding[kshift]] for ki in range(nkpts): for kj in range(nkpts): kb = kconserv[ki, kshift, kj] idx = np.ix_([ki], [kj], nonzero_opadding[ki], nonzero_opadding[kj], nonzero_vpadding[kb]) new_r2[idx] = r2[idx] return eom.amplitudes_to_vector(new_r1, new_r2, kshift, kconserv)
[docs] class EOMIP(eom_rccsd.EOMIP): def __init__(self, cc): self.kpts = cc.kpts self.nonzero_opadding, self.nonzero_vpadding = self.get_padding_k_idx(cc) self.kconserv = cc.khelper.kconserv eom_rccsd.EOM.__init__(self, cc) kernel = ipccsd ipccsd = ipccsd ipccsd_star = ipccsd_star ccsd_star_contract = ipccsd_star_contract get_diag = ipccsd_diag matvec = ipccsd_matvec l_matvec = lipccsd_matvec mask_frozen = mask_frozen_ip get_padding_k_idx = get_padding_k_idx
[docs] def ipccsd_star_contract(self, ipccsd_evals, ipccsd_evecs, lipccsd_evecs, kshift, imds=None): return self.ccsd_star_contract(ipccsd_evals, ipccsd_evecs, lipccsd_evecs, kshift, imds=imds)
[docs] def get_init_guess(self, kshift, nroots=1, koopmans=False, diag=None): size = self.vector_size() dtype = getattr(diag, 'dtype', np.complex128) nroots = min(nroots, size) guess = [] if koopmans: for n in self.nonzero_opadding[kshift][::-1][:nroots]: g = np.zeros(int(size), dtype=dtype) g[n] = 1.0 g = self.mask_frozen(g, kshift, const=0.0) guess.append(g) else: idx = diag.argsort()[:nroots] for i in idx: g = np.zeros(int(size), dtype=dtype) g[i] = 1.0 g = self.mask_frozen(g, kshift, const=0.0) guess.append(g) return guess
@property def nkpts(self): return len(self.kpts)
[docs] def gen_matvec(self, kshift, imds=None, left=False, **kwargs): if imds is None: imds = self.make_imds() diag = self.get_diag(kshift, imds) if left: matvec = lambda xs: [self.l_matvec(x, kshift, imds, diag) for x in xs] else: matvec = lambda xs: [self.matvec(x, kshift, imds, diag) for x in xs] return matvec, diag
[docs] def vector_to_amplitudes(self, vector, kshift=None, nkpts=None, nmo=None, nocc=None, kconserv=None): if nmo is None: nmo = self.nmo if nocc is None: nocc = self.nocc if nkpts is None: nkpts = self.nkpts if kconserv is None: kconserv = self.kconserv return vector_to_amplitudes_ip(vector, kshift, nkpts, nmo, nocc, kconserv)
[docs] def amplitudes_to_vector(self, r1, r2, kshift, kconserv=None): if kconserv is None: kconserv = self.kconserv return amplitudes_to_vector_ip(r1, r2, kshift, kconserv)
[docs] def vector_size(self): nocc = self.nocc nvir = self.nmo - nocc nkpts = self.nkpts return nocc + nkpts*nocc*(nkpts*nocc-1)*nvir//2
[docs] def make_imds(self, eris=None): imds = _IMDS(self._cc, eris=eris) imds.make_ip() return imds
[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) return imds
######################################## # EOM-EA-CCSD ########################################
[docs] def enforce_2p_spin_ea_doublet(r2, kconserv, kshift, orbspin): return enforce_2p_spin_doublet(r2, kconserv, kshift, orbspin, 'ea')
[docs] def spin2spatial_ea_doublet(r1, r2, kconserv, kshift, orbspin): '''Convert R1/R2 of spin orbital representation to R1/R2 of spatial orbital representation''' nkpts, nocc, nvir = np.array(r2.shape)[[1, 2, 3]] idxoa = [np.where(orbspin[k][:nocc] == 0)[0] for k in range(nkpts)] idxob = [np.where(orbspin[k][:nocc] == 1)[0] for k in range(nkpts)] idxva = [np.where(orbspin[k][nocc:] == 0)[0] for k in range(nkpts)] idxvb = [np.where(orbspin[k][nocc:] == 1)[0] for k in range(nkpts)] nocc_a = len(idxoa[0]) nocc_b = len(idxob[0]) nvir_a = len(idxva[0]) nvir_b = len(idxvb[0]) r1a = r1[idxva[kshift]] r1b = r1[idxvb[kshift]] r2aaa = np.zeros((nkpts,nkpts,nocc_a,nvir_a,nvir_a), dtype=r2.dtype) r2aba = np.zeros((nkpts,nkpts,nocc_a,nvir_b,nvir_a), dtype=r2.dtype) r2bab = np.zeros((nkpts,nkpts,nocc_b,nvir_a,nvir_b), dtype=r2.dtype) r2bbb = np.zeros((nkpts,nkpts,nocc_b,nvir_b,nvir_b), dtype=r2.dtype) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift, ka, kj] idxvaa = idxva[ka][:,None] * nvir + idxva[kb] idxvab = idxva[ka][:,None] * nvir + idxvb[kb] idxvba = idxvb[ka][:,None] * nvir + idxva[kb] idxvbb = idxvb[ka][:,None] * nvir + idxvb[kb] r2_tmp = r2[kj, ka].reshape(nocc, nvir**2) r2aaa_tmp = lib.take_2d(r2_tmp, idxoa[kj], idxvaa.ravel()) r2aba_tmp = lib.take_2d(r2_tmp, idxoa[kj], idxvba.ravel()) r2bab_tmp = lib.take_2d(r2_tmp, idxob[kj], idxvab.ravel()) r2bbb_tmp = lib.take_2d(r2_tmp, idxob[kj], idxvbb.ravel()) r2aaa[kj, ka] = r2aaa_tmp.reshape(nocc_a, nvir_a, nvir_a) r2aba[kj, ka] = r2aba_tmp.reshape(nocc_a, nvir_b, nvir_a) r2bab[kj, ka] = r2bab_tmp.reshape(nocc_b, nvir_a, nvir_b) r2bbb[kj, ka] = r2bbb_tmp.reshape(nocc_b, nvir_b, nvir_b) return [r1a, r1b], [r2aaa, r2aba, r2bab, r2bbb]
[docs] def spatial2spin_ea_doublet(r1, r2, kconserv, kshift, orbspin=None): '''Convert R1/R2 of spatial orbital representation to R1/R2 of spin orbital representation''' r1a, r1b = r1 r2aaa, r2aba, r2bab, r2bbb = r2 nkpts, nocc_a, nvir_a = np.array(r2aaa.shape)[[0, 2, 3]] nkpts, nocc_b, nvir_b = np.array(r2bbb.shape)[[0, 2, 3]] if orbspin is None: orbspin = np.zeros((nocc_a+nvir_a)*2, dtype=int) orbspin[1::2] = 1 nocc = nocc_a + nocc_b nvir = nvir_a + nvir_b idxoa = [np.where(orbspin[k][:nocc] == 0)[0] for k in range(nkpts)] idxob = [np.where(orbspin[k][:nocc] == 1)[0] for k in range(nkpts)] idxva = [np.where(orbspin[k][nocc:] == 0)[0] for k in range(nkpts)] idxvb = [np.where(orbspin[k][nocc:] == 1)[0] for k in range(nkpts)] r1 = np.zeros((nvir), dtype=r1a.dtype) r1[idxva[kshift]] = r1a r1[idxvb[kshift]] = r1b r2 = np.zeros((nkpts,nkpts,nocc,nvir**2), dtype=r2aaa.dtype) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift, ka, kj] idxvaa = idxva[ka][:,None] * nvir + idxva[kb] idxvab = idxva[ka][:,None] * nvir + idxvb[kb] idxvba = idxvb[ka][:,None] * nvir + idxva[kb] idxvbb = idxvb[ka][:,None] * nvir + idxvb[kb] r2aaa_tmp = r2aaa[kj,ka].reshape(nocc_a, nvir_a*nvir_a) r2aba_tmp = r2aba[kj,ka].reshape(nocc_a, nvir_b*nvir_a) r2bab_tmp = r2bab[kj,ka].reshape(nocc_b, nvir_a*nvir_b) r2bbb_tmp = r2bbb[kj,ka].reshape(nocc_b, nvir_b*nvir_b) lib.takebak_2d(r2[kj,ka], r2aaa_tmp, idxoa[kj], idxvaa.ravel()) lib.takebak_2d(r2[kj,ka], r2aba_tmp, idxoa[kj], idxvba.ravel()) lib.takebak_2d(r2[kj,ka], r2bab_tmp, idxob[kj], idxvab.ravel()) lib.takebak_2d(r2[kj,ka], r2bbb_tmp, idxob[kj], idxvbb.ravel()) r2aab_tmp = -r2aba[kj,kb].reshape(nocc_a, nvir_b*nvir_a) r2bba_tmp = -r2bab[kj,kb].reshape(nocc_b, nvir_a*nvir_b) lib.takebak_2d(r2[kj,ka], r2bba_tmp, idxob[kj], idxvba.T.ravel()) lib.takebak_2d(r2[kj,ka], r2aab_tmp, idxoa[kj], idxvab.T.ravel()) r2 = r2.reshape(nkpts, nkpts, nocc, nvir, nvir) return r1, r2
[docs] def amplitudes_to_vector_ea(r1, r2, kshift, kconserv): nkpts, nocc, nvir = np.asarray(r2.shape)[[0,2,3]] r2_tril = np.zeros((nocc*nkpts*nvir*(nkpts*nvir-1)//2), dtype=r2.dtype) index = 0 for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] if ka < kb: idx, idy = np.tril_indices(nvir, 0) else: idx, idy = np.tril_indices(nvir, -1) r2_tril[index:index + nocc*len(idy)] = r2[kj,ka,:,idx,idy].reshape(-1) index = index + nocc*len(idy) vector = np.hstack((r1, r2_tril)) return vector
[docs] def vector_to_amplitudes_ea(vector, kshift, nkpts, nmo, nocc, kconserv): nvir = nmo - nocc r1 = vector[:nvir].copy() r2_tril = vector[nvir:].copy().reshape(nocc*nkpts*nvir*(nkpts*nvir-1)//2) r2 = np.zeros((nkpts,nkpts,nocc,nvir,nvir), dtype=vector.dtype) index = 0 for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] if ka < kb: idx, idy = np.tril_indices(nvir, 0) else: idx, idy = np.tril_indices(nvir, -1) tmp = r2_tril[index:index + nocc*len(idy)].reshape(-1,nocc) r2[kj,ka,:,idx,idy] = tmp r2[kj,kb,:,idy,idx] = -tmp index = index + nocc*len(idy) return [r1,r2]
[docs] def eaccsd(eom, nroots=1, koopmans=False, guess=None, left=False, eris=None, imds=None, partition=None, kptlist=None, dtype=None): '''See `ipccsd()` for a description of arguments.''' return ipccsd(eom, nroots, koopmans, guess, left, eris, imds, partition, kptlist, dtype)
[docs] def eaccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''2hp operators are of the form s_{ j}^{ab}, i.e. 'jb' indices are coupled.''' if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nkpts = eom.nkpts kconserv = imds.kconserv r1, r2 = vector_to_amplitudes_ea(vector, kshift, nkpts, nmo, nocc, kconserv) Hr1 = np.einsum('ac,c->a', imds.Fvv[kshift], r1) for kl in range(nkpts): Hr1 += np.einsum('ld,lad->a', imds.Fov[kl], r2[kl, kshift]) for kc in range(nkpts): Hr1 += 0.5*np.einsum('alcd,lcd->a', imds.Wvovv[kshift,kl,kc], r2[kl,kc]) Hr2 = np.zeros_like(r2) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] += np.einsum('abcj,c->jab', imds.Wvvvo[ka,kb,kshift], r1) Hr2[kj,ka] += lib.einsum('ac,jcb->jab', imds.Fvv[ka], r2[kj,ka]) Hr2[kj,ka] -= lib.einsum('bc,jca->jab', imds.Fvv[kb], r2[kj,kb]) Hr2[kj,ka] -= lib.einsum('lj,lab->jab', imds.Foo[kj], r2[kj,ka]) for kd in range(nkpts): kl = kconserv[kj, kb, kd] Hr2[kj, ka] += lib.einsum('lbdj,lad->jab', imds.Wovvo[kl, kb, kd], r2[kl, ka]) # P(ab) kl = kconserv[kj, ka, kd] Hr2[kj, ka] -= lib.einsum('ladj,lbd->jab', imds.Wovvo[kl, ka, kd], r2[kl, kb]) kc = kconserv[ka, kd, kb] Hr2[kj, ka] += 0.5 * lib.einsum('abcd,jcd->jab', imds.Wvvvv[ka, kb, kc], r2[kj, kc]) tmp = lib.einsum('xyklcd,xylcd->k', imds.Woovv[kshift, :, :], r2[:, :]) # contract_{kl, kc} Hr2[:, :] -= 0.5*lib.einsum('k,xykjab->xyjab', tmp, imds.t2[kshift, :, :]) # sum_{kj, ka] vector = eom.amplitudes_to_vector(Hr1, Hr2, kshift) return vector
[docs] def leaccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''2hp operators are of the form s_{ j}^{ab}, i.e. 'jb' indices are coupled. See also `eaccsd_matvec`''' if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nkpts = eom.nkpts kconserv = imds.kconserv r1, r2 = vector_to_amplitudes_ea(vector, kshift, nkpts, nmo, nocc, kconserv) dtype = np.result_type(r1, r2) Hr1 = np.einsum('ca,c->a', imds.Fvv[kshift], r1) for kj, kb in itertools.product(range(nkpts), repeat=2): kc = kconserv[kshift,kb,kj] Hr1 += 0.5*lib.einsum('cbaj,jcb->a',imds.Wvvvo[kc,kb,kshift],r2[kj,kc]) Hr2 = np.zeros_like(r2) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] += lib.einsum('cjab,c->jab',imds.Wvovv[kshift,kj,ka],r1) Hr2[kj,kshift] += (kj==kb)*lib.einsum('jb,a->jab',imds.Fov[kj],r1) Hr2[kj,ka] -= (kj==ka)*lib.einsum('ja,b->jab',imds.Fov[kj],r1) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] tmp1 = lib.einsum('ca,jcb->jab',imds.Fvv[ka],r2[kj,ka]) tmp1T = lib.einsum('cb,jca->jab',imds.Fvv[kb],r2[kj,kb]) Hr2[kj,ka] += (tmp1 - tmp1T) Hr2[kj,ka] += -lib.einsum('jl,lab->jab',imds.Foo[kj],r2[kj,ka]) for kd in range(nkpts): km = kconserv[kj,kb,kd] tmp2 = lib.einsum('jdbm,mad->jab',imds.Wovvo[kj,kd,kb],r2[km,ka]) km = kconserv[kj,ka,kd] tmp2T = lib.einsum('jdam,mbd->jab',imds.Wovvo[kj,kd,ka],r2[km,kb]) Hr2[kj,ka] += (tmp2 - tmp2T) kc = kconserv[ka,kd,kb] Hr2[kj,ka] += 0.5*lib.einsum('cdab,jcd->jab',imds.Wvvvv[kc,kd,ka],r2[kj,kc]) tmp = np.zeros(nocc, dtype=dtype) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] tmp += lib.einsum('jab,kjab->k',r2[kj,ka],imds.t2[kshift,kj,ka]) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] += -0.5*lib.einsum('kjab,k->jab',imds.Woovv[kshift,kj,ka],tmp) vector = eom.amplitudes_to_vector(Hr1, Hr2, kshift) return vector
[docs] def eaccsd_diag(eom, kshift, imds=None): if imds is None: imds = eom.make_imds() t1 = imds.t1 nkpts, nocc, nvir = t1.shape kconserv = imds.kconserv Hr1 = np.diag(imds.Fvv[kshift]) Hr2 = np.zeros((nkpts,nkpts,nocc,nvir,nvir), dtype=t1.dtype) if eom.partition == 'mp': # This case is untested foo = eom.eris.fock[:,:nocc,:nocc] fvv = eom.eris.fock[:,nocc:,nocc:] for kj in range(nkpts): for ka in range(nkpts): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] -= foo[kj].diagonal()[:,None,None] Hr2[kj,ka] -= fvv[ka].diagonal()[None,:,None] Hr2[kj,ka] += fvv[kb].diagonal()[None,None,:] else: for kj in range(nkpts): for ka in range(nkpts): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] -= imds.Foo[kj].diagonal()[:,None,None] Hr2[kj,ka] += imds.Fvv[ka].diagonal()[None,:,None] Hr2[kj,ka] += imds.Fvv[kb].diagonal()[None,None,:] Hr2[kj,ka] += np.einsum('jbbj->jb', imds.Wovvo[kj,kb,kb])[:, None, :] Hr2[kj,ka] += np.einsum('jaaj->ja', imds.Wovvo[kj,ka,ka])[:, :, None] if ka == kconserv[ka,kb,kb]: Hr2[kj,ka] += np.einsum('abab->ab', imds.Wvvvv[ka,kb,ka])[None,:,:] Hr2[kj,ka] -= np.einsum('kjab,kjab->jab',imds.Woovv[kshift,kj,ka],imds.t2[kshift,kj,ka]) vector = amplitudes_to_vector_ea(Hr1, Hr2, kshift, kconserv) return vector
[docs] def eaccsd_star_contract(eom, eaccsd_evals, eaccsd_evecs, leaccsd_evecs, kshift, imds=None): """ Returns: e_star (list of float): The EA-CCSD* energy. Notes: The user should check to make sure the right and left eigenvalues before running the perturbative correction. Reference: Saeh, Stanton "...energy surfaces of radicals" JCP 111, 8275 (1999); DOI:10.1063/1.480171 """ assert (eom.partition is None) if imds is None: imds = eom.make_imds() t1, t2 = imds.t1, imds.t2 eris = imds.eris nkpts, nocc, nvir = t1.shape nmo = nocc + nvir dtype = np.result_type(t1, t2) kconserv = eom.kconserv mo_energy_occ = np.array([eris.mo_energy[ki][:nocc] for ki in range(nkpts)]) mo_energy_vir = np.array([eris.mo_energy[ki][nocc:] for ki in range(nkpts)]) mo_e_o = mo_energy_occ mo_e_v = mo_energy_vir eaccsd_evecs = np.array(eaccsd_evecs) leaccsd_evecs = np.array(leaccsd_evecs) e_star = [] eaccsd_evecs, leaccsd_evecs = [np.atleast_2d(x) for x in [eaccsd_evecs, leaccsd_evecs]] eaccsd_evals = np.atleast_1d(eaccsd_evals) for ea_eval, ea_evec, ea_levec in zip(eaccsd_evals, eaccsd_evecs, leaccsd_evecs): # Enforcing <L|R> = 1 l1, l2 = vector_to_amplitudes_ea(ea_levec, kshift, nkpts, nmo, nocc, kconserv) r1, r2 = vector_to_amplitudes_ea(ea_evec, kshift, nkpts, nmo, nocc, kconserv) ldotr = np.dot(l1, r1) + 0.5 * np.dot(l2.ravel(), r2.ravel()) logger.info(eom, 'Left-right amplitude overlap : %14.8e + 1j %14.8e', ldotr.real, ldotr.imag) if abs(ldotr) < 1e-7: logger.warn(eom, 'Small %s left-right amplitude overlap. Results ' 'may be inaccurate.', ldotr) l1 /= ldotr l2 /= ldotr deltaE = 0.0 + 1j*0.0 for ki, kj in itertools.product(range(nkpts), repeat=2): lijabc = np.zeros((nkpts,nkpts,nocc,nocc,nvir,nvir,nvir),dtype=dtype) rijabc = np.zeros((nkpts,nkpts,nocc,nocc,nvir,nvir,nvir),dtype=dtype) kklist = kpts_helper.get_kconserv3(eom._cc._scf.cell, eom._cc.kpts, [ki,kj,kshift,range(nkpts),range(nkpts)]) for ka, kb in itertools.product(range(nkpts), repeat=2): #TODO: can reduce size of ijabc arrays since `kc` fixed from other k-points kc = kklist[ka,kb] # lijabc update if kc == kshift and kb == kconserv[ki,ka,kj]: lijabc[ka,kb] -= lib.einsum('ijab,c->ijabc', eris.oovv[ki,kj,ka], l1) km = kconserv[kj,ka,ki] lijabc[ka,kb] -= lib.einsum('jima,mbc->ijabc', eris.ooov[kj,ki,km], l2[km,kb]) ke = kconserv[ka,ki,kb] tmp = lib.einsum('ieab,jce->ijabc', eris.ovvv[ki,ke,ka], l2[kj,kc]) ke = kconserv[ka,kj,kb] tmpT = lib.einsum('jeab,ice->ijabc', eris.ovvv[kj,ke,ka], l2[ki,kc]) lijabc[ka,kb] -= (tmp - tmpT) # rijabc update ke = kconserv[kb,kshift,kc] tmp = lib.einsum('bcef,f->bce', eris.vvvv[kb,kc,ke], r1) tmp = lib.einsum('bce,ijae->ijabc', tmp, t2[ki,kj,ka]) rijabc[ka,kb] -= tmp km = kconserv[kj,kc,kshift] tmp = lib.einsum('mcje,e->mcj', eris.ovov[km,kc,kj], r1) tmp = lib.einsum('mcj,imab->ijabc', tmp, t2[ki,km,ka]) km = kconserv[ki,kc,kshift] tmpT = lib.einsum('mcie,e->mci', eris.ovov[km,kc,ki], r1) tmpT = lib.einsum('mci,jmab->ijabc', tmpT, t2[kj,km,ka]) rijabc[ka,kb] += (tmp - tmpT) km = kconserv[kj,ka,ki] rijabc[ka,kb] += lib.einsum('jima,mcb->ijabc', eris.ooov[kj,ki,km].conj(), r2[km,kc]) ke = kconserv[ka,ki,kb] tmp = lib.einsum('ieab,jce->ijabc', eris.ovvv[ki,ke,ka].conj(), r2[kj,kc]) ke = kconserv[ka,kj,kb] tmpT = lib.einsum('jeab,ice->ijabc', eris.ovvv[kj,ke,ka].conj(), r2[ki,kc]) rijabc[ka,kb] -= (tmp - tmpT) eabc = np.zeros((nkpts,nkpts,nvir,nvir,nvir), dtype=dtype) Plijabc = np.zeros_like(lijabc) Prijabc = np.zeros_like(rijabc) for ka, kb in itertools.product(range(nkpts), repeat=2): kc = kklist[ka,kb] # P(abc) Plijabc[ka,kb] = (lijabc[ka,kb] + lijabc[kb,kc].transpose(0,1,4,2,3) + lijabc[kc,ka].transpose(0,1,3,4,2)) Prijabc[ka,kb] = (rijabc[ka,kb] + rijabc[kb,kc].transpose(0,1,4,2,3) + rijabc[kc,ka].transpose(0,1,3,4,2)) eabc[ka,kb] = _get_epqr([0,nvir,ka,mo_e_v,eom.nonzero_vpadding], [0,nvir,kb,mo_e_v,eom.nonzero_vpadding], [0,nvir,kc,mo_e_v,eom.nonzero_vpadding], fac=[-1.,]*3) # Creating denominator eij = _get_epq([0,nocc,ki,mo_e_o,eom.nonzero_opadding], [0,nocc,kj,mo_e_o,eom.nonzero_opadding]) eijabc = (eij[None, None, :, :, None, None, None] + eabc[:, :, None, None, :, :, :]) denom = eijabc + ea_eval denom = 1. / denom deltaE += lib.einsum('xyijabc,xyijabc,xyijabc', Plijabc, Prijabc, denom) deltaE *= 1./12 deltaE = deltaE.real logger.info(eom, "Exc. energy, delta energy = %16.12f, %16.12f", ea_eval + deltaE, deltaE) e_star.append(ea_eval + deltaE) return e_star
[docs] def eaccsd_star(eom, nroots=1, koopmans=False, right_guess=None, left_guess=None, eris=None, imds=None, partition=None, kptlist=None, dtype=None, **kwargs): '''See `kernel()` for a description of arguments.''' if partition: raise NotImplementedError return perturbed_ccsd_kernel(eom, nroots=nroots, koopmans=koopmans, right_guess=right_guess, left_guess=left_guess, eris=eris, imds=imds, partition=partition, kptlist=kptlist, dtype=dtype)
[docs] def mask_frozen_ea(eom, vector, kshift, const=LARGE_DENOM): '''Replaces all frozen orbital indices of `vector` with the value `const`.''' r1, r2 = eom.vector_to_amplitudes(vector, kshift=kshift) kconserv = eom.kconserv nkpts = eom.nkpts # Get location of padded elements in occupied and virtual space nonzero_opadding, nonzero_vpadding = eom.nonzero_opadding, eom.nonzero_vpadding new_r1 = const * np.ones_like(r1) new_r2 = const * np.ones_like(r2) new_r1[nonzero_vpadding[kshift]] = r1[nonzero_vpadding[kshift]] for kj in range(nkpts): for ka in range(nkpts): kb = kconserv[kshift, ka, kj] idx = np.ix_([kj], [ka], nonzero_opadding[kj], nonzero_vpadding[ka], nonzero_vpadding[kb]) new_r2[idx] = r2[idx] return eom.amplitudes_to_vector(new_r1, new_r2, kshift, kconserv)
[docs] class EOMEA(eom_rccsd.EOMEA): def __init__(self, cc): self.kpts = cc.kpts self.nonzero_opadding, self.nonzero_vpadding = self.get_padding_k_idx(cc) self.kconserv = cc.khelper.kconserv eom_rccsd.EOM.__init__(self, cc) kernel = eaccsd eaccsd = eaccsd eaccsd_star = eaccsd_star ccsd_star_contract = eaccsd_star_contract get_diag = eaccsd_diag matvec = eaccsd_matvec l_matvec = leaccsd_matvec mask_frozen = mask_frozen_ea get_padding_k_idx = get_padding_k_idx
[docs] def eaccsd_star_contract(self, eaccsd_evals, eaccsd_evecs, leaccsd_evecs, kshift, imds=None): return self.ccsd_star_contract(eaccsd_evals, eaccsd_evecs, leaccsd_evecs, kshift, imds=imds)
[docs] def get_init_guess(self, kshift, nroots=1, koopmans=False, diag=None): size = self.vector_size() dtype = getattr(diag, 'dtype', np.complex128) nroots = min(nroots, size) guess = [] if koopmans: for n in self.nonzero_vpadding[kshift][:nroots]: g = np.zeros(int(size), dtype=dtype) g[n] = 1.0 g = self.mask_frozen(g, kshift, const=0.0) guess.append(g) else: idx = diag.argsort()[:nroots] for i in idx: g = np.zeros(int(size), dtype=dtype) g[i] = 1.0 g = self.mask_frozen(g, kshift, const=0.0) guess.append(g) return guess
@property def nkpts(self): return len(self.kpts)
[docs] def gen_matvec(self, kshift, imds=None, left=False, **kwargs): if imds is None: imds = self.make_imds() diag = self.get_diag(kshift, imds) if left: matvec = lambda xs: [self.l_matvec(x, kshift, imds, diag) for x in xs] else: matvec = lambda xs: [self.matvec(x, kshift, imds, diag) for x in xs] return matvec, diag
[docs] def vector_to_amplitudes(self, vector, kshift=None, nkpts=None, nmo=None, nocc=None, kconserv=None): if nmo is None: nmo = self.nmo if nocc is None: nocc = self.nocc if nkpts is None: nkpts = self.nkpts if kconserv is None: kconserv = self.kconserv return vector_to_amplitudes_ea(vector, kshift, nkpts, nmo, nocc, kconserv)
[docs] def amplitudes_to_vector(self, r1, r2, kshift, kconserv=None): if kconserv is None: kconserv = self.kconserv return amplitudes_to_vector_ea(r1, r2, kshift, kconserv)
[docs] def vector_size(self): nocc = self.nocc nvir = self.nmo - nocc nkpts = self.nkpts return nvir + nocc*nkpts*nvir*(nkpts*nvir-1)//2
[docs] def make_imds(self, eris=None): imds = _IMDS(self._cc, eris) imds.make_ea() return imds
[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) return imds
######################################## # EOM-EE-CCSD ########################################
[docs] def kernel_ee(eom, nroots=1, koopmans=False, guess=None, left=False, eris=None, imds=None, partition=None, kptlist=None, dtype=None, **kwargs): '''See `kernel()` for a description of arguments. This method is merely a simplified version of kernel() with a few parts removed, such as those involving `eom.mask_frozen()`. Slowly they will be added back for the completion of program. ''' 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=eris) nkpts = eom.nkpts if kptlist is None: kptlist = range(nkpts) # TODO mask frozen-orbital indices if dtype is None: dtype = np.result_type(*imds.t1) # Note that vector_size may change with kshift. Thus we do not fix # the length of each eval and evec evals = [None]*len(kptlist) evecs = [None]*len(kptlist) convs = [None]*len(kptlist) for k, kshift in enumerate(kptlist): print("\nkshift =", kshift) # vector size and thus, nroots depend on kshift in the case of even nkpts, size = eom.vector_size(kshift) nroots = min(nroots, size) matvec, diag = eom.gen_matvec(kshift, imds, left=left, **kwargs) if diag.size != size: raise ValueError("Number of diagonal elements in effective H does not match R vector size") # TODO update `diag` in case of frozen orbitals # TODO allow user provided guess vector # Since vector_size may change with kshift, it is difficult for users to # provide guesses. Similarly, `guess` from the previous `kshift` may not # work for the current `kshift` due to different vector_size. Thus for # now we keep `user_guess` false, and always compute `guess` on our own. user_guess = False guess = eom.get_init_guess(kshift, nroots, koopmans=koopmans, diag=diag, imds=imds) for ig, g in enumerate(guess): guess_norm = np.linalg.norm(g) guess_norm_tol = LOOSE_ZERO_TOL if guess_norm < guess_norm_tol: raise ValueError('Guess vector (id=%d) with norm %.4g is below threshold %.4g.\n' 'This could possibly be due to masking/freezing orbitals.\n' 'Check your guess vector to make sure it has sufficiently large norm.' % (ig, guess_norm, guess_norm_tol)) def precond(r, e0, x0): return r/(e0-diag+1e-12) eig = lib.davidson_nosym1 # TODO allow user provided guess vector or Koopmans if user_guess or koopmans: def pickeig(w, v, nr, envs): x0 = lib.linalg_helper._gen_x0(envs['v'], envs['xs']) idx = np.argmax( np.abs(np.dot(np.array(guess).conj(),np.array(x0).T)), axis=1 ) return lib.linalg_helper._eigs_cmplx2real(w, v, idx) conv_k, evals_k, evecs_k = eig(matvec, guess, precond, pick=pickeig, tol=eom.conv_tol, max_cycle=eom.max_cycle, max_space=eom.max_space, max_memory=eom.max_memory, nroots=nroots, verbose=eom.verbose) else: conv_k, evals_k, evecs_k = eig(matvec, guess, precond, tol=eom.conv_tol, max_cycle=eom.max_cycle, max_space=eom.max_space, max_memory=eom.max_memory, nroots=nroots, verbose=eom.verbose) evals_k = evals_k.real evals[k] = evals_k evecs[k] = evecs_k convs[k] = conv_k for n, en, vn in zip(range(nroots), evals_k, evecs_k): r1, r2 = eom.vector_to_amplitudes(vn, kshift=kshift) 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 = %0.6g', n, en, qp_weight) log.timer('EOM-CCSD', *cput0) return convs, evals, evecs
[docs] def eeccsd(eom, nroots=1, koopmans=False, guess=None, left=False, eris=None, imds=None, partition=None, kptlist=None, dtype=None): '''See `kernel_ee()` for a description of arguments.''' eom.converged, eom.e, eom.v \ = kernel_ee(eom, nroots, koopmans, guess, left, eris=eris, imds=imds, partition=partition, kptlist=kptlist, dtype=dtype) return eom.e, eom.v
[docs] def eeccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''Spin-orbital EOM-EE-CCSD equations with k points.''' # Ref: Wang, Tu, and Wang, J. Chem. Theory Comput. 10, 5567 (2014) Eqs.(9)-(10) # Note: Last line in Eq. (10) is superfluous. # See, e.g. Gwaltney, Nooijen, and Barlett, Chem. Phys. Lett. 248, 189 (1996) if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nvir = nmo - nocc nkpts = eom.nkpts kconserv = imds.kconserv kconserv_r1 = eom.get_kconserv_ee_r1(kshift) kconserv_r2 = eom.get_kconserv_ee_r2(kshift) r1, r2 = vector_to_amplitudes_ee(vector, kshift, nkpts, nmo, nocc, kconserv_r2) Hr1 = np.zeros_like(r1) for ki in range(nkpts): ka = kconserv_r1[ki] Hr1[ki] += np.einsum('ae,ie->ia', imds.Fvv[ka], r1[ki]) Hr1[ki] -= np.einsum('mi,ma->ia', imds.Foo[ki], r1[ki]) for km in range(nkpts): Hr1[ki] += np.einsum('me,imae->ia', imds.Fov[km], r2[ki, km, ka]) ke = kconserv_r1[km] Hr1[ki] += np.einsum('maei,me->ia', imds.Wovvo[km, ka, ke], r1[km]) for kn in range(nkpts): Hr1[ki] -= 0.5*np.einsum('mnie,mnae->ia', imds.Wooov[km, kn, ki], r2[km, kn, ka]) # Rename dummy index kn->ke Hr1[ki] += 0.5*np.einsum('amef,imef->ia', imds.Wvovv[ka, km, kn], r2[ki, km, kn]) Hr2 = np.zeros_like(r2) for ki, kj, ka in kpts_helper.loop_kkk(nkpts): kb = kconserv_r2[ki, ka, kj] # r_ijab <- P(ij) (-F_mj r_imab) # km - kj = G # => km = kj tmp_ij = np.einsum('mj,imab->ijab', -imds.Foo[kj], r2[ki, kj, ka]) # r_ijab <- P(ij) W_abej r_ie ke = kconserv_r1[ki] tmp_ij += np.einsum('abej,ie->ijab', imds.Wvvvo[ka, kb, ke], r1[ki]) Hr2[ki, kj, ka] += tmp_ij Hr2[kj, ki, ka] -= tmp_ij.transpose(1, 0, 2, 3) # r_ijab <- P(ab) F_be r_ijae tmp_ab = np.einsum('be,ijae->ijab', imds.Fvv[kb], r2[ki, kj, ka]) # r_ijab <- P(ab) (- W_mbij r_ma) # km + kb - ki - kj = G # => ki + kj - km - kb = G km = kconserv[ki, kb, kj] tmp_ab += np.einsum('mbij,ma->ijab', -imds.Wovoo[km, kb, ki], r1[km]) Hr2[ki, kj, ka] += tmp_ab Hr2[ki, kj, kb] -= tmp_ab.transpose(0, 1, 3, 2) # r_ijab <- 0.5 W_mnij r_mnab tmpoooo = np.zeros((nocc, nocc, nvir, nvir), dtype=r2.dtype) # r_ijab <- 0.5 W_abef r_ijef tmpvvvv = np.zeros((nocc, nocc, nvir, nvir), dtype=r2.dtype) for km in range(nkpts): # km + kn - ki - kj = G (as in W_mnij) kn = kconserv[ki, km, kj] tmpoooo += 0.5*np.einsum('mnij,mnab->ijab', imds.Woooo[km, kn, ki], r2[km, kn, ka]) # Rename dummy index km->ke tmpvvvv += 0.5*np.einsum('abef,ijef->ijab', imds.Wvvvv[ka, kb, km], r2[ki, kj, km]) Hr2[ki, kj, ka] += tmpoooo Hr2[ki, kj, ka] += tmpvvvv # r_ijab <- P(ij) P(ab) W_mbej r_imae for km in range(nkpts): # km + kb - ke - kj = G ke = kconserv[km, kj, kb] tmp = np.einsum('mbej,imae->ijab', imds.Wovvo[km, kb, ke], r2[ki, km, ka]) Hr2[ki, kj, ka] += tmp Hr2[kj, ki, ka] -= tmp.transpose(1, 0, 2, 3) Hr2[ki, kj, kb] -= tmp.transpose(0, 1, 3, 2) Hr2[kj, ki, kb] += tmp.transpose(1, 0, 3, 2) # # r_ijab <- P(ab) (-0.5 W_mnef t_ijae r_mnbf) # r_ijab <- P(ab) W_amfe t_ijfb r_me # r_ijab <- P(ij) (-0.5 W_mnef t_imab r_jnef) # r_ijab <- P(ij) W_mnie t_njab r_me # # Build intermediates M = W.r2 for the four terms above tmp_eb = np.zeros((nkpts, nvir, nvir), dtype=r2.dtype) tmp_fa = np.zeros_like(tmp_eb) tmp_jm = np.zeros((nkpts, nocc, nocc), dtype=r2.dtype) tmp_in = np.zeros_like(tmp_jm) for ke in range(nkpts): # M_eb = W_mnef r_mnbf (or equivalently, M_ea = W_mnef r_mnaf) # km + kn - ke - kf = G # km + kn - kb - kf = G + kshift # => ke - kb = G + kshift kb = kconserv_r1[ke] # x: km, y: kn tmp_eb[ke] += np.einsum('xymnef,xymnbf->eb', imds.Woovv[:, :, ke], r2[:, :, kb]) # M_fa = W_amfe r_me (or equivalently, M_fb = W_bmfe r_me) kf = ke # ki + kj - ka - kb = G + kshift # ki + kj - kf - kb = G # => kf - ka = G + kshift ka = kconserv_r1[kf] # x: km tmp_fa[kf] += np.einsum('xamfe,xme->fa', imds.Wvovv[ka, :, kf], r1) # M_jm = W_mnef r_jnef (or equivalently, M_im = W_mnef r_inef) kj = ke # km + kn - ke - kf = G # kj + kn - ke - kf = G + kshift # => kj - km = G + kshift km = kconserv_r1[kj] # x: kn, y: ke tmp_jm[kj] += np.einsum('xymnef,xyjnef->jm', imds.Woovv[km], r2[kj]) # M_in = W_mnie r_me (or equivalently, M_jn = W_mnje r_me) ki = ke # ki + kj - ka - kb = G + kshift # kn + kj - ka - kb = G # => ki - kn = G + kshift kn = kconserv_r1[ki] # x: km tmp_in[ki] += np.einsum('xmnie,xme->in', imds.Wooov[:, kn, ki], r1) for ki, kj, ka in kpts_helper.loop_kkk(nkpts): kb = kconserv_r2[ki, ka, kj] # r_ijab <- P(ab) (-0.5 M_eb t_ijae) # ki + kj - ka - ke = G ke = kconserv[ki, ka, kj] tmp_ab = 0.5*np.einsum('eb,ijae->ijab', -tmp_eb[ke], imds.t2[ki, kj, ka]) # r_ijab <- P(ab) M_fa t_ijfb # ki + kj - kf - kb = G kf = kconserv[ki, kb, kj] tmp_ab += np.einsum('fa,ijfb->ijab', tmp_fa[kf], imds.t2[ki, kj, kf]) Hr2[ki, kj, ka] += tmp_ab Hr2[ki, kj, kb] -= tmp_ab.transpose(0, 1, 3, 2) # r_ijab <- P(ij) (-0.5 M_jm t_imab) # kj - km = G + kshift km = kconserv_r1[kj] tmp_ij = 0.5*np.einsum('jm,imab->ijab', -tmp_jm[kj], imds.t2[ki, km, ka]) # r_ijab <- P(ij) M_in t_njab # ki - kn = G + kshift kn = kconserv_r1[ki] tmp_ij += np.einsum('in,njab->ijab', tmp_in[ki], imds.t2[kn, kj, ka]) Hr2[ki, kj, ka] += tmp_ij Hr2[kj, ki, ka] -= tmp_ij.transpose(1, 0, 2, 3) vector = amplitudes_to_vector_ee(Hr1, Hr2, kshift, kconserv_r2) return vector
[docs] def eeccsd_diag(eom, kshift, imds=None): '''Diagonal elements of similarity-transformed Hamiltonian''' if imds is None: imds = eom.make_imds() t1 = imds.t1 nkpts, nocc, nvir = t1.shape kconserv = eom.kconserv kconserv_r1 = eom.get_kconserv_ee_r1(kshift) kconserv_r2 = eom.get_kconserv_ee_r2(kshift) Hr1 = np.zeros((nkpts, nocc, nvir), dtype=t1.dtype) for ki in range(nkpts): ka = kconserv_r1[ki] Hr1[ki] -= imds.Foo[ki].diagonal()[:, None] Hr1[ki] += imds.Fvv[ka].diagonal()[None, :] Hr1[ki] += np.einsum('iaai->ia', imds.Wovvo[ki, ka, ka]) Hr2 = np.zeros((nkpts, nkpts, nkpts, nocc, nocc, nvir, nvir), dtype=t1.dtype) # TODO allow partition='mp' if eom.partition == "mp": raise NotImplementedError else: for ki, kj, ka in kpts_helper.loop_kkk(nkpts): kb = kconserv_r2[ki, ka, kj] Hr2[ki, kj, ka] -= imds.Foo[ki].diagonal()[:, None, None, None] Hr2[ki, kj, ka] -= imds.Foo[kj].diagonal()[None, :, None, None] Hr2[ki, kj, ka] += imds.Fvv[ka].diagonal()[None, None, :, None] Hr2[ki, kj, ka] += imds.Fvv[kb].diagonal()[None, None, None, :] Hr2[ki, kj, ka] += np.einsum('jbbj->jb', imds.Wovvo[kj, kb, kb])[None, :, None, :] Hr2[ki, kj, ka] += np.einsum('ibbi->ib', imds.Wovvo[ki, kb, kb])[:, None, None, :] Hr2[ki, kj, ka] += np.einsum('jaaj->ja', imds.Wovvo[kj, ka, ka])[None, :, :, None] Hr2[ki, kj, ka] += np.einsum('iaai->ia', imds.Wovvo[ki, ka, ka])[:, None, :, None] Hr2[ki, kj, ka] += np.einsum('ijij->ij', imds.Woooo[ki, kj, ki])[:, :, None, None] Hr2[ki, kj, ka] += np.einsum('abab->ab', imds.Wvvvv[ka, kb, ka])[None, None, :, :] # This is to make t2 are non-zero # Note that `kconserv` is used instead of `kconserv_r2` kk = kconserv[ka, kj, kb] Hr2[ki, kj, ka] -= np.einsum('kjab,kjab->jab', imds.Woovv[kk, kj, ka], imds.t2[kk, kj, ka])[None, :, :, :] kk = kconserv[ka, ki, kb] Hr2[ki, kj, ka] -= np.einsum('kiab,kiab->iab', imds.Woovv[kk, ki, ka], imds.t2[kk, ka, ka])[:, None, :, :] kc = kconserv[ki, kb, kj] Hr2[ki, kj, ka] -= np.einsum('ijcb,ijcb->ijb', imds.Woovv[ki, kj, kc], imds.t2[ki, kj, kc])[:, :, None, :] kc = kconserv[ki, ka, kj] Hr2[ki, kj, ka] -= np.einsum('ijca,ijca->ija', imds.Woovv[ki, kj, kc], imds.t2[ki, kj, kc])[:, :, :, None] # Make sure 4th argument you pass is `kconserv_r2` vector = amplitudes_to_vector_ee(Hr1, Hr2, kshift, kconserv_r2) return vector
[docs] def vector_to_amplitudes_ee(vector, kshift, nkpts, nmo, nocc, kconserv): '''Transform 1-dimensional array to 3- and 7-dimensional arrays, r1 and r2. For example: vector: a 1-d array with all r1 elements, and r2 elements whose indices satisfy (i k_i) > (j k_j) and (a k_a) > (b k_b) return: [r1, r2], where r1 = r_{i k_i}^{a k_a} is a 3-d array whose elements can be accessed via r1[k_i, i, a]. r2 = r_{i k_i, j k_j}^{a k_a, b k_b} is a 7-d array whose elements can be accessed via r2[k_i, k_j, k_a, i, j, a, b] ''' nvir = nmo - nocc r1 = vector[:nkpts*nocc*nvir].copy().reshape(nkpts, nocc, nvir) ki_i, kj_j = np.tril_indices(nkpts*nocc, -1) ida, idb = np.tril_indices(nvir, -1) r2 = np.zeros((nkpts*nocc, nkpts*nocc, nkpts, nvir, nvir), dtype=vector.dtype) r2_tril = vector[nkpts*nocc*nvir:].copy() offset = 0 nvir2_tril = nvir*(nvir-1)//2 nvir2 = nvir*nvir for ij in range(len(ki_i)): idx_ki_i = ki_i[ij] idx_kj_j = kj_j[ij] ki = idx_ki_i // nocc kj = idx_kj_j // nocc r2_ka_ab = np.zeros((nkpts, nvir, nvir), dtype=r2_tril.dtype) for ka in range(nkpts): kb = kconserv[ki, ka, kj] if ka == kb: tmp = r2_tril[offset:offset+nvir2_tril] r2_ka_ab[ka, ida, idb] = tmp r2_ka_ab[ka, idb, ida] = -tmp offset += nvir2_tril elif ka > kb: tmp = r2_tril[offset:offset+nvir2].reshape(nvir, nvir) r2_ka_ab[ka] = tmp r2_ka_ab[kb] = -tmp.transpose() offset += nvir2 r2[idx_ki_i, idx_kj_j] = r2_ka_ab r2[idx_kj_j, idx_ki_i] = -r2_ka_ab r2 = r2.reshape(nkpts, nocc, nkpts, nocc, nkpts, nvir, nvir).transpose(0, 2, 4, 1, 3, 5, 6) return [r1, r2]
[docs] def amplitudes_to_vector_ee(r1, r2, kshift, kconserv): '''Transform 3- and 7-dimensional arrays, r1 and r2, to a 1-dimensional array with unique indices. For example: r1: t_{i k_i}^{a k_a} r2: t_{i k_i, j k_j}^{a k_a, b k_b} return: a vector with all r1 elements, and r2 elements whose indices satisfy (i k_i) > (j k_j) and (a k_a) > (b k_b) ''' # r1 indices: k_i, i, a nkpts, nocc, nvir = np.asarray(r1.shape)[[0, 1, 2]] # r2 indices (old): k_i, k_j, k_a, i, j, a, b # r2 indices (new): (k_i, i), (k_j, j), (k_a, a, b) r2 = r2.transpose(0, 3, 1, 4, 2, 5, 6).reshape(nkpts*nocc, nkpts*nocc, nkpts, nvir, nvir) # Get (k_i, i) and (k_j, j) indices for the lower-triangle of r2 ki_i, kj_j = np.tril_indices(nkpts*nocc, -1) ida, idb = np.tril_indices(nvir, -1) vector = r1.ravel() for ij in range(len(ki_i)): ki = ki_i[ij] // nocc kj = kj_j[ij] // nocc r2ab = r2[ki_i[ij], kj_j[ij]] for ka in range(nkpts): kb = kconserv[ki, ka, kj] if ka == kb: vector = np.hstack((vector, r2ab[ka, ida, idb])) elif ka > kb: vector = np.hstack((vector, r2ab[ka].ravel())) return vector
[docs] class EOMEE(eom_rccsd.EOM): def __init__(self, cc): self.kpts = cc.kpts self.kconserv = cc.khelper.kconserv # debug self.debug_vals = np.array([None]*len(cc.kpts)) eom_rccsd.EOM.__init__(self, cc) kernel = eeccsd eeccsd = eeccsd matvec = eeccsd_matvec get_diag = eeccsd_diag @property def nkpts(self): return len(self.kpts)
[docs] def vector_size(self, kshift=0): '''Size of the linear excitation operator R vector based on spin-orbital basis. Kwargs: kshift : int index of kpt in R(k) Returns: size (int): number of unique elements in linear excitation operator R Notes: The vector size is kshift-dependent if nkpts is an even number ''' nocc = self.nocc # alpha+beta nvir = self.nmo-nocc # alpha+beta nkpts = self.nkpts size_r1 = nkpts*nocc*nvir if nkpts % 2 == 1: size_r2 = nkpts*nocc*(nkpts*nocc-1)//2*nvir*(nkpts*nvir-1)//2 else: size_oo = nocc*(nocc-1)//2 # When ki==kj, there are size_oo ways to create 2 holes size_vv = nvir*(nvir-1)//2 # When ka==kb, there are size_vv ways to create 2 particles size_r2 = 0 kconserv = self.get_kconserv_ee_r2(kshift) # TODO Optimize this 3-layer for loop, or find an elegant solution for ki, kj, ka in kpts_helper.loop_kkk(nkpts): kb = kconserv[ki, ka, kj] if ki == kj: if ka == kb: size_r2 += size_oo*size_vv elif ka > kb: size_r2 += size_oo*nvir**2 elif ki > kj: if ka == kb: size_r2 += nocc**2*size_vv elif ka > kb: size_r2 += nocc**2*nvir**2 return size_r1 + size_r2
[docs] def get_init_guess(self, kshift, nroots=1, koopmans=False, diag=None, **kwargs): """Initial guess vectors of R coefficients""" size = self.vector_size(kshift) dtype = getattr(diag, 'dtype', np.complex128) nroots = min(nroots, size) guess = [] # TODO do Koopmans later if koopmans: raise NotImplementedError else: idx = diag.argsort()[:nroots] for i in idx: g = np.zeros(int(size), dtype=dtype) g[i] = 1.0 # TODO do mask_frozen later guess.append(g) return guess
[docs] def gen_matvec(self, kshift, imds=None, left=False, **kwargs): if imds is None: imds = self.make_imds() diag = self.get_diag(kshift, imds) if left: # TODO allow left vectors to be computed raise NotImplementedError else: matvec = lambda xs: [self.matvec(x, kshift, imds, diag) for x in xs] return matvec, diag
[docs] def vector_to_amplitudes(self, vector, kshift=None, nkpts=None, nmo=None, nocc=None, kconserv=None): if nmo is None: nmo = self.nmo if nocc is None: nocc = self.nocc if nkpts is None: nkpts = self.nkpts if kconserv is None: kconserv = self.get_kconserv_ee_r2(kshift) return vector_to_amplitudes_ee(vector, kshift, nkpts, nmo, nocc, kconserv)
[docs] def amplitudes_to_vector(self, r1, r2, kshift, kconserv=None): if kconserv is None: kconserv = self.get_kconserv_ee_r2(kshift) return amplitudes_to_vector_ee(r1, r2, kshift, kconserv)
[docs] def get_kconserv_ee_r1(self, kshift=0): r'''Get the momentum conservation array for a set of k-points. Given k-point index m the array kconserv_r1[m] returns the index n that satisfies momentum conservation, (k(m) - k(n) - kshift) \dot a = 2n\pi This is used for symmetry of 1p-1h excitation operator vector R_{m k_m}^{n k_n} is zero unless n satisfies the above. Note that this method is adapted from `kpts_helper.get_kconserv()`. ''' kconserv_r1 = self.kconserv[:,kshift,0].copy() return kconserv_r1
# TODO merge it with `kpts_helper.get_kconserv()`
[docs] def get_kconserv_ee_r2(self, kshift=0): r'''Get the momentum conservation array for a set of k-points. Given k-point indices (k, l, m) the array kconserv_r2[k,l,m] returns the index n that satisfies momentum conservation, (k(k) - k(l) + k(m) - k(n) - kshift) \dot a = 2n\pi This is used for symmetry of 2p-2h excitation operator vector R_{k k_k, m k_m}^{l k_l n k_n} is zero unless n satisfies the above. Note that this method is adapted from `kpts_helper.get_kconserv()`. ''' cell = self._cc._scf.cell kpts = self.kpts nkpts = kpts.shape[0] a = cell.lattice_vectors() / (2 * np.pi) kconserv_r2 = np.zeros((nkpts, nkpts, nkpts), dtype=int) kvKLM = kpts[:, None, None, :] - kpts[:, None, :] + kpts # Apply k shift kvKLM = kvKLM - kpts[kshift] for N, kvN in enumerate(kpts): kvKLMN = np.einsum('wx,klmx->wklm', a, kvKLM - kvN) # check whether (1/(2pi) k_{KLMN} dot a) is an integer kvKLMN_int = np.rint(kvKLMN) mask = np.einsum('wklm->klm', abs(kvKLMN - kvKLMN_int)) < 1e-9 kconserv_r2[mask] = N return kconserv_r2
[docs] def make_imds(self, eris=None): imds = _IMDS(self._cc, eris=eris) imds.make_ee() return imds
class _IMDS: # Exactly the same as RCCSD IMDS except # -- rintermediates --> gintermediates # -- Loo, Lvv, cc_Fov --> Foo, Fvv, Fov # -- One less 2-virtual intermediate def __init__(self, cc, eris=None): self._cc = cc self.verbose = cc.verbose self.kconserv = kpts_helper.get_kconserv(cc._scf.cell, cc.kpts) self.stdout = cc.stdout self.t1, self.t2 = cc.t1, cc.t2 if eris is None: eris = cc.ao2mo() self.eris = eris self._made_shared = False self.made_ip_imds = False self.made_ea_imds = False self.made_ee_imds = False def _make_shared(self): cput0 = (logger.process_clock(), logger.perf_counter()) kconserv = self.kconserv t1, t2, eris = self.t1, self.t2, self.eris self.Foo = imd.Foo(self._cc, t1, t2, eris, kconserv) self.Fvv = imd.Fvv(self._cc, t1, t2, eris, kconserv) self.Fov = imd.Fov(self._cc, t1, t2, eris, kconserv) # 2 virtuals self.Wovvo = imd.Wovvo(self._cc, t1, t2, eris, kconserv) self.Woovv = eris.oovv self._made_shared = True logger.timer_debug1(self, 'EOM-CCSD shared intermediates', *cput0) return self def make_ip(self): if not self._made_shared: self._make_shared() cput0 = (logger.process_clock(), logger.perf_counter()) kconserv = self.kconserv t1, t2, eris = self.t1, self.t2, self.eris # 0 or 1 virtuals self.Woooo = imd.Woooo(self._cc, t1, t2, eris, kconserv) self.Wooov = imd.Wooov(self._cc, t1, t2, eris, kconserv) self.Wovoo = imd.Wovoo(self._cc, t1, t2, eris, kconserv) self.made_ip_imds = True logger.timer_debug1(self, 'EOM-CCSD IP intermediates', *cput0) return self def make_t3p2_ip(self, cc): 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 = False # Force update self.make_ip() # Make after t1/t2 updated self.Wovoo = self.Wovoo + Wovoo self.made_ip_imds = True logger.timer_debug1(self, 'EOM-CCSD(T)a IP intermediates', *cput0) return self def make_ea(self): if not self._made_shared: self._make_shared() cput0 = (logger.process_clock(), logger.perf_counter()) kconserv = self.kconserv t1, t2, eris = self.t1, self.t2, self.eris # FIXME DELETE WOOOO # 0 or 1 virtuals self.Woooo = imd.Woooo(self._cc, t1, t2, eris, kconserv) # 3 or 4 virtuals self.Wvovv = imd.Wvovv(self._cc, t1, t2, eris, kconserv) self.Wvvvv = imd.Wvvvv(self._cc, t1, t2, eris, kconserv) self.Wvvvo = imd.Wvvvo(self._cc, t1, t2, eris, kconserv) self.made_ea_imds = True logger.timer_debug1(self, 'EOM-CCSD EA intermediates', *cput0) return self def make_t3p2_ea(self, cc): 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 = False # Force update self.make_ea() # Make after t1/t2 updated self.Wvvvo = self.Wvvvo + Wvvvo self.made_ea_imds = True logger.timer_debug1(self, 'EOM-CCSD(T)a EA intermediates', *cput0) return self def make_ee(self): if not self._made_shared: self._make_shared() cput0 = (logger.process_clock(), logger.perf_counter()) kconserv = self.kconserv t1, t2, eris = self.t1, self.t2, self.eris if not self.made_ip_imds: # 0 or 1 virtuals self.Woooo = imd.Woooo(self._cc, t1, t2, eris, kconserv) self.Wooov = imd.Wooov(self._cc, t1, t2, eris, kconserv) self.Wovoo = imd.Wovoo(self._cc, t1, t2, eris, kconserv) if not self.made_ea_imds: # 3 or 4 virtuals self.Wvovv = imd.Wvovv(self._cc, t1, t2, eris, kconserv) self.Wvvvv = imd.Wvvvv(self._cc, t1, t2, eris, kconserv) self.Wvvvo = imd.Wvvvo(self._cc, t1, t2, eris, kconserv, self.Wvvvv) self.made_ee_imds = True logger.timer(self, 'EOM-CCSD EE intermediates', *cput0) return self if __name__ == '__main__': from pyscf.pbc import gto, scf, cc cell = gto.Cell() cell.atom=''' C 0.000000000000 0.000000000000 0.000000000000 C 1.685068664391 1.685068664391 1.685068664391 ''' cell.basis = { 'C': [[0, (0.8, 1.0)], [1, (1.0, 1.0)]]} cell.pseudo = 'gth-pade' cell.a = ''' 0.000000000, 3.370137329, 3.370137329 3.370137329, 0.000000000, 3.370137329 3.370137329, 3.370137329, 0.000000000''' cell.unit = 'B' cell.verbose = 5 cell.build() # Running HF and CCSD with 1x1x2 Monkhorst-Pack k-point mesh kmf = scf.KRHF(cell, kpts=cell.make_kpts([1,1,2]), exxdiv=None) kmf.conv_tol_grad = 1e-8 ehf = kmf.kernel() mycc = cc.KGCCSD(kmf) mycc.conv_tol = 1e-12 mycc.conv_tol_normt = 1e-10 eris = mycc.ao2mo(mycc.mo_coeff) ecc, t1, t2 = mycc.kernel() print(ecc - -0.155298393321855) eom = EOMIP(mycc) e, v = eom.ipccsd(nroots=2, kptlist=[0]) eom = EOMEA(mycc) eom.max_cycle = 100 e, v = eom.eaccsd(nroots=2, koopmans=True, kptlist=[0])