Source code for pyscf.pbc.gw.kgw_slow_supercell

"""
This module implements the G0W0 approximation on top of `pyscf.tdscf.rhf_slow` TDHF implementation. Unlike `gw.py`, all
integrals are stored in memory. Several variants of GW are available:

 * `pyscf.gw_slow`: the molecular implementation;
 * `pyscf.pbc.gw.gw_slow`: single-kpoint PBC (periodic boundary condition) implementation;
 * (this module) `pyscf.pbc.gw.kgw_slow_supercell`: a supercell approach to PBC implementation with multiple k-points.
   Runs the molecular code for a model with several k-points for the cost of discarding momentum conservation and using
   dense instead of sparse matrixes;
 * `pyscf.pbc.gw.kgw_slow`: a PBC implementation with multiple k-points;
"""

from pyscf.gw import gw_slow
from pyscf.lib import einsum

import numpy


# Convention for these modules:
# * IMDS contains routines for intermediates
# * kernel finds GW roots
# * GW provides a container


[docs] def corrected_moe(eri, k, p): """ Calculates the corrected orbital energy. Args: eri (PhysERI): a container with electron repulsion integrals; k (int): the k-point index; p (int): orbital; Returns: The corrected orbital energy. """ moe = eri.mo_energy[k][p] moc = eri.mo_coeff[k][:, p] nk = len(eri.mo_energy) vk = 0 for k2 in range(nk): vk -= eri.ao2mo_k(( moc[:, numpy.newaxis], eri.mo_coeff_full[k2][:, :eri.nocc_full[k2]], eri.mo_coeff_full[k2][:, :eri.nocc_full[k2]], moc[:, numpy.newaxis], ), (k, k2, k2, k)).squeeze().trace().real vk /= nk mf = eri.model v_mf = mf.get_veff()[k] - mf.get_j()[k] v_mf = einsum("i,ij,j", moc.conj(), v_mf, moc).real return moe + vk - v_mf
[docs] class IMDS(gw_slow.IMDS): orb_dims = 2 def __init__(self, td, eri=None): """ GW intermediates (k-version/supercell). Args: td: a container with TD solution; eri: a container with electron repulsion integrals; """ gw_slow.AbstractIMDS.__init__(self, td, eri=eri) self.nk = len(self.td._scf.mo_energy) # MF self.nocc = sum(self.eri.nocc) self.o = numpy.concatenate(tuple(e[:nocc] for e, nocc in zip(self.eri.mo_energy, self.eri.nocc))) self.v = numpy.concatenate(tuple(e[nocc:] for e, nocc in zip(self.eri.mo_energy, self.eri.nocc))) # TD nroots, _, k1, k2, o, v = self.td.xy.shape self.td_xy = self.td.xy.transpose(0, 1, 2, 4, 3, 5).reshape(nroots, 2, k1*o, k2*v) self.td_e = self.td.e self.tdm = self.construct_tdm()
[docs] def eri_ov(self, item): result = [] k = numpy.arange(self.nk) for k1 in k: result.append([]) for k2 in k: result[-1].append([]) for k3 in k: result[-1][-1].append([]) for k4 in k: result[-1][-1][-1].append(self.eri.eri_ov(item, (k1, k2, k3, k4))) r = numpy.block(result) return r / len(k)
__getitem__ = eri_ov def __plain_index__(self, p, spec=True): k, kp = p if kp < self.eri.nocc[k]: x = sum(self.eri.nocc[:k]) + kp if spec: return "o", x else: return x else: kp -= self.eri.nocc[k] x = sum(self.eri.nmo[:k]) - sum(self.eri.nocc[:k]) + kp if spec: return "v", x else: return x + self.nocc
[docs] def get_rhs(self, p, components=False): k, kp = p # return self.eri.mo_energy[k][kp] return corrected_moe(self.eri, k, kp)
[docs] def get_sigma_element(self, omega, p, eta, vir_sgn=1): return super(IMDS, self).get_sigma_element(omega, self.__plain_index__(p, spec=False), eta, vir_sgn=vir_sgn)
[docs] def initial_guess(self, p): k, kp = p return self.eri.mo_energy[k][kp]
@property def entire_space(self): assert all(i == self.eri.nmo[0] for i in self.eri.nmo) return [numpy.arange(self.nk), numpy.arange(self.eri.nmo[0])]
kernel = gw_slow.kernel
[docs] class GW(gw_slow.GW): base_imds = IMDS