From d05d6b0cf7ea1a21ae52e835f72a67e99bbd4913 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 29 Jan 2024 14:08:25 +0000 Subject: [PATCH 01/33] Cleaning up branches --- momentGW/pbc/gw.py | 35 +++++- momentGW/pbc/ints.py | 67 ++++++------ momentGW/pbc/tda.py | 255 ++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 321 insertions(+), 36 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 13630c3c..ebcc935c 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -175,6 +175,37 @@ def ao2mo(self, transform=True): return integrals + def build_se_moments(self, nmom_max, integrals, **kwargs): + """Build the moments of the self-energy. + + Parameters + ---------- + nmom_max : int + Maximum moment number to calculate. + integrals : KIntegrals + Density-fitted integrals. + + See functions in `momentGW.rpa` for `kwargs` options. + + Returns + ------- + se_moments_hole : numpy.ndarray + Moments of the hole self-energy at each k-point. If + `self.diagonal_se`, non-diagonal elements are set to zero. + se_moments_part : numpy.ndarray + Moments of the particle self-energy at each k-point. If + `self.diagonal_se`, non-diagonal elements are set to zero. + """ + + if self.polarizability.lower() == "dtda": + tda = dTDA(self, nmom_max, integrals, **kwargs) + return tda.kernel() + elif self.polarizability.lower() == "thc-dtda": + tda = thc.dTDA(self, nmom_max, integrals, **kwargs) + return tda.kernel() + else: + raise NotImplementedError + def solve_dyson(self, se_moments_hole, se_moments_part, se_static, integrals=None): """Solve the Dyson equation due to a self-energy resulting from a list of hole and particle moments, along with a static @@ -334,7 +365,7 @@ def make_rdm1(self, gf=None): @logging.with_timer("Energy") @logging.with_status("Calculating energy") - def energy_hf(self, gf=None, integrals=None): + def energy_hf(self, gf=None, integrals=None, **kwargs): """Calculate the one-body (Hartree--Fock) energy. Parameters @@ -367,7 +398,7 @@ def energy_hf(self, gf=None, integrals=None): "kpq,kpi,kqj->kij", self._scf.get_hcore(), self.mo_coeff.conj(), self.mo_coeff ) rdm1 = self.make_rdm1() - fock = integrals.get_fock(rdm1, h1e) + fock = integrals.get_fock(rdm1, h1e,**kwargs) # Calculate the Hartree--Fock energy at each k-point e_1b = sum(energy.hartree_fock(rdm1[k], fock[k], h1e[k]) for k in self.kpts.loop(1)) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index aca2ad09..43befcdc 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -10,6 +10,8 @@ from pyscf.ao2mo import _ao2mo from pyscf.pbc import tools from scipy.linalg import cholesky +from pyscf.pbc.dft.numint import eval_ao +from pyscf.pbc.dft.gen_grid import get_becke_grids from momentGW import logging, mpi_helper, util from momentGW.ints import Integrals, require_compression_metric @@ -417,38 +419,6 @@ def get_cderi_from_thc(self): self._blocks["Lia"] = Lia self._blocks["Lai"] = Lai - def update_coeffs(self, mo_coeff_g=None, mo_coeff_w=None, mo_occ_w=None): - """ - Update the MO coefficients in-place for the Green's function - and the screened Coulomb interaction. - - Parameters - ---------- - mo_coeff_g : numpy.ndarray, optional - Coefficients corresponding to the Green's function at each - k-point. Default value is `None`. - mo_coeff_w : numpy.ndarray, optional - Coefficients corresponding to the screened Coulomb - interaction at each k-point. Default value is `None`. - mo_occ_w : numpy.ndarray, optional - Occupations corresponding to the screened Coulomb - interaction at each k-point. Default value is `None`. - - Notes - ----- - If `mo_coeff_g` is `None`, the Green's function is assumed to - remain in the basis in which it was originally defined, and - vice-versa for `mo_coeff_w` and `mo_occ_w`. At least one of - `mo_coeff_g` and `mo_coeff_w` must be provided. - """ - return super().update_coeffs( - mo_coeff_g=mo_coeff_g, - mo_coeff_w=mo_coeff_w, - mo_occ_w=mo_occ_w, - ) - - @logging.with_timer("J matrix") - @logging.with_status("Building J matrix") def get_j(self, dm, basis="mo", other=None): """Build the J matrix. @@ -727,6 +697,39 @@ def get_fock(self, dm, h1e, **kwargs): """ return super().get_fock(dm, h1e, **kwargs) + def get_q_ij(self,q,mo_energy_w): + cell = self.with_df.cell + kpts = self.kpts + coords, weights = get_becke_grids(cell, level=5) + qij = np.zeros((len(kpts), self.nmo), dtype=complex) + for ki in kpts.loop(1): + psi_all = eval_ao(cell, coords, kpt=kpts[ki], deriv=1) + psi = psi_all[0] + psi_div = psi_all[1:4] + braket = lib.einsum( + "w,pw,dwq->dpq", weights, psi.T.conj(), psi_div + ) + + num_ao = -1.0j * lib.einsum( + "d,dpq->pq", q, braket + ) + num_mo = np.linalg.multi_dot( + (self.mo_coeff_w[ki][:, self.mo_occ_w[ki] > 0].T.conj(), + num_ao, + self.mo_coeff_w[ki][:, self.mo_occ_w[ki] == 0]) + ) + den = 1/(mo_energy_w[ki][self.mo_occ_w[ki] == 0, None] - mo_energy_w[ki][None,self.mo_occ_w[ki] > 0]) + qij[ki] = (den.T * num_mo).flatten() + + return qij + + + def reciprocal_lattice(self): + """ + Return the reciprocal lattice vectors. + """ + return self.with_df.cell.reciprocal_vectors() + @property def madelung(self): """ diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 0ade92d6..92ce5000 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -4,6 +4,10 @@ import numpy as np import scipy.special +from pyscf import lib +from pyscf.agf2 import mpi_helper +from pyscf.pbc import dft +from pyscf.pbc.gw.krgw_ac import get_qij from momentGW import logging, mpi_helper, util from momentGW.tda import dTDA as MoldTDA @@ -33,6 +37,21 @@ class dTDA(MoldTDA): interaction, respectively. If `None`, use `gw.mo_occ` for both. Default value is `None`. """ + def __init__( + self, + gw, + nmom_max, + integrals, + mo_energy=None, + mo_occ=None, + head_wings=False, + ): + super().__init__(gw, nmom_max, integrals, mo_energy, mo_occ) + self.head_wings = head_wings + if self.head_wings: + q = np.array([1e-3, 0, 0]).reshape(1, 3) + self.q_abs = self.kpts.cell.get_abs_kpts(q) + self.qij = self.integrals.get_q_ij(self.q_abs[0], self.mo_energy_w) @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -49,8 +68,21 @@ def build_dd_moments(self): kpts = self.kpts moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) + + if self.head_wings: + head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + M = np.zeros((self.nkpts), dtype=object) + norm_q_abs = np.linalg.norm(self.q_abs[0]) + + for kb in kpts.loop(1): + M[kb] = self.integrals.Lia[kb, kb] #self.qij[kb]+self.integrals.Lia[kb, kb] + # Get the zeroth order moment for q in kpts.loop(1): + if self.head_wings: + head[q, 0] += (np.sqrt(4.*np.pi) / norm_q_abs)*self.qij[q] + wing[q, 0] += self.integrals.Lia[q, q] for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts @@ -58,6 +90,11 @@ def build_dd_moments(self): # Get the higher order moments for i in range(1, self.nmom_max + 1): for q in kpts.loop(1): + tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) + tmp_head = np.zeros((self.nmo, self.nmo), dtype=complex) + tmp_wings = np.zeros((self.naux[q], self.naux[q]), dtype=complex) + + for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) @@ -73,16 +110,39 @@ def build_dd_moments(self): tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) + if q == 0 and self.head_wings: + head[kb, i] += head[kb, i - 1] * d.ravel() + wing[kb, i] += wing[kb, i - 1] * d.ravel() + + tmp_head += np.dot(M[kb].T.conj(), M[kb]) + tmp_wings += np.dot(wing[kb, i - 1], M[kb].T.conj()) + + tmp = mpi_helper.allreduce(tmp) tmp *= 2.0 tmp /= self.nkpts + tmp_head = mpi_helper.allreduce(tmp_head) + tmp_head *= 2.0 + tmp_head /= self.nkpts + tmp_wings = mpi_helper.allreduce(tmp_wings) + tmp_wings *= 2.0 + tmp_wings /= self.nkpts + for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) - return moments + if q == 0 and self.head_wings: + head[kb, i] += np.dot(head[kb, i - 1], tmp_head) + wing[kb, i] += np.dot(tmp_wings, M[kb]) + + + if self.head_wings: + return {"moments":moments, "head":head, "wing":wing, "M":M} + else: + return moments def kernel(self, exact=False): """ @@ -220,13 +280,34 @@ def build_se_moments(self, moments_dd): eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) eta = np.zeros((self.nkpts, self.nkpts), dtype=object) + if self.head_wings: + cell = self.kpts.cell + cell_vol = cell.vol + total_vol = cell_vol * self.nkpts + q0 = (6*np.pi/total_vol)**(1/3) + ewald = -(2 * q0)/ np.pi + norm_q_abs = np.linalg.norm(self.q_abs[0]) + eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + # Get the moments in (aux|aux) and rotate to (mo|mo) for n in range(self.nmom_max + 1): for q in kpts.loop(1): eta_aux = 0 for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) + if self.head_wings: + eta_aux += np.dot(moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj()) + if q==0: + eta_head[kb, n] += -(np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["head"][kb, n], + self.qij[kb]) + tmp_wing = (np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["wing"][kb, n], + self.qij[kb].conj()) + eta_wings[kb, n] += -(np.sqrt(cell_vol / np.pi ** 3) * q0 ** (2 / 3) * + lib.einsum("Pa,P->a", self.integrals.Lia[kb, kb], tmp_wing)) + # For the wing we have incorporated a factor of 2 for two wings + else: + eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) eta_aux = mpi_helper.allreduce(eta_aux) eta_aux *= 2.0 @@ -242,12 +323,182 @@ def build_se_moments(self, moments_dd): Lp = self.integrals.Lpx[kp, kx][:, :, x] subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) + if self.head_wings and q == 0: + original = eta[kp, q][x, n] + eta[kp, q][x, n] += eta_head[kp, n]*original + eta[kp, q][x, n] += lib.einsum("n, nm-> nm",eta_wings[kp, n], original) + eta[kp, q][x, n] += ewald + + cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) + # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) + + if self.gw.fc: + moments_dd_fc = self.build_dd_moments_fc() + + moments_occ_fc, moments_vir_fc = self.build_se_moments_fc(*moments_dd_fc) + + moments_occ += moments_occ_fc + moments_vir += moments_vir_fc + + cput1 = lib.logger.timer(self.gw, "fc correction", *cput1) + return moments_occ, moments_vir + def build_dd_moments_fc(self): + """ + Build the moments of the "head" (G=0, G'=0) and "wing" + (G=P, G'=0) density-density response. + """ + + kpts = self.kpts + integrals = self.integrals + + # q-point for q->0 finite size correction + qpt = np.array([1e-3, 0.0, 0.0]) + qpt = self.kpts.get_abs_kpts(qpt) + + # 1/sqrt(Ω) * ⟨Ψ_{ik}|e^{iqr}|Ψ_{ak-q}⟩ + qij = self.build_pert_term(qpt) + + # Build the DD moments for the "head" (G=0, G'=0) correction + moments_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=complex) + for k in kpts.loop(1): + d = lib.direct_sum( + "a-i->ia", + self.mo_energy_w[k][self.mo_occ_w[k] == 0], + self.mo_energy_w[k][self.mo_occ_w[k] > 0], + ) + dn = np.ones_like(d) + for n in range(self.nmom_max + 1): + moments_head[k, n] = lib.einsum("ia,ia,ia->", qij[k], qij[k].conj(), dn) + dn *= d + + # Build the DD moments for the "wing" (G=P, G'=0) correction + moments_wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + for k in kpts.loop(1): + d = lib.direct_sum( + "a-i->ia", + self.mo_energy_w[k][self.mo_occ_w[k] == 0], + self.mo_energy_w[k][self.mo_occ_w[k] > 0], + ) + dn = np.ones_like(d) + for n in range(self.nmom_max + 1): + moments_wing[k, n] = lib.einsum( + "Lx,x,x->L", + integrals.Lia[k, k], + qij[k].conj().ravel(), + dn.ravel(), + ) + dn *= d + + moments_head *= -4.0 * np.pi + moments_wing *= -4.0 * np.pi + + return moments_head, moments_wing + + def build_se_moments_fc(self, moments_head, moments_wing): + """ + Build the moments of the self-energy corresponding to the + "wing" (G=P, G'=0) and "head" (G=0, G'=0) density-density + response via convolution. + """ + + kpts = self.kpts + integrals = self.integrals + moms = np.arange(self.nmom_max + 1) + + # Construct the self-energy moments for the "head" (G=0, G'=0) + moments_occ_h = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) + moments_vir_h = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) + for n in moms: + fp = scipy.special.binom(n, moms) + fh = fp * (-1) ** moms + for k in kpts.loop(1): + eo = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] > 0], n - moms) + to = lib.einsum("t,kt,t->", fh, eo, moments_head[k]) + moments_occ_h[k, n] = np.diag([to] * self.nmo) + + ev = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] == 0], n - moms) + tv = lib.einsum("t,kt,t->", fp, ev, moments_head[k]) + moments_vir_h[k, n] = np.diag([tv] * self.nmo) + + # Construct the self-energy moments for the "wing" (G=P, G'=0) + moments_occ_w = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) + moments_vir_w = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) + for n in moms: + fp = scipy.special.binom(n, moms) + fh = fp * (-1) ** moms + for k in kpts.loop(1): + eta = np.zeros( + (self.integrals.nmo_g[k], self.nmom_max + 1, self.nmo), dtype=complex + ) + for t in moms: + eta[:, t] = lib.einsum("L,Lpx->xp", moments_wing[k, t], integrals.Lpx[k, k]) + eta[:, t] += lib.einsum("xpL,L->xp", integrals.Lpx[k, k].conj().T,moments_wing[k, t].conj()) + + eo = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] > 0], n - moms) + to = lib.einsum("t,kt,ktp->p", fh, eo, eta[self.mo_occ_g[k] > 0]) + moments_occ_w[k, n] = np.diag(to) + + ev = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] == 0], n - moms) + tv = lib.einsum("t,kt,ktp->p", fp, ev, eta[self.mo_occ_g[k] == 0]) + moments_vir_w[k, n] = np.diag(tv) + + moments_occ = moments_occ_h + moments_occ_w + moments_vir = moments_vir_h + moments_vir_w + + factor = -2.0 / np.pi * (6.0 * np.pi**2 / self.gw.cell.vol / self.nkpts) ** (1.0 / 3.0) + moments_occ *= factor + moments_vir *= factor + + return moments_occ, moments_vir + + def build_pert_term(self, qpt): + """ + Compute 1/sqrt(Ω) * ⟨Ψ_{ik}|e^{iqr}|Ψ_{ak-q}⟩ at q-point index + q using perturbation theory. + """ + + coords, weights = dft.gen_grid.get_becke_grids(self.gw.cell, level=5) + ngrid = len(coords) + + qij = np.zeros((self.nkpts,), dtype=object) + for k in self.kpts.loop(1): + ao_p = dft.numint.eval_ao(self.gw.cell, coords, kpt=self.kpts[k], deriv=1) + ao, ao_grad = ao_p[0], ao_p[1:4] + + ao_ao_grad = lib.einsum("g,gm,xgn->xmn", weights, ao.conj(), ao_grad) + q_ao_ao_grad = lib.einsum("x,xmn->mn", qpt, ao_ao_grad) * -1.0j + q_mo_mo_grad = lib.einsum( + "mn,mi,na->ia", + q_ao_ao_grad, + self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] > 0].conj(), + self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] == 0], + ) + + d = lib.direct_sum( + "a-i->ia", + self.mo_energy_w[k][self.mo_occ_w[k] == 0], + self.mo_energy_w[k][self.mo_occ_w[k] > 0], + ) + dens = q_mo_mo_grad / d + qij[k] = dens / np.sqrt(self.gw.cell.vol) + + return qij + + + + build_dd_moments_exact = build_dd_moments + + @property + def naux(self): + """Number of auxiliaries.""" + return self.integrals.naux + @property def nov(self): """Get the number of ov states in W.""" From 83e77412dabd302ac2f765eb947b92d37218b1fd Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 8 Feb 2024 18:43:50 +0000 Subject: [PATCH 02/33] Aligning with dRPA and ewald correction added. --- momentGW/pbc/base.py | 2 ++ momentGW/pbc/gw.py | 58 ++++++++++++++++++++++++++++++++++++++++---- momentGW/pbc/ints.py | 8 +++--- momentGW/pbc/tda.py | 16 ++++++++++++ 4 files changed, 75 insertions(+), 9 deletions(-) diff --git a/momentGW/pbc/base.py b/momentGW/pbc/base.py index 62d74d54..757e1cd6 100644 --- a/momentGW/pbc/base.py +++ b/momentGW/pbc/base.py @@ -62,9 +62,11 @@ class BaseKGW(BaseGW): # --- Extra PBC options fc = False + head_wings = False _opts = BaseGW._opts + [ "fc", + "head_wings", ] get_nmo = get_nmo diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index ebcc935c..5c7e8875 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -5,6 +5,8 @@ import numpy as np from dyson import MBLSE, Lehmann, MixedMBLSE +from pyscf.pbc import tools +from functools import reduce from momentGW import energy, logging, util from momentGW.gw import GW @@ -71,23 +73,69 @@ def name(self): @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, mo_coeff=None, mo_energy=None): """ Build the static part of the self-energy, including the Fock matrix. Parameters ---------- - integrals : KIntegrals + integrals : Integrals Integrals object. Returns ------- se_static : numpy.ndarray - Static part of the self-energy at each k-point. If - `self.diagonal_se`, non-diagonal elements are set to zero. + Static part of the self-energy. If `self.diagonal_se`, + non-diagonal elements are set to zero. """ - return super().build_se_static(integrals) + + if mo_energy is None: + mo_energy = self.mo_energy + + if getattr(self._scf, "xc", "hf") == "hf": + se_static = np.zeros_like(self._scf.make_rdm1(mo_coeff=self.mo_coeff)) + if self.fc: + with util.SilentSCF(self._scf): + vmf = self._scf.get_j() - self._scf.get_veff() + dm = self._scf.make_rdm1(mo_coeff=self.mo_coeff) + vk = integrals.get_k(dm, basis="ao") + + s = self.cell.pbc_intor('int1e_ovlp', hermi=1, kpts=self.kpts) + madelung = tools.pbc.madelung(self.cell, self.kpts) + for k in range(len(self.kpts)): + vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) + + se_static = vmf - vk * 0.5 + se_static = util.einsum( + "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff + ) + + else: + with util.SilentSCF(self._scf): + vmf = self._scf.get_j() - self._scf.get_veff() + dm = self._scf.make_rdm1(mo_coeff=self.mo_coeff) + vk = integrals.get_k(dm, basis="ao") + + if self.fc: + s = self.cell.pbc_intor('int1e_ovlp', hermi=1, kpts=self.kpts) + madelung = tools.pbc.madelung(self.cell, self.kpts) + for k in range(len(self.kpts)): + vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) + se_static = vmf - vk * 0.5 + + se_static = util.einsum( + "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff + ) + + + + if self.diagonal_se: + se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) + + se_static += util.einsum("...p,...pq->...pq", mo_energy, np.eye(se_static.shape[-1])) + + return se_static def build_se_moments(self, nmom_max, integrals, **kwargs): """Build the moments of the self-energy. diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 43befcdc..1623a92e 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -400,7 +400,7 @@ def get_cderi_from_thc(self): Lpx[ki, kj] = Lpx_k Lia[ki, kj] = Lia_k - q = self.kpts.member(self.kpts.wrap_around(-self.kpts[q])) + invq = self.kpts.member(self.kpts.wrap_around(-self.kpts[q])) block_switch = util.einsum("Pp,Pq,PQ->Qpq", coll_kj.conj(), coll_ki, cholesky_cou) @@ -410,7 +410,7 @@ def get_cderi_from_thc(self): ) tmp = util.einsum("Lpq,pi,qj->Lij", block_switch, coeffs[0].conj(), coeffs[1]) tmp = tmp.swapaxes(1, 2) - tmp = tmp.reshape(self.naux[q], -1) + tmp = tmp.reshape(self.naux[invq], -1) Lai_k += tmp Lai[ki, kj] = Lai_k @@ -482,7 +482,7 @@ def get_j(self, dm, basis="mo", other=None): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") block = block[0] + block[1] * 1.0j - block = block.reshape(block.shape[0], self.nmo, self.nmo) + block = block.reshape(self.naux_full, self.nmo, self.nmo) b0, b1 = b1, b1 + block.shape[0] buf[b0:b1] += util.einsum("Lpq,pq->L", block, dm[kk].conj()) @@ -494,7 +494,7 @@ def get_j(self, dm, basis="mo", other=None): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") block = block[0] + block[1] * 1.0j - block = block.reshape(block.shape[0], self.nmo, self.nmo) + block = block.reshape(self.naux_full, self.nmo, self.nmo) b0, b1 = b1, b1 + block.shape[0] vj[ki] += util.einsum("Lpq,L->pq", block, buf[b0:b1]) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 92ce5000..aa7ce49d 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -331,6 +331,22 @@ def build_se_moments(self, moments_dd): cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) + for n in range(self.nmom_max + 1): + a = 0 + for q in kpts.loop(1): + for kj in kpts.loop(1, mpi=True): + kx = kpts.member(kpts.wrap_around(kpts[kj] - kpts[q])) + for x in range(self.mo_energy_g[kx].size): + kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) + a += np.sum(eta[kj, kb][x,n]) + print(a.real) + # print("break") + # for q in kpts.loop(1): + # for kj in kpts.loop(1, mpi=True): + # kb = kpts.member(kpts.wrap_around(kpts[kj] - kpts[q])) + # print(np.sum(np.dot(self.integrals.Lia[kj, kb].T.conj(), self.integrals.Lia[kj, kb]))) + + # print(self.integrals.Lia[0,0]) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) From c4dcf82b4831ee957f8f2c93ed48dabf14ec7a20 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 21 Feb 2024 23:18:15 +0000 Subject: [PATCH 03/33] Updates to Head term --- momentGW/pbc/gw.py | 7 +++++ momentGW/pbc/ints.py | 2 +- momentGW/pbc/tda.py | 63 ++++++++++++++++++++++++++------------------ 3 files changed, 45 insertions(+), 27 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 5c7e8875..5456daf1 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -244,6 +244,13 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): Moments of the particle self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ + if self.fc or self.head_wings: + hw = True + else: + hw = False + kwargs = dict( + head_wings=hw, + ) if self.polarizability.lower() == "dtda": tda = dTDA(self, nmom_max, integrals, **kwargs) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 1623a92e..ab81e82d 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -701,7 +701,7 @@ def get_q_ij(self,q,mo_energy_w): cell = self.with_df.cell kpts = self.kpts coords, weights = get_becke_grids(cell, level=5) - qij = np.zeros((len(kpts), self.nmo), dtype=complex) + qij = np.zeros((len(kpts), self.nocc_w[0] * self.nvir_w[0]), dtype=complex) for ki in kpts.loop(1): psi_all = eval_ao(cell, coords, kpt=kpts[ki], deriv=1) psi = psi_all[0] diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index aa7ce49d..0f2939fe 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -51,7 +51,14 @@ def __init__( if self.head_wings: q = np.array([1e-3, 0, 0]).reshape(1, 3) self.q_abs = self.kpts.cell.get_abs_kpts(q) - self.qij = self.integrals.get_q_ij(self.q_abs[0], self.mo_energy_w) + self.qij = self.build_pert_term(self.q_abs[0]) + zeroth = 0j + for q in self.kpts.loop(1): + norm_q_abs = np.linalg.norm(self.q_abs[0]) + zeroth += -(4. * np.pi / norm_q_abs**2) * np.dot(self.qij[q].conj(),self.qij[q]) + a = np.sum(self.qij[q].conj()*self.qij[q]) + print(np.allclose(a, np.dot(self.qij[q].conj(),self.qij[q]))) + @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -74,27 +81,22 @@ def build_dd_moments(self): wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) M = np.zeros((self.nkpts), dtype=object) norm_q_abs = np.linalg.norm(self.q_abs[0]) - for kb in kpts.loop(1): M[kb] = self.integrals.Lia[kb, kb] #self.qij[kb]+self.integrals.Lia[kb, kb] # Get the zeroth order moment for q in kpts.loop(1): - if self.head_wings: - head[q, 0] += (np.sqrt(4.*np.pi) / norm_q_abs)*self.qij[q] - wing[q, 0] += self.integrals.Lia[q, q] for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts + if self.head_wings: + head[q, 0] += (np.sqrt(4. * np.pi) / norm_q_abs) * self.qij[q].conj() + # wing[q, 0] += self.integrals.Lia[q, q] + # Get the higher order moments for i in range(1, self.nmom_max + 1): for q in kpts.loop(1): - tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) - tmp_head = np.zeros((self.nmo, self.nmo), dtype=complex) - tmp_wings = np.zeros((self.naux[q], self.naux[q]), dtype=complex) - - for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) @@ -103,19 +105,21 @@ def build_dd_moments(self): (self.mo_occ_w[kj], self.mo_occ_w[kb]), ) moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] + if self.head_wings and q==0: + head[kb, i] += head[kb, i - 1] * d.ravel() + # wing[kb, i] += wing[kb, i - 1] * d.ravel() tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) + tmp_head = np.zeros((self.naux[q]), dtype=complex) + tmp_wings = np.zeros((self.naux[q], self.naux[q]), dtype=complex) for ki in kpts.loop(1, mpi=True): ka = kpts.member(kpts.wrap_around(kpts[q] + kpts[ki])) tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) if q == 0 and self.head_wings: - head[kb, i] += head[kb, i - 1] * d.ravel() - wing[kb, i] += wing[kb, i - 1] * d.ravel() - - tmp_head += np.dot(M[kb].T.conj(), M[kb]) - tmp_wings += np.dot(wing[kb, i - 1], M[kb].T.conj()) + tmp_head += lib.einsum("a,aP->P",head[kb, i - 1], M[kb].T.conj()) + # tmp_wings += np.dot(wing[kb, i - 1], M[kb].T.conj()) tmp = mpi_helper.allreduce(tmp) @@ -125,9 +129,10 @@ def build_dd_moments(self): tmp_head = mpi_helper.allreduce(tmp_head) tmp_head *= 2.0 tmp_head /= self.nkpts - tmp_wings = mpi_helper.allreduce(tmp_wings) - tmp_wings *= 2.0 - tmp_wings /= self.nkpts + + # tmp_wings = mpi_helper.allreduce(tmp_wings) + # tmp_wings *= 2.0 + # tmp_wings /= self.nkpts for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) @@ -135,8 +140,8 @@ def build_dd_moments(self): moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) if q == 0 and self.head_wings: - head[kb, i] += np.dot(head[kb, i - 1], tmp_head) - wing[kb, i] += np.dot(tmp_wings, M[kb]) + head[kb, i] += lib.einsum("P,Pa->a",tmp_head, M[kb]) + # wing[kb, i] += np.dot(tmp_wings, M[kb]) if self.head_wings: @@ -299,12 +304,12 @@ def build_se_moments(self, moments_dd): if self.head_wings: eta_aux += np.dot(moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj()) if q==0: - eta_head[kb, n] += -(np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["head"][kb, n], + eta_head[kb, n] += -(np.sqrt(4. * np.pi) / norm_q_abs) * np.sum(moments_dd["head"][kb, n]* self.qij[kb]) - tmp_wing = (np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["wing"][kb, n], - self.qij[kb].conj()) - eta_wings[kb, n] += -(np.sqrt(cell_vol / np.pi ** 3) * q0 ** (2 / 3) * - lib.einsum("Pa,P->a", self.integrals.Lia[kb, kb], tmp_wing)) + # tmp_wing = (np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["wing"][kb, n], + # self.qij[kb].conj()) + # eta_wings[kb, n] += -(np.sqrt(cell_vol / np.pi ** 3) * q0 ** (2 / 3) * + # lib.einsum("Pa,P->a", self.integrals.Lia[kb, kb], tmp_wing)) # For the wing we have incorporated a factor of 2 for two wings else: eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) @@ -325,6 +330,12 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) if self.head_wings and q == 0: original = eta[kp, q][x, n] + # eta[kp, q][x, n] += eta_head[kp, n]*original + # print(eta_head[kp, n]/np.max(original)) + # if n==0: + # print("head",eta_head[kp, n]) + # print("body",eta[kp, q][x, n]) + # eta[kp, q][x, n] += lib.einsum("n, nm-> nm",eta_wings[kp, n], original) eta[kp, q][x, n] += eta_head[kp, n]*original eta[kp, q][x, n] += lib.einsum("n, nm-> nm",eta_wings[kp, n], original) eta[kp, q][x, n] += ewald @@ -502,7 +513,7 @@ def build_pert_term(self, qpt): self.mo_energy_w[k][self.mo_occ_w[k] > 0], ) dens = q_mo_mo_grad / d - qij[k] = dens / np.sqrt(self.gw.cell.vol) + qij[k] = dens.flatten() / np.sqrt(self.gw.cell.vol) return qij From 2b0bb39c958c58c26912d3ad4690111784621c95 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 29 Feb 2024 17:11:32 +0000 Subject: [PATCH 04/33] kpt changes for rounding --- momentGW/pbc/tda.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 0f2939fe..c6608fc0 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -52,12 +52,12 @@ def __init__( q = np.array([1e-3, 0, 0]).reshape(1, 3) self.q_abs = self.kpts.cell.get_abs_kpts(q) self.qij = self.build_pert_term(self.q_abs[0]) - zeroth = 0j - for q in self.kpts.loop(1): - norm_q_abs = np.linalg.norm(self.q_abs[0]) - zeroth += -(4. * np.pi / norm_q_abs**2) * np.dot(self.qij[q].conj(),self.qij[q]) - a = np.sum(self.qij[q].conj()*self.qij[q]) - print(np.allclose(a, np.dot(self.qij[q].conj(),self.qij[q]))) + # zeroth = 0j + # for q in self.kpts.loop(1): + # norm_q_abs = np.linalg.norm(self.q_abs[0]) + # zeroth += -(4. * np.pi / norm_q_abs**2) * np.dot(self.qij[q].conj(),self.qij[q]) + # a = np.sum(self.qij[q].conj()*self.qij[q]) + # print(np.allclose(a, np.dot(self.qij[q].conj(),self.qij[q]))) @logging.with_timer("Density-density moments") @@ -330,7 +330,7 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) if self.head_wings and q == 0: original = eta[kp, q][x, n] - # eta[kp, q][x, n] += eta_head[kp, n]*original + eta[kp, q][x, n] += eta_head[kp, n]*original # print(eta_head[kp, n]/np.max(original)) # if n==0: # print("head",eta_head[kp, n]) From 779a964456025393c4056a6b2f7b9c09abf97d8e Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 11 Mar 2024 17:17:37 +0000 Subject: [PATCH 05/33] Changes hopefully fixing the naux issues and subsequent Hermiticity issue. --- momentGW/pbc/ints.py | 2 -- momentGW/pbc/tda.py | 11 +++-------- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index ab81e82d..c2458f50 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -569,7 +569,6 @@ def get_k(self, dm, basis="mo", ewald=False): ) for ki in self.kpts.loop(1, mpi=True): b1 = 0 - kk = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) for block in self.with_df.sr_loop((ki, kk), compact=False): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") @@ -582,7 +581,6 @@ def get_k(self, dm, basis="mo", ewald=False): for ki in self.kpts.loop(1, mpi=True): b1 = 0 - kk = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) for block in self.with_df.sr_loop((kk, ki), compact=False): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index c6608fc0..d1c4afb8 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -52,12 +52,6 @@ def __init__( q = np.array([1e-3, 0, 0]).reshape(1, 3) self.q_abs = self.kpts.cell.get_abs_kpts(q) self.qij = self.build_pert_term(self.q_abs[0]) - # zeroth = 0j - # for q in self.kpts.loop(1): - # norm_q_abs = np.linalg.norm(self.q_abs[0]) - # zeroth += -(4. * np.pi / norm_q_abs**2) * np.dot(self.qij[q].conj(),self.qij[q]) - # a = np.sum(self.qij[q].conj()*self.qij[q]) - # print(np.allclose(a, np.dot(self.qij[q].conj(),self.qij[q]))) @logging.with_timer("Density-density moments") @@ -81,8 +75,9 @@ def build_dd_moments(self): wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) M = np.zeros((self.nkpts), dtype=object) norm_q_abs = np.linalg.norm(self.q_abs[0]) - for kb in kpts.loop(1): - M[kb] = self.integrals.Lia[kb, kb] #self.qij[kb]+self.integrals.Lia[kb, kb] + for kj in kpts.loop(1, mpi=True): + kb = kpts.member(kpts.wrap_around(kpts[0] + kpts[kj])) + M[kb] = self.integrals.Lia[kb, kb] # Get the zeroth order moment for q in kpts.loop(1): From 12929dc828b5f6e7d8fc53385f3cf95a416d00cd Mon Sep 17 00:00:00 2001 From: mkakcl Date: Tue, 19 Mar 2024 18:21:50 +0000 Subject: [PATCH 06/33] Initial wings implementation. --- momentGW/pbc/tda.py | 70 +++++++++++++-------------------------------- 1 file changed, 20 insertions(+), 50 deletions(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index d1c4afb8..856397ca 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -72,7 +72,6 @@ def build_dd_moments(self): if self.head_wings: head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) M = np.zeros((self.nkpts), dtype=object) norm_q_abs = np.linalg.norm(self.q_abs[0]) for kj in kpts.loop(1, mpi=True): @@ -87,7 +86,6 @@ def build_dd_moments(self): if self.head_wings: head[q, 0] += (np.sqrt(4. * np.pi) / norm_q_abs) * self.qij[q].conj() - # wing[q, 0] += self.integrals.Lia[q, q] # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -102,11 +100,10 @@ def build_dd_moments(self): moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] if self.head_wings and q==0: head[kb, i] += head[kb, i - 1] * d.ravel() - # wing[kb, i] += wing[kb, i - 1] * d.ravel() tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) tmp_head = np.zeros((self.naux[q]), dtype=complex) - tmp_wings = np.zeros((self.naux[q], self.naux[q]), dtype=complex) + for ki in kpts.loop(1, mpi=True): ka = kpts.member(kpts.wrap_around(kpts[q] + kpts[ki])) @@ -114,7 +111,6 @@ def build_dd_moments(self): if q == 0 and self.head_wings: tmp_head += lib.einsum("a,aP->P",head[kb, i - 1], M[kb].T.conj()) - # tmp_wings += np.dot(wing[kb, i - 1], M[kb].T.conj()) tmp = mpi_helper.allreduce(tmp) @@ -125,10 +121,6 @@ def build_dd_moments(self): tmp_head *= 2.0 tmp_head /= self.nkpts - # tmp_wings = mpi_helper.allreduce(tmp_wings) - # tmp_wings *= 2.0 - # tmp_wings /= self.nkpts - for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) @@ -136,11 +128,10 @@ def build_dd_moments(self): if q == 0 and self.head_wings: head[kb, i] += lib.einsum("P,Pa->a",tmp_head, M[kb]) - # wing[kb, i] += np.dot(tmp_wings, M[kb]) if self.head_wings: - return {"moments":moments, "head":head, "wing":wing, "M":M} + return {"moments":moments, "head":head} else: return moments @@ -285,7 +276,6 @@ def build_se_moments(self, moments_dd): cell_vol = cell.vol total_vol = cell_vol * self.nkpts q0 = (6*np.pi/total_vol)**(1/3) - ewald = -(2 * q0)/ np.pi norm_q_abs = np.linalg.norm(self.q_abs[0]) eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) @@ -301,11 +291,7 @@ def build_se_moments(self, moments_dd): if q==0: eta_head[kb, n] += -(np.sqrt(4. * np.pi) / norm_q_abs) * np.sum(moments_dd["head"][kb, n]* self.qij[kb]) - # tmp_wing = (np.sqrt(4. * np.pi) / norm_q_abs) * np.dot(moments_dd["wing"][kb, n], - # self.qij[kb].conj()) - # eta_wings[kb, n] += -(np.sqrt(cell_vol / np.pi ** 3) * q0 ** (2 / 3) * - # lib.einsum("Pa,P->a", self.integrals.Lia[kb, kb], tmp_wing)) - # For the wing we have incorporated a factor of 2 for two wings + eta_wings[kb,n] += (np.sqrt(4. * np.pi) / norm_q_abs) * lib.einsum("Pa,a->P", moments_dd["moments"][q, kb, n], self.qij[kb]) else: eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) @@ -325,48 +311,32 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) if self.head_wings and q == 0: original = eta[kp, q][x, n] - eta[kp, q][x, n] += eta_head[kp, n]*original - # print(eta_head[kp, n]/np.max(original)) - # if n==0: - # print("head",eta_head[kp, n]) - # print("body",eta[kp, q][x, n]) - # eta[kp, q][x, n] += lib.einsum("n, nm-> nm",eta_wings[kp, n], original) - eta[kp, q][x, n] += eta_head[kp, n]*original - eta[kp, q][x, n] += lib.einsum("n, nm-> nm",eta_wings[kp, n], original) - eta[kp, q][x, n] += ewald - cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) + eta[kp, q][x, n] += (2/np.pi)*(q0)*eta_head[kp, n]*original - for n in range(self.nmom_max + 1): - a = 0 - for q in kpts.loop(1): - for kj in kpts.loop(1, mpi=True): - kx = kpts.member(kpts.wrap_around(kpts[kj] - kpts[q])) - for x in range(self.mo_energy_g[kx].size): - kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - a += np.sum(eta[kj, kb][x,n]) - print(a.real) - # print("break") - # for q in kpts.loop(1): - # for kj in kpts.loop(1, mpi=True): - # kb = kpts.member(kpts.wrap_around(kpts[kj] - kpts[q])) - # print(np.sum(np.dot(self.integrals.Lia[kj, kb].T.conj(), self.integrals.Lia[kj, kb]))) + wing_tmp = lib.einsum("Pp,P->p", Lp, eta_wings[kp, n]) + wing_tmp = wing_tmp.real*2 + wing_tmp *= -(np.sqrt(cell_vol/ (4*(np.pi** 3))) * q0**2) - # print(self.integrals.Lia[0,0]) + eta[kp, q][x, n] += lib.einsum("p,pq->pq", wing_tmp, original) - # Construct the self-energy moments - moments_occ, moments_vir = self.convolve(eta) + cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) - if self.gw.fc: - moments_dd_fc = self.build_dd_moments_fc() + # Construct the self-energy moments + moments_occ, moments_vir = self.convolve(eta) - moments_occ_fc, moments_vir_fc = self.build_se_moments_fc(*moments_dd_fc) - moments_occ += moments_occ_fc - moments_vir += moments_vir_fc - cput1 = lib.logger.timer(self.gw, "fc correction", *cput1) + # if self.gw.fc: + # moments_dd_fc = self.build_dd_moments_fc() + # + # moments_occ_fc, moments_vir_fc = self.build_se_moments_fc(*moments_dd_fc) + # + # moments_occ += moments_occ_fc + # moments_vir += moments_vir_fc + # + # cput1 = lib.logger.timer(self.gw, "fc correction", *cput1) return moments_occ, moments_vir From 068374803d1dcf1fd7bbefbdca254887417b53b7 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Tue, 19 Mar 2024 19:07:37 +0000 Subject: [PATCH 07/33] Small fix --- momentGW/pbc/tda.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 856397ca..0e9b306d 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -275,7 +275,7 @@ def build_se_moments(self, moments_dd): cell = self.kpts.cell cell_vol = cell.vol total_vol = cell_vol * self.nkpts - q0 = (6*np.pi/total_vol)**(1/3) + q0 = (6*np.pi**2/total_vol)**(1/3) norm_q_abs = np.linalg.norm(self.q_abs[0]) eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) From f7a844e3d99a4a3b344a0155dad407f303f9a9b3 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 21 Mar 2024 17:28:24 +0000 Subject: [PATCH 08/33] Head and wings cleanup --- momentGW/pbc/base.py | 2 - momentGW/pbc/gw.py | 8 +- momentGW/pbc/tda.py | 172 ++++++------------------------------------- 3 files changed, 28 insertions(+), 154 deletions(-) diff --git a/momentGW/pbc/base.py b/momentGW/pbc/base.py index 757e1cd6..62d74d54 100644 --- a/momentGW/pbc/base.py +++ b/momentGW/pbc/base.py @@ -62,11 +62,9 @@ class BaseKGW(BaseGW): # --- Extra PBC options fc = False - head_wings = False _opts = BaseGW._opts + [ "fc", - "head_wings", ] get_nmo = get_nmo diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 5456daf1..2fcabb75 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -244,12 +244,12 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): Moments of the particle self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - if self.fc or self.head_wings: - hw = True + if self.fc: + fc = True else: - hw = False + fc = False kwargs = dict( - head_wings=hw, + fc=fc, ) if self.polarizability.lower() == "dtda": diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 0e9b306d..ec2b9951 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -44,11 +44,11 @@ def __init__( integrals, mo_energy=None, mo_occ=None, - head_wings=False, + fc=False, ): super().__init__(gw, nmom_max, integrals, mo_energy, mo_occ) - self.head_wings = head_wings - if self.head_wings: + self.fc = fc + if self.fc: q = np.array([1e-3, 0, 0]).reshape(1, 3) self.q_abs = self.kpts.cell.get_abs_kpts(q) self.qij = self.build_pert_term(self.q_abs[0]) @@ -70,13 +70,8 @@ def build_dd_moments(self): moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) - if self.head_wings: + if self.fc: head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - M = np.zeros((self.nkpts), dtype=object) - norm_q_abs = np.linalg.norm(self.q_abs[0]) - for kj in kpts.loop(1, mpi=True): - kb = kpts.member(kpts.wrap_around(kpts[0] + kpts[kj])) - M[kb] = self.integrals.Lia[kb, kb] # Get the zeroth order moment for q in kpts.loop(1): @@ -84,8 +79,8 @@ def build_dd_moments(self): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts - if self.head_wings: - head[q, 0] += (np.sqrt(4. * np.pi) / norm_q_abs) * self.qij[q].conj() + if self.fc: + head[q, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[q].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -98,7 +93,7 @@ def build_dd_moments(self): (self.mo_occ_w[kj], self.mo_occ_w[kb]), ) moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] - if self.head_wings and q==0: + if self.fc and q==0: head[kb, i] += head[kb, i - 1] * d.ravel() tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) @@ -109,8 +104,9 @@ def build_dd_moments(self): tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) - if q == 0 and self.head_wings: - tmp_head += lib.einsum("a,aP->P",head[kb, i - 1], M[kb].T.conj()) + if q == 0 and self.fc: + tmp_head += lib.einsum("a,aP->P",head[kb, i - 1], + self.integrals.Lia[kb, kb].T.conj()) tmp = mpi_helper.allreduce(tmp) @@ -126,11 +122,12 @@ def build_dd_moments(self): moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) - if q == 0 and self.head_wings: - head[kb, i] += lib.einsum("P,Pa->a",tmp_head, M[kb]) + if q == 0 and self.fc: + head[kb, i] += lib.einsum("P,Pa->a", tmp_head, + self.integrals.Lia[kb, kb]) - if self.head_wings: + if self.fc: return {"moments":moments, "head":head} else: return moments @@ -271,9 +268,8 @@ def build_se_moments(self, moments_dd): eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) eta = np.zeros((self.nkpts, self.nkpts), dtype=object) - if self.head_wings: - cell = self.kpts.cell - cell_vol = cell.vol + if self.fc: + cell_vol = self.kpts.cell.vol total_vol = cell_vol * self.nkpts q0 = (6*np.pi**2/total_vol)**(1/3) norm_q_abs = np.linalg.norm(self.q_abs[0]) @@ -286,12 +282,15 @@ def build_se_moments(self, moments_dd): eta_aux = 0 for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - if self.head_wings: + if self.fc: eta_aux += np.dot(moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj()) if q==0: - eta_head[kb, n] += -(np.sqrt(4. * np.pi) / norm_q_abs) * np.sum(moments_dd["head"][kb, n]* - self.qij[kb]) - eta_wings[kb,n] += (np.sqrt(4. * np.pi) / norm_q_abs) * lib.einsum("Pa,a->P", moments_dd["moments"][q, kb, n], self.qij[kb]) + eta_head[kb, n] += (-(np.sqrt(4. * np.pi) / norm_q_abs) * + np.sum(moments_dd["head"][kb, n] * + self.qij[kb])) + eta_wings[kb, n] += ((np.sqrt(4. * np.pi) / norm_q_abs) * + lib.einsum("Pa,a->P", moments_dd["moments"][q, kb, n], + self.qij[kb])) else: eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) @@ -309,7 +308,7 @@ def build_se_moments(self, moments_dd): Lp = self.integrals.Lpx[kp, kx][:, :, x] subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) - if self.head_wings and q == 0: + if self.fc and q == 0: original = eta[kp, q][x, n] eta[kp, q][x, n] += (2/np.pi)*(q0)*eta_head[kp, n]*original @@ -322,131 +321,10 @@ def build_se_moments(self, moments_dd): cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) - # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) - - # if self.gw.fc: - # moments_dd_fc = self.build_dd_moments_fc() - # - # moments_occ_fc, moments_vir_fc = self.build_se_moments_fc(*moments_dd_fc) - # - # moments_occ += moments_occ_fc - # moments_vir += moments_vir_fc - # - # cput1 = lib.logger.timer(self.gw, "fc correction", *cput1) - - return moments_occ, moments_vir - - def build_dd_moments_fc(self): - """ - Build the moments of the "head" (G=0, G'=0) and "wing" - (G=P, G'=0) density-density response. - """ - - kpts = self.kpts - integrals = self.integrals - - # q-point for q->0 finite size correction - qpt = np.array([1e-3, 0.0, 0.0]) - qpt = self.kpts.get_abs_kpts(qpt) - - # 1/sqrt(Ω) * ⟨Ψ_{ik}|e^{iqr}|Ψ_{ak-q}⟩ - qij = self.build_pert_term(qpt) - - # Build the DD moments for the "head" (G=0, G'=0) correction - moments_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=complex) - for k in kpts.loop(1): - d = lib.direct_sum( - "a-i->ia", - self.mo_energy_w[k][self.mo_occ_w[k] == 0], - self.mo_energy_w[k][self.mo_occ_w[k] > 0], - ) - dn = np.ones_like(d) - for n in range(self.nmom_max + 1): - moments_head[k, n] = lib.einsum("ia,ia,ia->", qij[k], qij[k].conj(), dn) - dn *= d - - # Build the DD moments for the "wing" (G=P, G'=0) correction - moments_wing = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - for k in kpts.loop(1): - d = lib.direct_sum( - "a-i->ia", - self.mo_energy_w[k][self.mo_occ_w[k] == 0], - self.mo_energy_w[k][self.mo_occ_w[k] > 0], - ) - dn = np.ones_like(d) - for n in range(self.nmom_max + 1): - moments_wing[k, n] = lib.einsum( - "Lx,x,x->L", - integrals.Lia[k, k], - qij[k].conj().ravel(), - dn.ravel(), - ) - dn *= d - - moments_head *= -4.0 * np.pi - moments_wing *= -4.0 * np.pi - - return moments_head, moments_wing - - def build_se_moments_fc(self, moments_head, moments_wing): - """ - Build the moments of the self-energy corresponding to the - "wing" (G=P, G'=0) and "head" (G=0, G'=0) density-density - response via convolution. - """ - - kpts = self.kpts - integrals = self.integrals - moms = np.arange(self.nmom_max + 1) - - # Construct the self-energy moments for the "head" (G=0, G'=0) - moments_occ_h = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) - moments_vir_h = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) - for n in moms: - fp = scipy.special.binom(n, moms) - fh = fp * (-1) ** moms - for k in kpts.loop(1): - eo = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] > 0], n - moms) - to = lib.einsum("t,kt,t->", fh, eo, moments_head[k]) - moments_occ_h[k, n] = np.diag([to] * self.nmo) - - ev = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] == 0], n - moms) - tv = lib.einsum("t,kt,t->", fp, ev, moments_head[k]) - moments_vir_h[k, n] = np.diag([tv] * self.nmo) - - # Construct the self-energy moments for the "wing" (G=P, G'=0) - moments_occ_w = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) - moments_vir_w = np.zeros((self.nkpts, self.nmom_max + 1, self.nmo, self.nmo), dtype=complex) - for n in moms: - fp = scipy.special.binom(n, moms) - fh = fp * (-1) ** moms - for k in kpts.loop(1): - eta = np.zeros( - (self.integrals.nmo_g[k], self.nmom_max + 1, self.nmo), dtype=complex - ) - for t in moms: - eta[:, t] = lib.einsum("L,Lpx->xp", moments_wing[k, t], integrals.Lpx[k, k]) - eta[:, t] += lib.einsum("xpL,L->xp", integrals.Lpx[k, k].conj().T,moments_wing[k, t].conj()) - - eo = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] > 0], n - moms) - to = lib.einsum("t,kt,ktp->p", fh, eo, eta[self.mo_occ_g[k] > 0]) - moments_occ_w[k, n] = np.diag(to) - - ev = np.power.outer(self.mo_energy_g[k][self.mo_occ_g[k] == 0], n - moms) - tv = lib.einsum("t,kt,ktp->p", fp, ev, eta[self.mo_occ_g[k] == 0]) - moments_vir_w[k, n] = np.diag(tv) - - moments_occ = moments_occ_h + moments_occ_w - moments_vir = moments_vir_h + moments_vir_w - - factor = -2.0 / np.pi * (6.0 * np.pi**2 / self.gw.cell.vol / self.nkpts) ** (1.0 / 3.0) - moments_occ *= factor - moments_vir *= factor - return moments_occ, moments_vir def build_pert_term(self, qpt): @@ -482,8 +360,6 @@ def build_pert_term(self, qpt): return qij - - build_dd_moments_exact = build_dd_moments @property From 3431d181585d5d6180af30e98e5f4b3bdaa446e1 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 3 Apr 2024 15:42:54 +0100 Subject: [PATCH 09/33] Potential MPI fix --- momentGW/pbc/tda.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index ec2b9951..0404c71c 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -80,7 +80,7 @@ def build_dd_moments(self): moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts if self.fc: - head[q, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[q].conj() + head[kj, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[kj].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -94,7 +94,7 @@ def build_dd_moments(self): ) moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] if self.fc and q==0: - head[kb, i] += head[kb, i - 1] * d.ravel() + head[kj, i] += head[kj, i - 1] * d.ravel() tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) tmp_head = np.zeros((self.naux[q]), dtype=complex) @@ -105,7 +105,7 @@ def build_dd_moments(self): tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) if q == 0 and self.fc: - tmp_head += lib.einsum("a,aP->P",head[kb, i - 1], + tmp_head += lib.einsum("a,aP->P",head[kj, i - 1], self.integrals.Lia[kb, kb].T.conj()) @@ -123,8 +123,8 @@ def build_dd_moments(self): moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) if q == 0 and self.fc: - head[kb, i] += lib.einsum("P,Pa->a", tmp_head, - self.integrals.Lia[kb, kb]) + head[kj, i] += lib.einsum("P,Pa->a", tmp_head, + self.integrals.Lia[kj, kj]) if self.fc: From fdd0b30fabf1fecb2b3075f234e03a8ebb71cb85 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 3 Apr 2024 19:03:42 +0100 Subject: [PATCH 10/33] More cleanup --- momentGW/pbc/gw.py | 7 ++----- momentGW/pbc/tda.py | 2 +- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 2fcabb75..6baf70d8 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -73,7 +73,7 @@ def name(self): @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals, mo_coeff=None, mo_energy=None): + def build_se_static(self, integrals): """ Build the static part of the self-energy, including the Fock matrix. @@ -90,9 +90,6 @@ def build_se_static(self, integrals, mo_coeff=None, mo_energy=None): non-diagonal elements are set to zero. """ - if mo_energy is None: - mo_energy = self.mo_energy - if getattr(self._scf, "xc", "hf") == "hf": se_static = np.zeros_like(self._scf.make_rdm1(mo_coeff=self.mo_coeff)) if self.fc: @@ -133,7 +130,7 @@ def build_se_static(self, integrals, mo_coeff=None, mo_energy=None): if self.diagonal_se: se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) - se_static += util.einsum("...p,...pq->...pq", mo_energy, np.eye(se_static.shape[-1])) + se_static += util.einsum("...p,...pq->...pq", self.mo_energy, np.eye(se_static.shape[-1])) return se_static diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 0404c71c..d346b4f5 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -80,7 +80,7 @@ def build_dd_moments(self): moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts if self.fc: - head[kj, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[kj].conj() + head[q, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[q].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): From aff4dc55855ce09f88085bef0c4721588495b3e6 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 17:32:08 +0100 Subject: [PATCH 11/33] A few merge issues resolved --- momentGW/pbc/ints.py | 7 +++++-- momentGW/pbc/tda.py | 3 --- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index c2458f50..ceef29b5 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -482,7 +482,7 @@ def get_j(self, dm, basis="mo", other=None): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") block = block[0] + block[1] * 1.0j - block = block.reshape(self.naux_full, self.nmo, self.nmo) + block = block.reshape(block.shape[0], self.nmo, self.nmo) b0, b1 = b1, b1 + block.shape[0] buf[b0:b1] += util.einsum("Lpq,pq->L", block, dm[kk].conj()) @@ -494,7 +494,7 @@ def get_j(self, dm, basis="mo", other=None): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") block = block[0] + block[1] * 1.0j - block = block.reshape(self.naux_full, self.nmo, self.nmo) + block = block.reshape(block.shape[0], self.nmo, self.nmo) b0, b1 = b1, b1 + block.shape[0] vj[ki] += util.einsum("Lpq,L->pq", block, buf[b0:b1]) @@ -569,6 +569,8 @@ def get_k(self, dm, basis="mo", ewald=False): ) for ki in self.kpts.loop(1, mpi=True): b1 = 0 + kk = self.kpts.member(self.kpts.wrap_around( + self.kpts[q] + self.kpts[ki])) for block in self.with_df.sr_loop((ki, kk), compact=False): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") @@ -581,6 +583,7 @@ def get_k(self, dm, basis="mo", ewald=False): for ki in self.kpts.loop(1, mpi=True): b1 = 0 + kk = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) for block in self.with_df.sr_loop((kk, ki), compact=False): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index d346b4f5..a4d7c4ff 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -53,7 +53,6 @@ def __init__( self.q_abs = self.kpts.cell.get_abs_kpts(q) self.qij = self.build_pert_term(self.q_abs[0]) - @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") def build_dd_moments(self): @@ -319,8 +318,6 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += lib.einsum("p,pq->pq", wing_tmp, original) - cput1 = lib.logger.timer(self.gw, "rotating DD moments", *cput0) - # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) From ac1ad7effc97e8ff4de89362ba5cf8d4bdce9443 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 17:38:29 +0100 Subject: [PATCH 12/33] Linting --- momentGW/pbc/gw.py | 11 +++++------ momentGW/pbc/ints.py | 31 ++++++++++++++--------------- momentGW/pbc/tda.py | 47 ++++++++++++++++++++++---------------------- 3 files changed, 44 insertions(+), 45 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 6baf70d8..92a7e401 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -3,10 +3,11 @@ periodic systems. """ +from functools import reduce + import numpy as np from dyson import MBLSE, Lehmann, MixedMBLSE from pyscf.pbc import tools -from functools import reduce from momentGW import energy, logging, util from momentGW.gw import GW @@ -98,7 +99,7 @@ def build_se_static(self, integrals): dm = self._scf.make_rdm1(mo_coeff=self.mo_coeff) vk = integrals.get_k(dm, basis="ao") - s = self.cell.pbc_intor('int1e_ovlp', hermi=1, kpts=self.kpts) + s = self.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts) madelung = tools.pbc.madelung(self.cell, self.kpts) for k in range(len(self.kpts)): vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) @@ -115,7 +116,7 @@ def build_se_static(self, integrals): vk = integrals.get_k(dm, basis="ao") if self.fc: - s = self.cell.pbc_intor('int1e_ovlp', hermi=1, kpts=self.kpts) + s = self.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts) madelung = tools.pbc.madelung(self.cell, self.kpts) for k in range(len(self.kpts)): vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) @@ -125,8 +126,6 @@ def build_se_static(self, integrals): "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff ) - - if self.diagonal_se: se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) @@ -450,7 +449,7 @@ def energy_hf(self, gf=None, integrals=None, **kwargs): "kpq,kpi,kqj->kij", self._scf.get_hcore(), self.mo_coeff.conj(), self.mo_coeff ) rdm1 = self.make_rdm1() - fock = integrals.get_fock(rdm1, h1e,**kwargs) + fock = integrals.get_fock(rdm1, h1e, **kwargs) # Calculate the Hartree--Fock energy at each k-point e_1b = sum(energy.hartree_fock(rdm1[k], fock[k], h1e[k]) for k in self.kpts.loop(1)) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index ceef29b5..62b09ddd 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -9,9 +9,9 @@ from pyscf import lib from pyscf.ao2mo import _ao2mo from pyscf.pbc import tools -from scipy.linalg import cholesky -from pyscf.pbc.dft.numint import eval_ao from pyscf.pbc.dft.gen_grid import get_becke_grids +from pyscf.pbc.dft.numint import eval_ao +from scipy.linalg import cholesky from momentGW import logging, mpi_helper, util from momentGW.ints import Integrals, require_compression_metric @@ -569,8 +569,7 @@ def get_k(self, dm, basis="mo", ewald=False): ) for ki in self.kpts.loop(1, mpi=True): b1 = 0 - kk = self.kpts.member(self.kpts.wrap_around( - self.kpts[q] + self.kpts[ki])) + kk = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) for block in self.with_df.sr_loop((ki, kk), compact=False): if block[2] == -1: raise NotImplementedError("Low dimensional integrals") @@ -698,7 +697,7 @@ def get_fock(self, dm, h1e, **kwargs): """ return super().get_fock(dm, h1e, **kwargs) - def get_q_ij(self,q,mo_energy_w): + def get_q_ij(self, q, mo_energy_w): cell = self.with_df.cell kpts = self.kpts coords, weights = get_becke_grids(cell, level=5) @@ -707,24 +706,24 @@ def get_q_ij(self,q,mo_energy_w): psi_all = eval_ao(cell, coords, kpt=kpts[ki], deriv=1) psi = psi_all[0] psi_div = psi_all[1:4] - braket = lib.einsum( - "w,pw,dwq->dpq", weights, psi.T.conj(), psi_div - ) + braket = lib.einsum("w,pw,dwq->dpq", weights, psi.T.conj(), psi_div) - num_ao = -1.0j * lib.einsum( - "d,dpq->pq", q, braket - ) + num_ao = -1.0j * lib.einsum("d,dpq->pq", q, braket) num_mo = np.linalg.multi_dot( - (self.mo_coeff_w[ki][:, self.mo_occ_w[ki] > 0].T.conj(), - num_ao, - self.mo_coeff_w[ki][:, self.mo_occ_w[ki] == 0]) + ( + self.mo_coeff_w[ki][:, self.mo_occ_w[ki] > 0].T.conj(), + num_ao, + self.mo_coeff_w[ki][:, self.mo_occ_w[ki] == 0], + ) + ) + den = 1 / ( + mo_energy_w[ki][self.mo_occ_w[ki] == 0, None] + - mo_energy_w[ki][None, self.mo_occ_w[ki] > 0] ) - den = 1/(mo_energy_w[ki][self.mo_occ_w[ki] == 0, None] - mo_energy_w[ki][None,self.mo_occ_w[ki] > 0]) qij[ki] = (den.T * num_mo).flatten() return qij - def reciprocal_lattice(self): """ Return the reciprocal lattice vectors. diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index a4d7c4ff..030a7eab 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -37,6 +37,7 @@ class dTDA(MoldTDA): interaction, respectively. If `None`, use `gw.mo_occ` for both. Default value is `None`. """ + def __init__( self, gw, @@ -68,7 +69,6 @@ def build_dd_moments(self): kpts = self.kpts moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) - if self.fc: head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) @@ -79,7 +79,9 @@ def build_dd_moments(self): moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts if self.fc: - head[q, 0] += (np.sqrt(4. * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[q].conj() + head[q, 0] += (np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[ + q + ].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -92,7 +94,7 @@ def build_dd_moments(self): (self.mo_occ_w[kj], self.mo_occ_w[kb]), ) moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] - if self.fc and q==0: + if self.fc and q == 0: head[kj, i] += head[kj, i - 1] * d.ravel() tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) @@ -104,9 +106,9 @@ def build_dd_moments(self): tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) if q == 0 and self.fc: - tmp_head += lib.einsum("a,aP->P",head[kj, i - 1], - self.integrals.Lia[kb, kb].T.conj()) - + tmp_head += lib.einsum( + "a,aP->P", head[kj, i - 1], self.integrals.Lia[kb, kb].T.conj() + ) tmp = mpi_helper.allreduce(tmp) tmp *= 2.0 @@ -122,12 +124,10 @@ def build_dd_moments(self): moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) if q == 0 and self.fc: - head[kj, i] += lib.einsum("P,Pa->a", tmp_head, - self.integrals.Lia[kj, kj]) - + head[kj, i] += lib.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) if self.fc: - return {"moments":moments, "head":head} + return {"moments": moments, "head": head} else: return moments @@ -270,7 +270,7 @@ def build_se_moments(self, moments_dd): if self.fc: cell_vol = self.kpts.cell.vol total_vol = cell_vol * self.nkpts - q0 = (6*np.pi**2/total_vol)**(1/3) + q0 = (6 * np.pi**2 / total_vol) ** (1 / 3) norm_q_abs = np.linalg.norm(self.q_abs[0]) eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) @@ -282,14 +282,16 @@ def build_se_moments(self, moments_dd): for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) if self.fc: - eta_aux += np.dot(moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj()) - if q==0: - eta_head[kb, n] += (-(np.sqrt(4. * np.pi) / norm_q_abs) * - np.sum(moments_dd["head"][kb, n] * - self.qij[kb])) - eta_wings[kb, n] += ((np.sqrt(4. * np.pi) / norm_q_abs) * - lib.einsum("Pa,a->P", moments_dd["moments"][q, kb, n], - self.qij[kb])) + eta_aux += np.dot( + moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj() + ) + if q == 0: + eta_head[kb, n] += -(np.sqrt(4.0 * np.pi) / norm_q_abs) * np.sum( + moments_dd["head"][kb, n] * self.qij[kb] + ) + eta_wings[kb, n] += (np.sqrt(4.0 * np.pi) / norm_q_abs) * lib.einsum( + "Pa,a->P", moments_dd["moments"][q, kb, n], self.qij[kb] + ) else: eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) @@ -310,18 +312,17 @@ def build_se_moments(self, moments_dd): if self.fc and q == 0: original = eta[kp, q][x, n] - eta[kp, q][x, n] += (2/np.pi)*(q0)*eta_head[kp, n]*original + eta[kp, q][x, n] += (2 / np.pi) * (q0) * eta_head[kp, n] * original wing_tmp = lib.einsum("Pp,P->p", Lp, eta_wings[kp, n]) - wing_tmp = wing_tmp.real*2 - wing_tmp *= -(np.sqrt(cell_vol/ (4*(np.pi** 3))) * q0**2) + wing_tmp = wing_tmp.real * 2 + wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi**3))) * q0**2) eta[kp, q][x, n] += lib.einsum("p,pq->pq", wing_tmp, original) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) - return moments_occ, moments_vir def build_pert_term(self, qpt): From c7beb4fc40095a8f97f258e237836d38c532e0e8 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 17:59:02 +0100 Subject: [PATCH 13/33] More linting --- momentGW/pbc/gw.py | 45 +++++++------------------------------------- momentGW/pbc/ints.py | 27 -------------------------- momentGW/pbc/tda.py | 3 --- 3 files changed, 7 insertions(+), 68 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 92a7e401..7e4b2a1b 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -154,6 +154,13 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): Moments of the particle self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ + if self.fc: + fc = True + else: + fc = False + kwargs = dict( + fc=fc, + ) if self.polarizability.lower() == "dtda": tda = dTDA(self, nmom_max, integrals, **kwargs) @@ -219,44 +226,6 @@ def ao2mo(self, transform=True): return integrals - def build_se_moments(self, nmom_max, integrals, **kwargs): - """Build the moments of the self-energy. - - Parameters - ---------- - nmom_max : int - Maximum moment number to calculate. - integrals : KIntegrals - Density-fitted integrals. - - See functions in `momentGW.rpa` for `kwargs` options. - - Returns - ------- - se_moments_hole : numpy.ndarray - Moments of the hole self-energy at each k-point. If - `self.diagonal_se`, non-diagonal elements are set to zero. - se_moments_part : numpy.ndarray - Moments of the particle self-energy at each k-point. If - `self.diagonal_se`, non-diagonal elements are set to zero. - """ - if self.fc: - fc = True - else: - fc = False - kwargs = dict( - fc=fc, - ) - - if self.polarizability.lower() == "dtda": - tda = dTDA(self, nmom_max, integrals, **kwargs) - return tda.kernel() - elif self.polarizability.lower() == "thc-dtda": - tda = thc.dTDA(self, nmom_max, integrals, **kwargs) - return tda.kernel() - else: - raise NotImplementedError - def solve_dyson(self, se_moments_hole, se_moments_part, se_static, integrals=None): """Solve the Dyson equation due to a self-energy resulting from a list of hole and particle moments, along with a static diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 62b09ddd..58c31643 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -697,33 +697,6 @@ def get_fock(self, dm, h1e, **kwargs): """ return super().get_fock(dm, h1e, **kwargs) - def get_q_ij(self, q, mo_energy_w): - cell = self.with_df.cell - kpts = self.kpts - coords, weights = get_becke_grids(cell, level=5) - qij = np.zeros((len(kpts), self.nocc_w[0] * self.nvir_w[0]), dtype=complex) - for ki in kpts.loop(1): - psi_all = eval_ao(cell, coords, kpt=kpts[ki], deriv=1) - psi = psi_all[0] - psi_div = psi_all[1:4] - braket = lib.einsum("w,pw,dwq->dpq", weights, psi.T.conj(), psi_div) - - num_ao = -1.0j * lib.einsum("d,dpq->pq", q, braket) - num_mo = np.linalg.multi_dot( - ( - self.mo_coeff_w[ki][:, self.mo_occ_w[ki] > 0].T.conj(), - num_ao, - self.mo_coeff_w[ki][:, self.mo_occ_w[ki] == 0], - ) - ) - den = 1 / ( - mo_energy_w[ki][self.mo_occ_w[ki] == 0, None] - - mo_energy_w[ki][None, self.mo_occ_w[ki] > 0] - ) - qij[ki] = (den.T * num_mo).flatten() - - return qij - def reciprocal_lattice(self): """ Return the reciprocal lattice vectors. diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 030a7eab..40c41831 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -5,9 +5,7 @@ import numpy as np import scipy.special from pyscf import lib -from pyscf.agf2 import mpi_helper from pyscf.pbc import dft -from pyscf.pbc.gw.krgw_ac import get_qij from momentGW import logging, mpi_helper, util from momentGW.tda import dTDA as MoldTDA @@ -332,7 +330,6 @@ def build_pert_term(self, qpt): """ coords, weights = dft.gen_grid.get_becke_grids(self.gw.cell, level=5) - ngrid = len(coords) qij = np.zeros((self.nkpts,), dtype=object) for k in self.kpts.loop(1): From 26e873b0ff4210812aacf18723b385b6db1a0df5 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 18:00:58 +0100 Subject: [PATCH 14/33] Even more linting --- momentGW/pbc/ints.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 58c31643..ff1a6d6f 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -9,8 +9,6 @@ from pyscf import lib from pyscf.ao2mo import _ao2mo from pyscf.pbc import tools -from pyscf.pbc.dft.gen_grid import get_becke_grids -from pyscf.pbc.dft.numint import eval_ao from scipy.linalg import cholesky from momentGW import logging, mpi_helper, util From ad9cdeeaa1da1fcd08b335595de0347095db9783 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 18:40:26 +0100 Subject: [PATCH 15/33] Fix kwargs --- momentGW/pbc/gw.py | 6 ++---- momentGW/pbc/ints.py | 32 ++++++++++++++++++++++++++++++++ momentGW/pbc/tda.py | 4 ++-- 3 files changed, 36 insertions(+), 6 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 7e4b2a1b..e99b1464 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -158,12 +158,10 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): fc = True else: fc = False - kwargs = dict( - fc=fc, - ) + if self.polarizability.lower() == "dtda": - tda = dTDA(self, nmom_max, integrals, **kwargs) + tda = dTDA(self, nmom_max, integrals, fc, **kwargs) return tda.kernel() if self.polarizability.lower() == "drpa": rpa = dRPA(self, nmom_max, integrals, **kwargs) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index ff1a6d6f..47a22eea 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -417,6 +417,38 @@ def get_cderi_from_thc(self): self._blocks["Lia"] = Lia self._blocks["Lai"] = Lai + def update_coeffs(self, mo_coeff_g=None, mo_coeff_w=None, mo_occ_w=None): + """ + Update the MO coefficients in-place for the Green's function + and the screened Coulomb interaction. + + Parameters + ---------- + mo_coeff_g : numpy.ndarray, optional + Coefficients corresponding to the Green's function at each + k-point. Default value is `None`. + mo_coeff_w : numpy.ndarray, optional + Coefficients corresponding to the screened Coulomb + interaction at each k-point. Default value is `None`. + mo_occ_w : numpy.ndarray, optional + Occupations corresponding to the screened Coulomb + interaction at each k-point. Default value is `None`. + + Notes + ----- + If `mo_coeff_g` is `None`, the Green's function is assumed to + remain in the basis in which it was originally defined, and + vice-versa for `mo_coeff_w` and `mo_occ_w`. At least one of + `mo_coeff_g` and `mo_coeff_w` must be provided. + """ + return super().update_coeffs( + mo_coeff_g=mo_coeff_g, + mo_coeff_w=mo_coeff_w, + mo_occ_w=mo_occ_w, + ) + + @logging.with_timer("J matrix") + @logging.with_status("Building J matrix") def get_j(self, dm, basis="mo", other=None): """Build the J matrix. diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 40c41831..11cccf43 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -41,11 +41,11 @@ def __init__( gw, nmom_max, integrals, + fc=False, mo_energy=None, mo_occ=None, - fc=False, ): - super().__init__(gw, nmom_max, integrals, mo_energy, mo_occ) + super().__init__(gw, nmom_max, integrals, mo_energy=mo_energy, mo_occ=mo_occ) self.fc = fc if self.fc: q = np.array([1e-3, 0, 0]).reshape(1, 3) From 43370f5cde85e1ce3d2c57726cc48f2e424775e2 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 4 Apr 2024 18:42:22 +0100 Subject: [PATCH 16/33] Linting --- momentGW/pbc/gw.py | 1 - 1 file changed, 1 deletion(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index e99b1464..02a151d5 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -159,7 +159,6 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): else: fc = False - if self.polarizability.lower() == "dtda": tda = dTDA(self, nmom_max, integrals, fc, **kwargs) return tda.kernel() From d0acdee55e0ef42735295b1133803f8da831e772 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Fri, 5 Apr 2024 17:43:40 +0100 Subject: [PATCH 17/33] einsum _contraction order correction --- momentGW/util.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/momentGW/util.py b/momentGW/util.py index be7ba390..07232b93 100644 --- a/momentGW/util.py +++ b/momentGW/util.py @@ -473,7 +473,12 @@ def _fallback(): # Reshape and transpose the output ct = ct.reshape(shape_ct, order="A") - c = ct.transpose(order_ct) + if ct.flags.f_contiguous: + c = np.asfortranarray(ct.transpose(order_ct)) + elif ct.flags.c_contiguous: + c = np.ascontiguousarray(ct.transpose(order_ct)) + else: + c = ct.transpose(order_ct) return c From fe621830edbaf6bf8ab88f6140d965a09331b005 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 8 Apr 2024 17:24:48 +0100 Subject: [PATCH 18/33] Additional MPI changes --- momentGW/pbc/tda.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 11cccf43..7154a44b 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -76,10 +76,10 @@ def build_dd_moments(self): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts - if self.fc: - head[q, 0] += (np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0])) * self.qij[ - q - ].conj() + if q == 0 and self.fc: + head[kj, 0] += ( + np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0]) + ) * self.qij[kj].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -105,7 +105,7 @@ def build_dd_moments(self): if q == 0 and self.fc: tmp_head += lib.einsum( - "a,aP->P", head[kj, i - 1], self.integrals.Lia[kb, kb].T.conj() + "a,aP->P", head[ki, i - 1], self.integrals.Lia[ki, ki].T.conj() ) tmp = mpi_helper.allreduce(tmp) From 62fd849031177d841eda28cfcb333ea23dbc94e3 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 10 Apr 2024 17:45:50 +0100 Subject: [PATCH 19/33] Change to util einsums --- momentGW/pbc/tda.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 7154a44b..c2a24e77 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -104,7 +104,7 @@ def build_dd_moments(self): tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) if q == 0 and self.fc: - tmp_head += lib.einsum( + tmp_head += util.einsum( "a,aP->P", head[ki, i - 1], self.integrals.Lia[ki, ki].T.conj() ) @@ -122,7 +122,7 @@ def build_dd_moments(self): moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) if q == 0 and self.fc: - head[kj, i] += lib.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) + head[kj, i] += util.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) if self.fc: return {"moments": moments, "head": head} @@ -287,7 +287,7 @@ def build_se_moments(self, moments_dd): eta_head[kb, n] += -(np.sqrt(4.0 * np.pi) / norm_q_abs) * np.sum( moments_dd["head"][kb, n] * self.qij[kb] ) - eta_wings[kb, n] += (np.sqrt(4.0 * np.pi) / norm_q_abs) * lib.einsum( + eta_wings[kb, n] += (np.sqrt(4.0 * np.pi) / norm_q_abs) * util.einsum( "Pa,a->P", moments_dd["moments"][q, kb, n], self.qij[kb] ) else: @@ -312,11 +312,11 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += (2 / np.pi) * (q0) * eta_head[kp, n] * original - wing_tmp = lib.einsum("Pp,P->p", Lp, eta_wings[kp, n]) + wing_tmp = util.einsum("Pp,P->p", Lp, eta_wings[kp, n]) wing_tmp = wing_tmp.real * 2 wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi**3))) * q0**2) - eta[kp, q][x, n] += lib.einsum("p,pq->pq", wing_tmp, original) + eta[kp, q][x, n] += util.einsum("p,pq->pq", wing_tmp, original) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) @@ -336,9 +336,9 @@ def build_pert_term(self, qpt): ao_p = dft.numint.eval_ao(self.gw.cell, coords, kpt=self.kpts[k], deriv=1) ao, ao_grad = ao_p[0], ao_p[1:4] - ao_ao_grad = lib.einsum("g,gm,xgn->xmn", weights, ao.conj(), ao_grad) - q_ao_ao_grad = lib.einsum("x,xmn->mn", qpt, ao_ao_grad) * -1.0j - q_mo_mo_grad = lib.einsum( + ao_ao_grad = util.einsum("g,gm,xgn->xmn", weights, ao.conj(), ao_grad) + q_ao_ao_grad = util.einsum("x,xmn->mn", qpt, ao_ao_grad) * -1.0j + q_mo_mo_grad = util.einsum( "mn,mi,na->ia", q_ao_ao_grad, self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] > 0].conj(), From cdc27568f7d10276ad14aca1fbd7412b811e0c85 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 17 Apr 2024 18:37:14 +0100 Subject: [PATCH 20/33] HW test for kTDA --- tests/test_kgw.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/test_kgw.py b/tests/test_kgw.py index 1b0080d6..baa72b34 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -168,6 +168,19 @@ def test_drpa_vs_supercell_compression(self): self._test_vs_supercell(gw, kgw, full=False, tol=1e-5) + def test_dtda_HW_regression(self): + nmom_max = 5 + kgw = KGW(self.mf) + kgw.polarizability = "dtda" + kgw.fc = True + conv, gf, se, _ = kgw.kernel(nmom_max) + gf_occ = gf[0].occupied().physical(weight=1e-1) + gf_vir = gf[0].virtual().physical(weight=1e-1) + self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7516926143, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925912088, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900808711, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839604877, 6) + if __name__ == "__main__": print("Running tests for KGW") From f16b1d9efefaca74a05913382d6348652d172768 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 17 Apr 2024 18:55:51 +0100 Subject: [PATCH 21/33] Update test --- tests/test_kgw.py | 42 ++++++++++++++---------------------------- 1 file changed, 14 insertions(+), 28 deletions(-) diff --git a/tests/test_kgw.py b/tests/test_kgw.py index 7499414d..dbbfd814 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -108,6 +108,20 @@ def test_drpa_vs_supercell(self): self._test_vs_supercell(gw, kgw, full=True) + def test_dtda_HW_regression(self): + nmom_max = 5 + kgw = KGW(self.mf) + kgw.polarizability = "dtda" + kgw.fc = True + conv, gf, se, _ = kgw.kernel(nmom_max) + gf_occ = gf[0].occupied().physical(weight=1e-1) + gf_vir = gf[0].virtual().physical(weight=1e-1) + self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7516810595, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925912088, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900808711, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839604877, 6) + + def test_dtda_vs_supercell_fock_loop(self): nmom_max = 5 @@ -168,34 +182,6 @@ def test_drpa_vs_supercell_compression(self): self._test_vs_supercell(gw, kgw, full=False, tol=1e-5) - def test_dtda_vs_supercell_frozen(self): - nmom_max = 3 - - kgw = KGW(self.mf) - kgw.polarizability = "dtda" - kgw.frozen = [0] - kgw.kernel(nmom_max) - - gw = GW(self.smf) - gw.__dict__.update({opt: getattr(kgw, opt) for opt in kgw._opts}) - gw.frozen = list(range(len(self.kpts))) - gw.kernel(nmom_max) - - self._test_vs_supercell(gw, kgw, full=True) - - def test_dtda_HW_regression(self): - nmom_max = 5 - kgw = KGW(self.mf) - kgw.polarizability = "dtda" - kgw.fc = True - conv, gf, se, _ = kgw.kernel(nmom_max) - gf_occ = gf[0].occupied().physical(weight=1e-1) - gf_vir = gf[0].virtual().physical(weight=1e-1) - self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7516926143, 6) - self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925912088, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900808711, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839604877, 6) - if __name__ == "__main__": print("Running tests for KGW") unittest.main() From 1fe4a2c7c2b6387a502996cd86814755a3972256 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 18 Apr 2024 15:42:10 +0100 Subject: [PATCH 22/33] Frozen static field fix --- momentGW/pbc/gw.py | 7 ++++++- tests/test_kgw.py | 6 +++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 2e4e59c3..977448b6 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -91,8 +91,12 @@ def build_se_static(self, integrals): non-diagonal elements are set to zero. """ + mask = self.active + dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) + if getattr(self._scf, "xc", "hf") == "hf": - se_static = np.zeros_like(self._scf.make_rdm1(mo_coeff=self.mo_coeff)) + se_static = np.zeros_like(dm) + se_static = se_static[..., mask, :][..., :, mask] if self.fc: with util.SilentSCF(self._scf): vmf = self._scf.get_j() - self._scf.get_veff() @@ -121,6 +125,7 @@ def build_se_static(self, integrals): for k in range(len(self.kpts)): vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) se_static = vmf - vk * 0.5 + se_static = se_static[..., mask, :][..., :, mask] se_static = util.einsum( "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff diff --git a/tests/test_kgw.py b/tests/test_kgw.py index dbbfd814..533cad0c 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -117,9 +117,9 @@ def test_dtda_HW_regression(self): gf_occ = gf[0].occupied().physical(weight=1e-1) gf_vir = gf[0].virtual().physical(weight=1e-1) self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7516810595, 6) - self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925912088, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900808711, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839604877, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925838943, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900606812, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839090190, 6) def test_dtda_vs_supercell_fock_loop(self): From 79d86961ea9062e75cc3cec7450685aef6734f1b Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 1 May 2024 18:54:16 +0100 Subject: [PATCH 23/33] Changes pre discussion in PR --- momentGW/pbc/gw.py | 73 ++++++---------- momentGW/pbc/ints.py | 5 +- momentGW/pbc/tda.py | 204 +++++++++++++++++++++++++++---------------- tests/test_kgw.py | 15 ++++ 4 files changed, 173 insertions(+), 124 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 977448b6..b0a5e4c6 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -81,62 +81,45 @@ def build_se_static(self, integrals): Parameters ---------- - integrals : Integrals + integrals : KIntegrals Integrals object. Returns ------- se_static : numpy.ndarray - Static part of the self-energy. If `self.diagonal_se`, - non-diagonal elements are set to zero. + Static part of the self-energy at each k-point. If + `self.diagonal_se`, non-diagonal elements are set to zero. """ + if self.fc: + # Get intermediates + mask = self.active + dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) - mask = self.active - dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) - - if getattr(self._scf, "xc", "hf") == "hf": - se_static = np.zeros_like(dm) - se_static = se_static[..., mask, :][..., :, mask] - if self.fc: - with util.SilentSCF(self._scf): - vmf = self._scf.get_j() - self._scf.get_veff() - dm = self._scf.make_rdm1(mo_coeff=self.mo_coeff) - vk = integrals.get_k(dm, basis="ao") - - s = self.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts) - madelung = tools.pbc.madelung(self.cell, self.kpts) - for k in range(len(self.kpts)): - vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) - - se_static = vmf - vk * 0.5 - se_static = util.einsum( - "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff - ) - - else: + # Get the contribution from the exchange-correlation potential with util.SilentSCF(self._scf): - vmf = self._scf.get_j() - self._scf.get_veff() - dm = self._scf.make_rdm1(mo_coeff=self.mo_coeff) - vk = integrals.get_k(dm, basis="ao") - - if self.fc: - s = self.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts) - madelung = tools.pbc.madelung(self.cell, self.kpts) - for k in range(len(self.kpts)): - vk[k] += madelung * reduce(np.dot, (s[k], dm[k], s[k])) - se_static = vmf - vk * 0.5 - se_static = se_static[..., mask, :][..., :, mask] + veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, + mask] + vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] + vhf = integrals.get_veff(dm, j=vj, basis="ao", ewald=self.fc) + se_static = vhf - veff se_static = util.einsum( - "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff + "...pq,...pi,...qj->...ij", se_static, + np.conj(self.mo_coeff), self.mo_coeff ) + # If diagonal approximation, set non-diagonal elements to zero + if self.diagonal_se: + se_static = util.einsum("...pq,pq->...pq", se_static, + np.eye(se_static.shape[-1])) - if self.diagonal_se: - se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) + # Add the Fock matrix contribution + se_static += util.einsum("...p,...pq->...pq", self.mo_energy, + np.eye(se_static.shape[-1])) - se_static += util.einsum("...p,...pq->...pq", self.mo_energy, np.eye(se_static.shape[-1])) + return se_static - return se_static + else: + return super().build_se_static(integrals) def build_se_moments(self, nmom_max, integrals, **kwargs): """Build the moments of the self-energy. @@ -159,13 +142,9 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): Moments of the particle self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - if self.fc: - fc = True - else: - fc = False if self.polarizability.lower() == "dtda": - tda = dTDA(self, nmom_max, integrals, fc, **kwargs) + tda = dTDA(self, nmom_max, integrals, fc=self.fc, **kwargs) return tda.kernel() if self.polarizability.lower() == "drpa": rpa = dRPA(self, nmom_max, integrals, **kwargs) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index e198a29b..a0dae9a5 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -667,8 +667,11 @@ def get_ewald(self, dm, basis="mo"): else: ovlp = self.with_df.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts._kpts) + madelung = tools.pbc.madelung(self.with_df.cell, self.kpts._kpts) + # Initialise the Ewald matrix - ew = util.einsum("kpq,kpi,kqj->kij", dm, ovlp.conj(), ovlp) + ovlp = np.asarray(ovlp) + ew = madelung*util.einsum("kpq,kpi,kqj->kij", dm, ovlp.conj(), ovlp) return ew diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index c2a24e77..475f69fe 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -41,16 +41,14 @@ def __init__( gw, nmom_max, integrals, - fc=False, mo_energy=None, mo_occ=None, + fc=False, ): super().__init__(gw, nmom_max, integrals, mo_energy=mo_energy, mo_occ=mo_occ) self.fc = fc if self.fc: - q = np.array([1e-3, 0, 0]).reshape(1, 3) - self.q_abs = self.kpts.cell.get_abs_kpts(q) - self.qij = self.build_pert_term(self.q_abs[0]) + self.pw_hw = self.build_pert_term(self.q_abs[0]) @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -67,23 +65,16 @@ def build_dd_moments(self): kpts = self.kpts moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) - if self.fc: - head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - # Get the zeroth order moment for q in kpts.loop(1): for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, 0] += self.integrals.Lia[kj, kb] / self.nkpts - if q == 0 and self.fc: - head[kj, 0] += ( - np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0]) - ) * self.qij[kj].conj() - # Get the higher order moments for i in range(1, self.nmom_max + 1): for q in kpts.loop(1): + tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) @@ -92,42 +83,116 @@ def build_dd_moments(self): (self.mo_occ_w[kj], self.mo_occ_w[kb]), ) moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None] - if self.fc and q == 0: - head[kj, i] += head[kj, i - 1] * d.ravel() - - tmp = np.zeros((self.naux[q], self.naux[q]), dtype=complex) - tmp_head = np.zeros((self.naux[q]), dtype=complex) - - for ki in kpts.loop(1, mpi=True): - ka = kpts.member(kpts.wrap_around(kpts[q] + kpts[ki])) - - tmp += np.dot(moments[q, ka, i - 1], self.integrals.Lia[ki, ka].T.conj()) - if q == 0 and self.fc: - tmp_head += util.einsum( - "a,aP->P", head[ki, i - 1], self.integrals.Lia[ki, ki].T.conj() - ) + tmp += np.dot(moments[q, kb, i - 1], self.integrals.Lia[kj, kb].T.conj()) tmp = mpi_helper.allreduce(tmp) tmp *= 2.0 tmp /= self.nkpts - tmp_head = mpi_helper.allreduce(tmp_head) - tmp_head *= 2.0 - tmp_head /= self.nkpts - for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, i] += np.dot(tmp, self.integrals.Lai[kj, kb].conj()) - if q == 0 and self.fc: - head[kj, i] += util.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) + return moments + + def build_head_dd_moments(self): + """Build the head for the moments of the density-density response. + + Returns + ------- + head_moments : numpy.ndarray + Head moments of the density-density response at each k-point. + """ + kpts = self.kpts + head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + HW_const = np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0]) + + # Get the zeroth order moment + for kj in kpts.loop(1, mpi=True): + head[kj, 0] += HW_const* self.pw_hw[kj].conj() + + # Get the higher order moments + for i in range(1, self.nmom_max + 1): + tmp_head = np.zeros((self.naux[0]), dtype=complex) + for kj in kpts.loop(1, mpi=True): + d = util.build_1h1p_energies( + (self.mo_energy_w[kj], self.mo_energy_w[kj]), + (self.mo_occ_w[kj], self.mo_occ_w[kj]), + ) + head[kj, i] += head[kj, i - 1] * d.ravel() + + tmp_head += util.einsum( + "a,aP->P", head[kj, i - 1], self.integrals.Lia[kj, kj].T.conj() + ) + + tmp_head = mpi_helper.allreduce(tmp_head) + tmp_head *= 2.0 + tmp_head /= self.nkpts + + for kj in kpts.loop(1, mpi=True): + head[kj, i] += util.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) + + return head + + def _add_eta_fc_correction(self,moments_dd, eta): + """Correct the moments of the self-energy using the Head + and wings corrections. + + Parameters + ---------- + moments_dd : numpy.ndarray + Moments of the density-density response at each k-point. + + eta : numpy.ndarray + Moments of the density-density response partly transformed + into moments of the screened Coulomb interaction at each + k-point. + + Returns + ------- + eta : numpy.ndarray + Corrected moments of the density-density response partly + transformed into moments of the screened Coulomb interaction + at each k-point. + + """ + head = self.build_head_dd_moments() + + kpts = self.kpts + + cell_vol = self.kpts.cell.vol + total_vol = self.kpts.cell.vol * self.nkpts + q0 = (6 * np.pi ** 2 / total_vol) ** (1 / 3) + norm_q_abs = np.linalg.norm(self.q_abs[0]) + HW_const = np.sqrt(4.0 * np.pi) / norm_q_abs + + eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + + for n in range(self.nmom_max + 1): + for kj in kpts.loop(1, mpi=True): + eta_head[kj, n] += -HW_const * np.sum(head[kj, n] * self.pw_hw[kj]) + eta_wings[kj, n] += HW_const * util.einsum( + "Pa,a->P", moments_dd[0, kj, n], self.pw_hw[kj] + ) + for kx in kpts.loop(1, mpi=True): + for x in range(self.mo_energy_g[kx].size): + Lp = self.integrals.Lpx[kx, kx][:, :, x] + + original = eta[kx, 0][x, n] + + eta[kx, 0][x, n] += (2 / np.pi) * (q0) * eta_head[kx, n] * original + + wing_tmp = util.einsum("Pp,P->p", Lp, eta_wings[kx, n]) + wing_tmp = wing_tmp.real * 2 + wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi ** 3))) * q0 ** 2) + + eta[kx, 0][x, n] += util.einsum("p,pq->pq", wing_tmp, original) + return eta + - if self.fc: - return {"moments": moments, "head": head} - else: - return moments def kernel(self, exact=False): """ @@ -265,33 +330,13 @@ def build_se_moments(self, moments_dd): eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) eta = np.zeros((self.nkpts, self.nkpts), dtype=object) - if self.fc: - cell_vol = self.kpts.cell.vol - total_vol = cell_vol * self.nkpts - q0 = (6 * np.pi**2 / total_vol) ** (1 / 3) - norm_q_abs = np.linalg.norm(self.q_abs[0]) - eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - # Get the moments in (aux|aux) and rotate to (mo|mo) for n in range(self.nmom_max + 1): for q in kpts.loop(1): eta_aux = 0 for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - if self.fc: - eta_aux += np.dot( - moments_dd["moments"][q, kb, n], self.integrals.Lia[kj, kb].T.conj() - ) - if q == 0: - eta_head[kb, n] += -(np.sqrt(4.0 * np.pi) / norm_q_abs) * np.sum( - moments_dd["head"][kb, n] * self.qij[kb] - ) - eta_wings[kb, n] += (np.sqrt(4.0 * np.pi) / norm_q_abs) * util.einsum( - "Pa,a->P", moments_dd["moments"][q, kb, n], self.qij[kb] - ) - else: - eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) + eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) eta_aux = mpi_helper.allreduce(eta_aux) eta_aux *= 2.0 @@ -307,16 +352,9 @@ def build_se_moments(self, moments_dd): Lp = self.integrals.Lpx[kp, kx][:, :, x] subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) - if self.fc and q == 0: - original = eta[kp, q][x, n] - - eta[kp, q][x, n] += (2 / np.pi) * (q0) * eta_head[kp, n] * original - wing_tmp = util.einsum("Pp,P->p", Lp, eta_wings[kp, n]) - wing_tmp = wing_tmp.real * 2 - wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi**3))) * q0**2) - - eta[kp, q][x, n] += util.einsum("p,pq->pq", wing_tmp, original) + if self.fc: + eta = self._add_eta_fc_correction(moments_dd,eta) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) @@ -325,13 +363,23 @@ def build_se_moments(self, moments_dd): def build_pert_term(self, qpt): """ - Compute 1/sqrt(Ω) * ⟨Ψ_{ik}|e^{iqr}|Ψ_{ak-q}⟩ at q-point index - q using perturbation theory. + Build the charge-density density matrix at q-point index qpt + using perturbation theory. + + Parameters + ---------- + qpt : numpy.ndarray + q-point index representing the limit of our plane waves. + + Returns + ------- + pw_hw : numpy.ndarray + Charge-density density matrix in the long wavelength limit. """ coords, weights = dft.gen_grid.get_becke_grids(self.gw.cell, level=5) - qij = np.zeros((self.nkpts,), dtype=object) + pw_hw = np.zeros((self.nkpts,), dtype=object) for k in self.kpts.loop(1): ao_p = dft.numint.eval_ao(self.gw.cell, coords, kpt=self.kpts[k], deriv=1) ao, ao_grad = ao_p[0], ao_p[1:4] @@ -345,17 +393,16 @@ def build_pert_term(self, qpt): self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] == 0], ) - d = lib.direct_sum( - "a-i->ia", - self.mo_energy_w[k][self.mo_occ_w[k] == 0], - self.mo_energy_w[k][self.mo_occ_w[k] > 0], + d = util.build_1h1p_energies( + (self.mo_energy_w[k], self.mo_energy_w[k]), + (self.mo_occ_w[k], self.mo_occ_w[k]), ) + dens = q_mo_mo_grad / d - qij[k] = dens.flatten() / np.sqrt(self.gw.cell.vol) + pw_hw[k] = dens.flatten() / np.sqrt(self.gw.cell.vol) - return qij + return pw_hw - build_dd_moments_exact = build_dd_moments @property def naux(self): @@ -379,3 +426,8 @@ def kpts(self): def nkpts(self): """Get the number of k-points.""" return self.gw.nkpts + + @property + def q_abs(self): + """Get the absolute value of a region near the origin (q=0).""" + return self.kpts.cell.get_abs_kpts(np.array([1e-3, 0, 0]).reshape(1, 3)) diff --git a/tests/test_kgw.py b/tests/test_kgw.py index 533cad0c..0e4c710b 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -108,6 +108,21 @@ def test_drpa_vs_supercell(self): self._test_vs_supercell(gw, kgw, full=True) + def test_dtda_vs_supercell_frozen(self): + nmom_max = 3 + + kgw = KGW(self.mf) + kgw.polarizability = "dtda" + kgw.frozen = [0] + kgw.kernel(nmom_max) + + gw = GW(self.smf) + gw.__dict__.update({opt: getattr(kgw, opt) for opt in kgw._opts}) + gw.frozen = list(range(len(self.kpts))) + gw.kernel(nmom_max) + + self._test_vs_supercell(gw, kgw, full=True) + def test_dtda_HW_regression(self): nmom_max = 5 kgw = KGW(self.mf) From f00bbfc1ef4ab2603d2e91fb04a159f837cf5454 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Wed, 1 May 2024 19:04:41 +0100 Subject: [PATCH 24/33] Linting --- momentGW/pbc/gw.py | 14 ++++++-------- momentGW/pbc/tda.py | 13 +++++-------- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index b0a5e4c6..d98d744a 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -97,24 +97,22 @@ def build_se_static(self, integrals): # Get the contribution from the exchange-correlation potential with util.SilentSCF(self._scf): - veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, - mask] + veff = self._scf.get_veff(None, dm)[..., mask, :][..., :,mask] vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] vhf = integrals.get_veff(dm, j=vj, basis="ao", ewald=self.fc) se_static = vhf - veff se_static = util.einsum( - "...pq,...pi,...qj->...ij", se_static, - np.conj(self.mo_coeff), self.mo_coeff + "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff ) # If diagonal approximation, set non-diagonal elements to zero if self.diagonal_se: - se_static = util.einsum("...pq,pq->...pq", se_static, - np.eye(se_static.shape[-1])) + se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) # Add the Fock matrix contribution - se_static += util.einsum("...p,...pq->...pq", self.mo_energy, - np.eye(se_static.shape[-1])) + se_static += util.einsum( + "...p,...pq->...pq", self.mo_energy, np.eye(se_static.shape[-1]) + ) return se_static diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 475f69fe..0a1be4d6 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -111,7 +111,7 @@ def build_head_dd_moments(self): # Get the zeroth order moment for kj in kpts.loop(1, mpi=True): - head[kj, 0] += HW_const* self.pw_hw[kj].conj() + head[kj, 0] += HW_const * self.pw_hw[kj].conj() # Get the higher order moments for i in range(1, self.nmom_max + 1): @@ -136,7 +136,7 @@ def build_head_dd_moments(self): return head - def _add_eta_fc_correction(self,moments_dd, eta): + def _add_eta_fc_correction(self, moments_dd, eta): """Correct the moments of the self-energy using the Head and wings corrections. @@ -164,7 +164,7 @@ def _add_eta_fc_correction(self,moments_dd, eta): cell_vol = self.kpts.cell.vol total_vol = self.kpts.cell.vol * self.nkpts - q0 = (6 * np.pi ** 2 / total_vol) ** (1 / 3) + q0 = (6 * np.pi**2 / total_vol) ** (1 / 3) norm_q_abs = np.linalg.norm(self.q_abs[0]) HW_const = np.sqrt(4.0 * np.pi) / norm_q_abs @@ -187,13 +187,11 @@ def _add_eta_fc_correction(self,moments_dd, eta): wing_tmp = util.einsum("Pp,P->p", Lp, eta_wings[kx, n]) wing_tmp = wing_tmp.real * 2 - wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi ** 3))) * q0 ** 2) + wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi**3))) * q0**2) eta[kx, 0][x, n] += util.einsum("p,pq->pq", wing_tmp, original) return eta - - def kernel(self, exact=False): """ Run the polarizability calculation to compute moments of the @@ -354,7 +352,7 @@ def build_se_moments(self, moments_dd): eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) if self.fc: - eta = self._add_eta_fc_correction(moments_dd,eta) + eta = self._add_eta_fc_correction(moments_dd, eta) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) @@ -403,7 +401,6 @@ def build_pert_term(self, qpt): return pw_hw - @property def naux(self): """Number of auxiliaries.""" From cbe49f82782aba7358a0bb24c5dcf1a739e988b0 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 2 May 2024 02:46:41 +0100 Subject: [PATCH 25/33] Further linting --- momentGW/pbc/gw.py | 2 +- momentGW/pbc/ints.py | 2 +- tests/test_kgw.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index d98d744a..40443377 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -97,7 +97,7 @@ def build_se_static(self, integrals): # Get the contribution from the exchange-correlation potential with util.SilentSCF(self._scf): - veff = self._scf.get_veff(None, dm)[..., mask, :][..., :,mask] + veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask] vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] vhf = integrals.get_veff(dm, j=vj, basis="ao", ewald=self.fc) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index a0dae9a5..eb77ad1c 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -671,7 +671,7 @@ def get_ewald(self, dm, basis="mo"): # Initialise the Ewald matrix ovlp = np.asarray(ovlp) - ew = madelung*util.einsum("kpq,kpi,kqj->kij", dm, ovlp.conj(), ovlp) + ew = madelung * util.einsum("kpq,kpi,kqj->kij", dm, ovlp.conj(), ovlp) return ew diff --git a/tests/test_kgw.py b/tests/test_kgw.py index 0e4c710b..72b70656 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -136,7 +136,6 @@ def test_dtda_HW_regression(self): self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900606812, 6) self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839090190, 6) - def test_dtda_vs_supercell_fock_loop(self): nmom_max = 5 @@ -197,6 +196,7 @@ def test_drpa_vs_supercell_compression(self): self._test_vs_supercell(gw, kgw, full=False, tol=1e-5) + if __name__ == "__main__": print("Running tests for KGW") unittest.main() From 5d845b4f1c1f6cbb9025628f3c5ea442d95cd37e Mon Sep 17 00:00:00 2001 From: mkakcl Date: Thu, 2 May 2024 02:49:47 +0100 Subject: [PATCH 26/33] More linting --- momentGW/pbc/gw.py | 3 --- momentGW/pbc/tda.py | 1 - 2 files changed, 4 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 40443377..ad531463 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -3,11 +3,8 @@ periodic systems. """ -from functools import reduce - import numpy as np from dyson import MBLSE, Lehmann, MixedMBLSE -from pyscf.pbc import tools from momentGW import energy, logging, util from momentGW.gw import GW diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 0a1be4d6..7b03d2de 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -4,7 +4,6 @@ import numpy as np import scipy.special -from pyscf import lib from pyscf.pbc import dft from momentGW import logging, mpi_helper, util From e5e2b1a0dca8bb5578a4ccd6081eec6288e8fa23 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 28 May 2024 17:10:19 +0100 Subject: [PATCH 27/33] Support inheritance for static SE with finite size corrections --- momentGW/gw.py | 38 ++++++++++++++++++++++-- momentGW/pbc/gw.py | 67 ++++++++++++++++++++++++------------------ momentGW/pbc/uhf/gw.py | 8 +++-- momentGW/uhf/gw.py | 6 ++-- 4 files changed, 83 insertions(+), 36 deletions(-) diff --git a/momentGW/gw.py b/momentGW/gw.py index e1d715fc..e4ebffbc 100644 --- a/momentGW/gw.py +++ b/momentGW/gw.py @@ -144,9 +144,37 @@ def name(self): polarizability = self.polarizability.upper().replace("DTDA", "dTDA").replace("DRPA", "dRPA") return f"{polarizability}-G0W0" + def get_veff(self, integrals, dm=None, **kwargs): + """Get the effective potential. + + Parameters + ---------- + integrals : Integrals + Integrals object. + dm : numpy.ndarray, optional + Density matrix. If `None`, determine using `self.make_rdm1`. + Default value is `None`. + **kwargs : dict, optional + Additional keyword arguments passed to the integrals object. + + Returns + ------- + veff : numpy.ndarray + Effective potential. + """ + + # Get the density matrix + if dm is None: + dm = self.make_rdm1() + + # Get the effective potential + veff = integrals.get_veff(dm, **kwargs) + + return veff + @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, force_build=False): """ Build the static part of the self-energy, including the Fock matrix. @@ -155,6 +183,10 @@ def build_se_static(self, integrals): ---------- integrals : Integrals Integrals object. + force_build : bool, optional + If `True`, force the static self-energy to be built, even if + the underlying mean-field object is Hartree--Fock. Default + value is `False`. Returns ------- @@ -168,7 +200,7 @@ def build_se_static(self, integrals): dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) # Get the contribution from the exchange-correlation potential - if getattr(self._scf, "xc", "hf") == "hf": + if getattr(self._scf, "xc", "hf") == "hf" and not force_build: se_static = np.zeros_like(dm) se_static = se_static[..., mask, :][..., :, mask] else: @@ -176,7 +208,7 @@ def build_se_static(self, integrals): veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask] vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] - vhf = integrals.get_veff(dm, j=vj, basis="ao") + vhf = self.get_veff(integrals, dm=dm, j=vj, basis="ao") se_static = vhf - veff se_static = util.einsum( "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index ad531463..58581e25 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -69,9 +69,41 @@ def name(self): polarizability = self.polarizability.upper().replace("DTDA", "dTDA").replace("DRPA", "dRPA") return f"{polarizability}-KG0W0" + def get_veff(self, integrals, dm=None, **kwargs): + """Get the effective potential. + + Parameters + ---------- + integrals : KIntegrals + Integrals object. + dm : numpy.ndarray, optional + Density matrix at each k-point. If `None`, determine using + `self.make_rdm1`. Default value is `None`. + **kwargs : dict, optional + Additional keyword arguments passed to the integrals object. + + Returns + ------- + veff : numpy.ndarray + Effective potential at each k-point. + """ + + # Get the density matrix + if dm is None: + dm = self.make_rdm1() + + # Set the default options + if "ewald" not in kwargs: + kwargs = {**kwargs, "ewald": self.fc} + + # Get the effective potential + veff = integrals.get_veff(dm, **kwargs) + + return veff + @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, **kwargs): """ Build the static part of the self-energy, including the Fock matrix. @@ -80,6 +112,8 @@ def build_se_static(self, integrals): ---------- integrals : KIntegrals Integrals object. + **kwargs : dict, optional + Additional keyword arguments. Returns ------- @@ -87,34 +121,9 @@ def build_se_static(self, integrals): Static part of the self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - if self.fc: - # Get intermediates - mask = self.active - dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) - - # Get the contribution from the exchange-correlation potential - with util.SilentSCF(self._scf): - veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask] - vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] - - vhf = integrals.get_veff(dm, j=vj, basis="ao", ewald=self.fc) - se_static = vhf - veff - se_static = util.einsum( - "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff - ) - # If diagonal approximation, set non-diagonal elements to zero - if self.diagonal_se: - se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) - - # Add the Fock matrix contribution - se_static += util.einsum( - "...p,...pq->...pq", self.mo_energy, np.eye(se_static.shape[-1]) - ) - - return se_static - - else: - return super().build_se_static(integrals) + if "force_build" not in kwargs: + kwargs = {**kwargs, "force_build": self.fc} + return super().build_se_static(integrals, **kwargs) def build_se_moments(self, nmom_max, integrals, **kwargs): """Build the moments of the self-energy. diff --git a/momentGW/pbc/uhf/gw.py b/momentGW/pbc/uhf/gw.py index 646afde9..1f03b5d9 100644 --- a/momentGW/pbc/uhf/gw.py +++ b/momentGW/pbc/uhf/gw.py @@ -71,7 +71,7 @@ def name(self): @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, **kwargs): """ Build the static part of the self-energy, including the Fock matrix. @@ -80,6 +80,8 @@ def build_se_static(self, integrals): ---------- integrals : KUIntegrals Integrals object. + **kwargs : dict, optional + Additional keyword arguments. Returns ------- @@ -88,7 +90,9 @@ def build_se_static(self, integrals): channel. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - return super().build_se_static(integrals) + if "force_build" not in kwargs: + kwargs = {**kwargs, "force_build": self.fc} + return super().build_se_static(integrals, **kwargs) @logging.with_timer("Integral construction") @logging.with_status("Constructing integrals") diff --git a/momentGW/uhf/gw.py b/momentGW/uhf/gw.py index 48d1d979..5a4f83b3 100644 --- a/momentGW/uhf/gw.py +++ b/momentGW/uhf/gw.py @@ -75,7 +75,7 @@ def name(self): @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, **kwargs): """ Build the static part of the self-energy, including the Fock matrix. @@ -84,6 +84,8 @@ def build_se_static(self, integrals): ---------- integrals : UIntegrals Integrals object. + **kwargs : dict, optional + Additional keyword arguments. Returns ------- @@ -91,7 +93,7 @@ def build_se_static(self, integrals): Static part of the self-energy for each spin channel. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - return super().build_se_static(integrals) + return super().build_se_static(integrals, **kwargs) def build_se_moments(self, nmom_max, integrals, **kwargs): """Build the moments of the self-energy. From 1e6a5bdde5aa9cbc5134ab90c88ebe3961db02ac Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 28 May 2024 17:15:14 +0100 Subject: [PATCH 28/33] Fix madelung attribute --- momentGW/pbc/ints.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index eb77ad1c..3c056f9d 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -2,8 +2,6 @@ Integral helpers with periodic boundary conditions. """ -from collections import defaultdict - import h5py import numpy as np from pyscf import lib @@ -663,15 +661,12 @@ def get_ewald(self, dm, basis="mo"): # Get the overlap matrix if basis == "mo": - ovlp = defaultdict(lambda: np.eye(self.nmo)) + ovlp = np.concatenate([np.eye(self.nmo) for _ in self.kpts]) else: ovlp = self.with_df.cell.pbc_intor("int1e_ovlp", hermi=1, kpts=self.kpts._kpts) - madelung = tools.pbc.madelung(self.with_df.cell, self.kpts._kpts) - # Initialise the Ewald matrix - ovlp = np.asarray(ovlp) - ew = madelung * util.einsum("kpq,kpi,kqj->kij", dm, ovlp.conj(), ovlp) + ew = self.madelung * util.einsum("kpq,kpi,kqj->kij", dm, np.conj(ovlp), ovlp) return ew @@ -747,7 +742,7 @@ def madelung(self): Return the Madelung constant for the lattice. """ if self._madelung is None: - self._madeling = tools.pbc.madelung(self.with_df.cell, self.kpts._kpts) + self._madelung = tools.pbc.madelung(self.with_df.cell, self.kpts._kpts) return self._madelung @property From 7534d99b368c983346d01c1ad5e3cb128aea0ebe Mon Sep 17 00:00:00 2001 From: mkakcl Date: Fri, 27 Sep 2024 15:57:34 +0100 Subject: [PATCH 29/33] HWB implementation --- momentGW/gw.py | 4 +- momentGW/pbc/base.py | 6 +- momentGW/pbc/gw.py | 50 ++++------ momentGW/pbc/ints.py | 68 +++++++++++++ momentGW/pbc/rpa.py | 222 ++++++++++++++++++++++++++++++++++++++----- momentGW/pbc/tda.py | 221 +++++++++++++++--------------------------- 6 files changed, 368 insertions(+), 203 deletions(-) diff --git a/momentGW/gw.py b/momentGW/gw.py index e1d715fc..8893481e 100644 --- a/momentGW/gw.py +++ b/momentGW/gw.py @@ -146,7 +146,7 @@ def name(self): @logging.with_timer("Static self-energy") @logging.with_status("Building static self-energy") - def build_se_static(self, integrals): + def build_se_static(self, integrals, **kwargs): """ Build the static part of the self-energy, including the Fock matrix. @@ -176,7 +176,7 @@ def build_se_static(self, integrals): veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask] vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] - vhf = integrals.get_veff(dm, j=vj, basis="ao") + vhf = integrals.get_veff(dm, j=vj, basis="ao", **kwargs) se_static = vhf - veff se_static = util.einsum( "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff diff --git a/momentGW/pbc/base.py b/momentGW/pbc/base.py index 016b8669..7f28e165 100644 --- a/momentGW/pbc/base.py +++ b/momentGW/pbc/base.py @@ -52,14 +52,14 @@ class BaseKGW(BaseGW): thc_opts : dict, optional Dictionary of options to be used for THC calculations. Current implementation requires a filepath to import the THC integrals. - fc : bool, optional + fsc : bool, optional If `True`, apply finite size corrections. Default value is `False`. """ _defaults = OrderedDict( **BaseGW._defaults, - fc=False, + fsc=None, ) _defaults["compression"] = None @@ -70,7 +70,7 @@ def __init__(self, mf, **kwargs): super().__init__(mf, **kwargs) # Options - self.fc = False + self.fsc = None # Attributes self._kpts = KPoints(self.cell, getattr(mf, "kpts", np.zeros((1, 3)))) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index ad531463..987f79d9 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -56,7 +56,7 @@ class KGW(BaseKGW, GW): thc_opts : dict, optional Dictionary of options to be used for THC calculations. Current implementation requires a filepath to import the THC integrals. - fc : bool, optional + fsc : bool, optional If `True`, apply finite size corrections. Default value is `False`. """ @@ -87,34 +87,24 @@ def build_se_static(self, integrals): Static part of the self-energy at each k-point. If `self.diagonal_se`, non-diagonal elements are set to zero. """ - if self.fc: - # Get intermediates - mask = self.active - dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff) - - # Get the contribution from the exchange-correlation potential - with util.SilentSCF(self._scf): - veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask] - vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask] - - vhf = integrals.get_veff(dm, j=vj, basis="ao", ewald=self.fc) - se_static = vhf - veff - se_static = util.einsum( - "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff - ) - # If diagonal approximation, set non-diagonal elements to zero - if self.diagonal_se: - se_static = util.einsum("...pq,pq->...pq", se_static, np.eye(se_static.shape[-1])) - - # Add the Fock matrix contribution - se_static += util.einsum( - "...p,...pq->...pq", self.mo_energy, np.eye(se_static.shape[-1]) + if self.fsc is not None: + if len(list(self.fsc))>3: + raise ValueError( + "Finite size corrections require as an input a combination of H, W and B " + "for the different finite size corrections (H - Head, W - Wing, B - Body)") + for i,letter in enumerate(list(self.fsc)): + if letter not in ["H","W","B"]: + raise ValueError( + "Finite size corrections require as an input a combination of H, W and B " + "for the different finite size corrections (H - Head, W - Wing, B - Body)") + kwargs = dict( + ewald=True, ) - - return se_static - else: - return super().build_se_static(integrals) + kwargs = dict( + ewald=False, + ) + return super().build_se_static(integrals, **kwargs) def build_se_moments(self, nmom_max, integrals, **kwargs): """Build the moments of the self-energy. @@ -139,10 +129,10 @@ def build_se_moments(self, nmom_max, integrals, **kwargs): """ if self.polarizability.lower() == "dtda": - tda = dTDA(self, nmom_max, integrals, fc=self.fc, **kwargs) + tda = dTDA(self, nmom_max, integrals, fsc=self.fsc, **kwargs) return tda.kernel() if self.polarizability.lower() == "drpa": - rpa = dRPA(self, nmom_max, integrals, **kwargs) + rpa = dRPA(self, nmom_max, integrals, fsc=self.fsc, **kwargs) return rpa.kernel() elif self.polarizability.lower() == "thc-dtda": tda = thc.dTDA(self, nmom_max, integrals, **kwargs) @@ -181,6 +171,8 @@ def ao2mo(self, transform=True): compression=self.compression, compression_tol=self.compression_tol, store_full=self.fock_loop, + mo_energy=self.mo_energy, + fc=self.fc, input_path=self.thc_opts["file_path"], ) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index eb77ad1c..e68ee911 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -10,6 +10,7 @@ from pyscf.ao2mo import _ao2mo from pyscf.pbc import tools from scipy.linalg import cholesky +from pyscf.pbc import dft from momentGW import logging, mpi_helper, util from momentGW.ints import Integrals, require_compression_metric @@ -47,6 +48,8 @@ def __init__( compression="ia", compression_tol=1e-10, store_full=False, + mo_energy=None, + fsc=None, input_path=None, ): Integrals.__init__( @@ -61,12 +64,15 @@ def __init__( # Options self.input_path = input_path + self.fsc = fsc # Attributes self.kpts = kpts self._madelung = None self._naux_full = None self._naux = None + if mo_energy is not None: + self.mo_energy_w = mo_energy["w"] @logging.with_status("Computing compression metric") def get_compression_metric(self): @@ -202,6 +208,12 @@ def transform(self, do_Lpq=None, do_Lpx=True, do_Lia=True): # ao2mo function for both real and complex integrals tao = np.empty([], dtype=np.int32) + # Building plane waves for finite size corrections + if self.fsc is not None: + q_abs = self.kpts.cell.get_abs_kpts(np.array([1e-3, 0, 0]).reshape(1, 3)) + HW_const = np.sqrt(4.0 * np.pi) / np.linalg.norm(q_abs[0]) + pw = HW_const * self.build_pert_term(q_abs[0]) + def _ao2mo_e2(Lpq, mo_coeff, orb_slice, out=None): mo_coeff = np.asarray(mo_coeff, order="F") if np.iscomplexobj(Lpq): @@ -215,6 +227,7 @@ def _ao2mo_e2(Lpq, mo_coeff, orb_slice, out=None): Lpx = {} Lia = {} Lai = {} + Mia = {} for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): @@ -338,6 +351,9 @@ def _ao2mo_e2(Lpq, mo_coeff, orb_slice, out=None): Lai[ki, kj] = Lai_k + if q == 0 and self.fsc is not None: + Mia[ki] = np.vstack([pw[ki], Lia[ki, ki]]) + # Store the arrays if do_Lpq: self._blocks["Lpq"] = Lpq @@ -346,6 +362,8 @@ def _ao2mo_e2(Lpq, mo_coeff, orb_slice, out=None): if do_Lia: self._blocks["Lia"] = Lia self._blocks["Lai"] = Lai + if self.fsc is not None: + self._blocks["Mia"] = Mia def get_cderi_from_thc(self): """ @@ -552,6 +570,9 @@ def get_k(self, dm, basis="mo", ewald=False): basis : str, optional Basis in which to build the K matrix. One of `("ao", "mo")`. Default value is `"mo"`. + ewald : bool, optional + Whether to include the Ewald exchange divergence. Default + value is `False`. Returns ------- @@ -675,6 +696,48 @@ def get_ewald(self, dm, basis="mo"): return ew + def build_pert_term(self, qpt): + """ + Build the charge-density density matrix at q-point index qpt + using perturbation theory. + + Parameters + ---------- + qpt : numpy.ndarray + q-point index representing the limit of our plane waves. + + Returns + ------- + pw_hw : numpy.ndarray + Charge-density density matrix in the long wavelength limit. + """ + + coords, weights = dft.gen_grid.get_becke_grids(self.kpts.cell, level=5) + + pw_hw = np.zeros((len(self.kpts),), dtype=object) + for k in self.kpts.loop(1): + ao_p = dft.numint.eval_ao(self.kpts.cell, coords, kpt=self.kpts[k], deriv=1) + ao, ao_grad = ao_p[0], ao_p[1:4] + + ao_ao_grad = util.einsum("g,gm,xgn->xmn", weights, ao.conj(), ao_grad) + q_ao_ao_grad = util.einsum("x,xmn->mn", qpt, ao_ao_grad) * -1.0j + q_mo_mo_grad = util.einsum( + "mn,mi,na->ia", + q_ao_ao_grad, + self.mo_coeff_w[k][:, self.mo_occ_w[k] > 0].conj(), + self.mo_coeff_w[k][:, self.mo_occ_w[k] == 0], + ) + + d = util.build_1h1p_energies( + (self.mo_energy_w[k], self.mo_energy_w[k]), + (self.mo_occ_w[k], self.mo_occ_w[k]), + ) + + dens = q_mo_mo_grad / d + pw_hw[k] = dens.flatten() / np.sqrt(self.kpts.cell.vol) + + return pw_hw + def get_jk(self, dm, **kwargs): """Build the J and K matrices. @@ -755,6 +818,11 @@ def Lai(self): """Get the full uncompressed ``(aux, MO, MO)`` integrals.""" return self._blocks["Lai"] + @property + def Mia(self): + """Get the finite size corrected uncompressed ``(1+ aux, MO, MO)`` integrals.""" + return self._blocks["Mia"] + @property def nmo(self): """Get the number of MOs.""" diff --git a/momentGW/pbc/rpa.py b/momentGW/pbc/rpa.py index 332a7944..3c460db7 100644 --- a/momentGW/pbc/rpa.py +++ b/momentGW/pbc/rpa.py @@ -71,14 +71,19 @@ def _build_diag_eri(self): for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): - kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) - diag_eri[q, kb] = ( - np.sum(np.abs(self.integrals.Lia[ki, kb]) ** 2, axis=0) / self.nkpts - ) + if q == 0 and self.fsc is not None: + diag_eri[q, ki] = ( + np.sum(np.abs(self.integrals.Mia[ki, ki]) ** 2, axis=0) / self.nkpts + ) + else: + kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) + diag_eri[q, kb] = ( + np.sum(np.abs(self.integrals.Lia[ki, kb]) ** 2, axis=0) / self.nkpts + ) return diag_eri - def _build_Liad(self, Lia, d): + def _build_Liad(self, Lia, d, corrected=False): """Construct the ``Liad`` array. Returns @@ -92,12 +97,15 @@ def _build_Liad(self, Lia, d): for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): - kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) - Liad[q, kb] = Lia[ki, kb] * d[q, kb] + if q == 0 and corrected: + Liad[q, ki] = self.integrals.Mia[ki] * d[q, ki] + else: + kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) + Liad[q, kb] = Lia[ki, kb] * d[q, kb] return Liad - def _build_Liadinv(self, Lia, d): + def _build_Liadinv(self, Lia, d, corrected=False): """Construct the ``Liadinv`` array. Returns @@ -111,8 +119,11 @@ def _build_Liadinv(self, Lia, d): for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): - kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) - Liadinv[q, kb] = Lia[ki, kb] / d[q, kb] + if q == 0 and corrected: + Liadinv[q, ki] = self.integrals.Mia[ki] / d[q, ki] + else: + kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) + Liadinv[q, kb] = Lia[ki, kb] / d[q, kb] return Liadinv @@ -136,16 +147,16 @@ def integrate(self): diag_eri = self._build_diag_eri() # Get the offset integral quadrature - quad = self.optimise_offset_quad(d, diag_eri) + quad_off = self.optimise_offset_quad(d, diag_eri) # Perform the offset integral - offset = self.eval_offset_integral(quad, d) + offset = self.eval_offset_integral(quad_off, d) # Get the main integral quadrature - quad = self.optimise_main_quad(d, diag_eri) + quad_main = self.optimise_main_quad(d, diag_eri) # Perform the main integral - integral = self.eval_main_integral(quad, d) + integral = self.eval_main_integral(quad_main, d) # Report quadrature error if self.report_quadrature_error: @@ -166,7 +177,12 @@ def integrate(self): f"(half = [{style_half}]{a:.3e}[/], quarter = [{style_quar}]{b:.3e}[/])", ) - return integral[0] + offset + if self.fsc is not None: + integral_fsc = self.eval_main_integral_fsc(quad_main, d) + offset_fsc = self.eval_offset_integral_fsc(quad_off, d) + return integral[0] + offset, integral_fsc[0] + offset_fsc + else: + return integral[0] + offset @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -184,7 +200,22 @@ def build_dd_moments(self, integral=None): moments : numpy.ndarray Moments of the density-density response. """ + integral, integral_fsc = self.integrate() + moments = self.build_uncorrected_dd_moments(integral) + if self.fsc is not None: + corrected_moments = self.build_corrected_dd_moments(integral_fsc) + for i in range(self.nmom_max + 1): + for kj in self.kpts.loop(1, mpi=True): + if "B" in self.fsc: + moments[0, kj, i] = corrected_moments[kj, i] + else: + moments[0, kj, i] = np.concatenate([corrected_moments[kj, i]], moments[0,kj, i]) + + return moments + else: + return moments + def build_uncorrected_dd_moments(self, integral=None): if integral is None: integral = self.integrate() @@ -205,8 +236,9 @@ def build_dd_moments(self, integral=None): inter = 0.0 for kj in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - tmp += np.dot(Liadinv[q, kb], self.integrals.Lia[kj, kb].conj().T) + tmp += np.dot(Liadinv[q, kb], Lia[kj, kb].conj().T) inter += np.dot(integral[q, kb], Liadinv[q, kb].T.conj()) + tmp = mpi_helper.allreduce(tmp) inter = mpi_helper.allreduce(inter) tmp *= 2.0 @@ -225,10 +257,10 @@ def build_dd_moments(self, integral=None): for i in range(2, self.nmom_max + 1): for q in kpts.loop(1): tmp = 0.0 - for ka in kpts.loop(1, mpi=True): - kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[ka])) + for kj in kpts.loop(1, mpi=True): + kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) moments[q, kb, i] = moments[q, kb, i - 2] * d[q, kb] ** 2 - tmp += np.dot(moments[q, kb, i - 2], self.integrals.Lia[ka, kb].conj().T) + tmp += np.dot(moments[q, kb, i - 2], Lia[kj, kb].conj().T) tmp = mpi_helper.allreduce(tmp) tmp /= self.nkpts tmp *= 2 @@ -238,6 +270,55 @@ def build_dd_moments(self, integral=None): return moments + def build_corrected_dd_moments(self, integral_fsc=None): + if integral_fsc is None: + integral, integral_fsc = self.integrate() + + kpts = self.kpts + Mia = self.integrals.Mia + moments = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + + # Construct the energy differences + d = self._build_d() + + # Calculate (L|ia) D_{ia} and (L|ia) D_{ia}^{-1} intermediates + Liad = self._build_Liad(self.integrals.Lia, d, corrected=True)[0] + Liadinv = self._build_Liadinv(self.integrals.Lia, d, corrected=True)[0] + + # Get the zeroth order moment + tmp = np.zeros((self.naux[0]+1, self.naux[0]+1), dtype=complex) + inter = 0.0 + for ka in kpts.loop(1, mpi=True): + tmp += np.dot(Liadinv[ka], Mia.T.conj()) + inter += np.dot(integral_fsc[ka], Liadinv[ka].T.conj()) + + tmp = mpi_helper.allreduce(tmp) + inter = mpi_helper.allreduce(inter) + tmp *= 2.0 + u = np.linalg.inv(np.eye(tmp.shape[0]) * self.nkpts / 2 + tmp) + + rest = np.dot(inter, u) * self.nkpts / 2 + for ki in kpts.loop(1, mpi=True): + moments[ki, 0] = integral_fsc[ki] / d[0, ki] * self.nkpts / 2 + moments[ki, 0] -= np.dot(rest, Liadinv[ki]) * 2 + + # Get the first order moment + moments[:, 1] = Liad / self.nkpts + + # Get the higher order moments + for i in range(2, self.nmom_max + 1): + tmp = 0.0 + for ka in kpts.loop(1, mpi=True): + moments[ka, i] = moments[ka, i - 2] * d[0, ka] ** 2 + tmp += np.dot(moments[ka, i - 2], Mia.conj().T) + tmp = mpi_helper.allreduce(tmp) + tmp /= self.nkpts + tmp *= 2 + for ki in kpts.loop(1, mpi=True): + moments[ki, i] += np.dot(tmp, Liad[ki]) * 2 + + return moments + def optimise_offset_quad(self, d, diag_eri, name="main"): """ Optimise the grid spacing of Gauss-Laguerre quadrature for the @@ -340,7 +421,7 @@ def eval_offset_integral(self, quad, d, Lia=None): if Lia is None: Lia = self.integrals.Lia Liad = self._build_Liad(Lia, d) - integrals = 2 * Liad / (self.nkpts**2) + integrals = 2 * Liad / (self.nkpts ** 2) kpts = self.kpts @@ -355,16 +436,55 @@ def eval_offset_integral(self, quad, d, Lia=None): lhs = mpi_helper.allreduce(lhs) lhs /= self.nkpts lhs *= 2 - for ka in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[ka])) rhs = self.integrals.Lia[ka, kb] * np.exp(-point * d[q, kb]) - rhs /= self.nkpts**2 + rhs /= self.nkpts ** 2 res = np.dot(lhs, rhs) integrals[q, kb] += res * weight * 4 return integrals + def eval_offset_integral_fsc(self, quad, d): + """Evaluate the offset integral. + + Parameters + ---------- + quad : tuple + The quadrature points and weights. + d : numpy.ndarray + Orbital energy differences at each k-point. + + Returns + ------- + integral : numpy.ndarray + Offset integral. + """ + + # Get the integral intermediates + Mia = self.integrals.Mia + Liad = self._build_Liad(self.integrals.Lia, d, corrected=True)[0] + integrals = 2 * Liad/ (self.nkpts**2) + + kpts = self.kpts + + # Calculate the integral for each point + for point, weight in zip(*quad): + lhs = 0.0 + for ka in kpts.loop(1, mpi=True): + expval = np.exp(-point * d[0, ka]) + lhs += np.dot(Liad[ka] * expval[None], Mia[ka].T.conj()) + lhs = mpi_helper.allreduce(lhs) + lhs /= self.nkpts + lhs *= 2 + for ka in kpts.loop(1, mpi=True): + rhs = Mia[ka] * np.exp(-point * d[0, ka]) + rhs /= self.nkpts ** 2 + res = np.dot(lhs,rhs) + integrals[ka] += res * weight * 4 + + return integrals + def optimise_main_quad(self, d, diag_eri, name="main"): """ Optimise the grid spacing of Clenshaw-Curtis quadrature for the @@ -511,7 +631,7 @@ def eval_main_integral(self, quad, d, Lia=None): qz = 0.0 for ki in kpts.loop(1, mpi=True): kj = kpts.member(kpts.wrap_around(kpts[q] + kpts[ki])) - f[kj] = 1.0 / (d[q, kj] ** 2 + point**2) + f[kj] = 1.0 / (d[q, kj] ** 2 + point ** 2) pre = (Lia[ki, kj] * f[kj]) * (4 / self.nkpts) qz += np.dot(pre, Liad[q, kj].T.conj()) qz = mpi_helper.allreduce(qz) @@ -521,8 +641,8 @@ def eval_main_integral(self, quad, d, Lia=None): for ka in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[ka])) - contrib[q, kb] = 2 * np.dot(inner, Lia[ka, kb]) / (self.nkpts**2) - value = weight * (contrib[q, kb] * f[kb] * (point**2 / np.pi)) + contrib[q, kb] = 2 * np.dot(inner, Lia[ka, kb]) / (self.nkpts ** 2) + value = weight * (contrib[q, kb] * f[kb] * (point ** 2 / np.pi)) integral[0, q, kb] += value if i % 2 == 0 and self.report_quadrature_error: @@ -531,3 +651,55 @@ def eval_main_integral(self, quad, d, Lia=None): integral[2, q, kb] += 4 * value return integral + + def eval_main_integral_fsc(self, quad, d): + """Evaluate the main integral. + + Parameters + ---------- + quad : tuple + The quadrature points and weights. + + Variables + ---------- + d : numpy.ndarray + Orbital energy differences at each k-point. + + Returns + ------- + integral : numpy.ndarray + Offset integral. + """ + + # Get the integral intermediates + Mia = self.integrals.Mia + Liad = self._build_Liad(self.integrals.Lia, d, corrected=True)[0] + + # Initialise the integral + dim = 3 if self.report_quadrature_error else 1 + integral = np.zeros((dim, self.nkpts), dtype=object) + + # Calculate the integral for each point + kpts = self.kpts + for i, (point, weight) in enumerate(zip(*quad)): + contrib = np.zeros_like(Liad, dtype=object) + f = np.zeros((self.nkpts), dtype=object) + qz = 0.0 + for ki in kpts.loop(1, mpi=True): + f[ki] = 1.0 / (d[0, ki] ** 2 + point ** 2) + pre = (Mia[ki] * f[ki]) * (4 / self.nkpts) + qz += np.dot(pre, Liad[ki].T.conj()) + qz = mpi_helper.allreduce(qz) + tmp = np.linalg.inv(np.eye(self.naux[0] + 1) + qz) - np.eye(self.naux[0] + 1) + inner = np.dot(qz, tmp) + for ki in kpts.loop(1, mpi=True): + contrib[ki] = 2 * np.dot(inner, Mia[ki]) / (self.nkpts ** 2) + value = weight * (contrib[ki] * f[ki] * (point ** 2 / np.pi)) + + integral[0, ki] += value + if i % 2 == 0 and self.report_quadrature_error: + integral[1, ki] += 2 * value + if i % 4 == 0 and self.report_quadrature_error: + integral[2, ki] += 4 * value + + return integral \ No newline at end of file diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index 7b03d2de..f1784253 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -4,7 +4,6 @@ import numpy as np import scipy.special -from pyscf.pbc import dft from momentGW import logging, mpi_helper, util from momentGW.tda import dTDA as MoldTDA @@ -42,12 +41,10 @@ def __init__( integrals, mo_energy=None, mo_occ=None, - fc=False, + fsc=False, ): super().__init__(gw, nmom_max, integrals, mo_energy=mo_energy, mo_occ=mo_occ) - self.fc = fc - if self.fc: - self.pw_hw = self.build_pert_term(self.q_abs[0]) + self.fsc = fsc @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -59,7 +56,21 @@ def build_dd_moments(self): moments : numpy.ndarray Moments of the density-density response at each k-point. """ + if self.fsc is not None: + corrected_moments = self.build_corrected_dd_moments() + moments = self.build_uncorrected_dd_moments() + for i in range(self.nmom_max + 1): + for kj in self.kpts.loop(1, mpi=True): + if "B" in self.fsc: + moments[0, kj, i] = corrected_moments[kj, i] + else: + moments[0, kj, i] = np.concatenate([corrected_moments[kj, i]], moments[0,kj, i]) + + return moments + else: + return self.build_uncorrected_dd_moments() + def build_uncorrected_dd_moments(self): # Initialise the moments kpts = self.kpts moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) @@ -96,7 +107,7 @@ def build_dd_moments(self): return moments - def build_head_dd_moments(self): + def build_corrected_dd_moments(self): """Build the head for the moments of the density-density response. Returns @@ -105,91 +116,89 @@ def build_head_dd_moments(self): Head moments of the density-density response at each k-point. """ kpts = self.kpts - head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - HW_const = np.sqrt(4.0 * np.pi) / np.linalg.norm(self.q_abs[0]) + moments = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) # Get the zeroth order moment for kj in kpts.loop(1, mpi=True): - head[kj, 0] += HW_const * self.pw_hw[kj].conj() + moments[kj, 0] += self.integrals.Mia/self.nkpts + # Get the higher order moments for i in range(1, self.nmom_max + 1): - tmp_head = np.zeros((self.naux[0]), dtype=complex) + tmp = np.zeros((self.naux[0]+1, self.naux[0]+1), dtype=complex) for kj in kpts.loop(1, mpi=True): d = util.build_1h1p_energies( (self.mo_energy_w[kj], self.mo_energy_w[kj]), (self.mo_occ_w[kj], self.mo_occ_w[kj]), ) - head[kj, i] += head[kj, i - 1] * d.ravel() - - tmp_head += util.einsum( - "a,aP->P", head[kj, i - 1], self.integrals.Lia[kj, kj].T.conj() - ) - - tmp_head = mpi_helper.allreduce(tmp_head) - tmp_head *= 2.0 - tmp_head /= self.nkpts + moments[kj, i] += moments[kj, i - 1] * d.ravel() + tmp += util.einsum("Pa,aQ->PQ", moments[kj, i - 1], self.integrals.Mia.T.conj()) + tmp = mpi_helper.allreduce(tmp) + tmp *= 2.0 + tmp /= self.nkpts for kj in kpts.loop(1, mpi=True): - head[kj, i] += util.einsum("P,Pa->a", tmp_head, self.integrals.Lia[kj, kj]) - - return head - - def _add_eta_fc_correction(self, moments_dd, eta): - """Correct the moments of the self-energy using the Head - and wings corrections. - - Parameters - ---------- - moments_dd : numpy.ndarray - Moments of the density-density response at each k-point. + moments[kj, i] += util.einsum("PQ,Qa->Pa", tmp, self.integrals.Mia) - eta : numpy.ndarray - Moments of the density-density response partly transformed - into moments of the screened Coulomb interaction at each - k-point. - - Returns - ------- - eta : numpy.ndarray - Corrected moments of the density-density response partly - transformed into moments of the screened Coulomb interaction - at each k-point. - - """ - head = self.build_head_dd_moments() + return moments + def build_eta(self, moments_dd): kpts = self.kpts - cell_vol = self.kpts.cell.vol - total_vol = self.kpts.cell.vol * self.nkpts - q0 = (6 * np.pi**2 / total_vol) ** (1 / 3) - norm_q_abs = np.linalg.norm(self.q_abs[0]) - HW_const = np.sqrt(4.0 * np.pi) / norm_q_abs + if self.fsc is not None: + cell_vol = kpts.cell.vol + q0 = (6 * np.pi ** 2 / (cell_vol * self.nkpts)) ** (1 / 3) - eta_head = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) - eta_wings = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) + # Setup dependent on diagonal SE + if self.gw.diagonal_se: + pqchar = pchar = qchar = "p" + eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo) + else: + pqchar, pchar, qchar = "pq", "p", "q" + eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) + eta = np.zeros((self.nkpts, self.nkpts), dtype=object) for n in range(self.nmom_max + 1): - for kj in kpts.loop(1, mpi=True): - eta_head[kj, n] += -HW_const * np.sum(head[kj, n] * self.pw_hw[kj]) - eta_wings[kj, n] += HW_const * util.einsum( - "Pa,a->P", moments_dd[0, kj, n], self.pw_hw[kj] - ) - for kx in kpts.loop(1, mpi=True): - for x in range(self.mo_energy_g[kx].size): - Lp = self.integrals.Lpx[kx, kx][:, :, x] + for q in kpts.loop(1): + eta_aux = 0 + for kj in kpts.loop(1, mpi=True): + kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) + eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) - original = eta[kx, 0][x, n] + eta_aux = mpi_helper.allreduce(eta_aux) + eta_aux *= 2.0 + #eta_aux /= self.nkpts - eta[kx, 0][x, n] += (2 / np.pi) * (q0) * eta_head[kx, n] * original + for kp in kpts.loop(1, mpi=True): + kx = kpts.member(kpts.wrap_around(kpts[kp] - kpts[q])) - wing_tmp = util.einsum("Pp,P->p", Lp, eta_wings[kx, n]) - wing_tmp = wing_tmp.real * 2 - wing_tmp *= -(np.sqrt(cell_vol / (4 * (np.pi**3))) * q0**2) + if not isinstance(eta[kp, q], np.ndarray): + eta[kp, q] = np.zeros(eta_shape(kx), dtype=eta_aux.dtype) - eta[kx, 0][x, n] += util.einsum("p,pq->pq", wing_tmp, original) - return eta + for x in range(self.mo_energy_g[kx].size): + Lp = self.integrals.Lpx[kp, kx][:, :, x] + subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" + if q==0 and self.fsc is not None: + eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux[1:,1:]/self.nkpts) + + if "H" in list(self.fsc): + if mpi_helper.rank == 0: + if self.gw.diagonal_se: + eta[kp, q][x, n][x] += (2/np.pi) * q0 * eta_aux[0,0] + else: + eta[kp, q][x, n][x, x] += (2/np.pi) * q0 * eta_aux[0,0] + if "W" in list(self.fsc): + wing_tmp = util.einsum(f"P,P{pchar}{qchar}->{pqchar}", eta_aux[0,1:], + self.integrals.Lpq[kp, kx]) + wing_tmp = ((q0 ** 2) * ( + (cell_vol / (4 * np.pi ** 3)) ** (1 / 2))) * wing_tmp.real + if self.gw.diagonal_se: # TODO: Check this + eta[kp, q][x, n][x] += 2*wing_tmp[x] + else: + eta[kp, q][x, n][x, :] += wing_tmp.T[x, :] + eta[kp, q][x, n][:, x] += wing_tmp[:, x] + else: + eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux/self.nkpts) def kernel(self, exact=False): """ @@ -315,90 +324,14 @@ def build_se_moments(self, moments_dd): moments_vir : numpy.ndarray Moments of the virtual self-energy at each k-point. """ - - kpts = self.kpts - - # Setup dependent on diagonal SE - if self.gw.diagonal_se: - pqchar = pchar = qchar = "p" - eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo) - else: - pqchar, pchar, qchar = "pq", "p", "q" - eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) - eta = np.zeros((self.nkpts, self.nkpts), dtype=object) - - # Get the moments in (aux|aux) and rotate to (mo|mo) - for n in range(self.nmom_max + 1): - for q in kpts.loop(1): - eta_aux = 0 - for kj in kpts.loop(1, mpi=True): - kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) - - eta_aux = mpi_helper.allreduce(eta_aux) - eta_aux *= 2.0 - eta_aux /= self.nkpts - - for kp in kpts.loop(1, mpi=True): - kx = kpts.member(kpts.wrap_around(kpts[kp] - kpts[q])) - - if not isinstance(eta[kp, q], np.ndarray): - eta[kp, q] = np.zeros(eta_shape(kx), dtype=eta_aux.dtype) - - for x in range(self.mo_energy_g[kx].size): - Lp = self.integrals.Lpx[kp, kx][:, :, x] - subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" - eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux) - - if self.fc: - eta = self._add_eta_fc_correction(moments_dd, eta) + eta = self.build_eta(moments_dd) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) return moments_occ, moments_vir - def build_pert_term(self, qpt): - """ - Build the charge-density density matrix at q-point index qpt - using perturbation theory. - - Parameters - ---------- - qpt : numpy.ndarray - q-point index representing the limit of our plane waves. - - Returns - ------- - pw_hw : numpy.ndarray - Charge-density density matrix in the long wavelength limit. - """ - - coords, weights = dft.gen_grid.get_becke_grids(self.gw.cell, level=5) - - pw_hw = np.zeros((self.nkpts,), dtype=object) - for k in self.kpts.loop(1): - ao_p = dft.numint.eval_ao(self.gw.cell, coords, kpt=self.kpts[k], deriv=1) - ao, ao_grad = ao_p[0], ao_p[1:4] - - ao_ao_grad = util.einsum("g,gm,xgn->xmn", weights, ao.conj(), ao_grad) - q_ao_ao_grad = util.einsum("x,xmn->mn", qpt, ao_ao_grad) * -1.0j - q_mo_mo_grad = util.einsum( - "mn,mi,na->ia", - q_ao_ao_grad, - self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] > 0].conj(), - self.integrals.mo_coeff_w[k][:, self.mo_occ_w[k] == 0], - ) - - d = util.build_1h1p_energies( - (self.mo_energy_w[k], self.mo_energy_w[k]), - (self.mo_occ_w[k], self.mo_occ_w[k]), - ) - - dens = q_mo_mo_grad / d - pw_hw[k] = dens.flatten() / np.sqrt(self.gw.cell.vol) - return pw_hw @property def naux(self): From e2941f734317e8de1ee0698362705bd5a5e70de9 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 30 Sep 2024 12:11:42 +0100 Subject: [PATCH 30/33] Corrections to HWB implementation --- momentGW/pbc/gw.py | 7 +++++-- momentGW/pbc/ints.py | 6 +++--- momentGW/pbc/rpa.py | 14 +++++++------- momentGW/pbc/tda.py | 28 +++++++++++++++++----------- 4 files changed, 32 insertions(+), 23 deletions(-) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index a58cc89c..650371d1 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -94,7 +94,10 @@ def get_veff(self, integrals, dm=None, **kwargs): # Set the default options if "ewald" not in kwargs: - kwargs = {**kwargs, "ewald": self.fc} + if self.fsc is not None: + kwargs = {**kwargs, "ewald": True} + else: + kwargs = {**kwargs, "ewald": False} # Get the effective potential veff = integrals.get_veff(dm, **kwargs) @@ -201,7 +204,7 @@ def ao2mo(self, transform=True): compression=self.compression, compression_tol=self.compression_tol, store_full=self.fock_loop, - mo_energy=self.mo_energy, + mo_energy_w=self.mo_energy, fsc=self.fsc, input_path=self.thc_opts["file_path"], ) diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 2452a3fd..3fc1cd08 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -46,7 +46,7 @@ def __init__( compression="ia", compression_tol=1e-10, store_full=False, - mo_energy=None, + mo_energy_w=None, fsc=None, input_path=None, ): @@ -69,8 +69,8 @@ def __init__( self._madelung = None self._naux_full = None self._naux = None - if mo_energy is not None: - self.mo_energy_w = mo_energy["w"] + if mo_energy_w is not None: + self.mo_energy_w = mo_energy_w @logging.with_status("Computing compression metric") def get_compression_metric(self): diff --git a/momentGW/pbc/rpa.py b/momentGW/pbc/rpa.py index 3c460db7..34ff44be 100644 --- a/momentGW/pbc/rpa.py +++ b/momentGW/pbc/rpa.py @@ -71,9 +71,9 @@ def _build_diag_eri(self): for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): - if q == 0 and self.fsc is not None: + if q == 0 and "B" in self.fsc: diag_eri[q, ki] = ( - np.sum(np.abs(self.integrals.Mia[ki, ki]) ** 2, axis=0) / self.nkpts + np.sum(np.abs(self.integrals.Mia[ki]) ** 2, axis=0) / self.nkpts ) else: kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) @@ -182,7 +182,7 @@ def integrate(self): offset_fsc = self.eval_offset_integral_fsc(quad_off, d) return integral[0] + offset, integral_fsc[0] + offset_fsc else: - return integral[0] + offset + return integral[0] + offset, None @logging.with_timer("Density-density moments") @logging.with_status("Constructing density-density moments") @@ -209,7 +209,7 @@ def build_dd_moments(self, integral=None): if "B" in self.fsc: moments[0, kj, i] = corrected_moments[kj, i] else: - moments[0, kj, i] = np.concatenate([corrected_moments[kj, i]], moments[0,kj, i]) + moments[0, kj, i] = np.concatenate((np.array([corrected_moments[kj, i][0,:]]), moments[0,kj, i])) return moments else: @@ -289,7 +289,7 @@ def build_corrected_dd_moments(self, integral_fsc=None): tmp = np.zeros((self.naux[0]+1, self.naux[0]+1), dtype=complex) inter = 0.0 for ka in kpts.loop(1, mpi=True): - tmp += np.dot(Liadinv[ka], Mia.T.conj()) + tmp += np.dot(Liadinv[ka], Mia[ka].T.conj()) inter += np.dot(integral_fsc[ka], Liadinv[ka].T.conj()) tmp = mpi_helper.allreduce(tmp) @@ -310,7 +310,7 @@ def build_corrected_dd_moments(self, integral_fsc=None): tmp = 0.0 for ka in kpts.loop(1, mpi=True): moments[ka, i] = moments[ka, i - 2] * d[0, ka] ** 2 - tmp += np.dot(moments[ka, i - 2], Mia.conj().T) + tmp += np.dot(moments[ka, i - 2], Mia[ka].conj().T) tmp = mpi_helper.allreduce(tmp) tmp /= self.nkpts tmp *= 2 @@ -702,4 +702,4 @@ def eval_main_integral_fsc(self, quad, d): if i % 4 == 0 and self.report_quadrature_error: integral[2, ki] += 4 * value - return integral \ No newline at end of file + return integral diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index f1784253..bc2dd810 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -64,7 +64,7 @@ def build_dd_moments(self): if "B" in self.fsc: moments[0, kj, i] = corrected_moments[kj, i] else: - moments[0, kj, i] = np.concatenate([corrected_moments[kj, i]], moments[0,kj, i]) + moments[0, kj, i] = np.concatenate((np.array([corrected_moments[kj, i][0,:]]), moments[0,kj, i])) return moments else: @@ -120,7 +120,7 @@ def build_corrected_dd_moments(self): # Get the zeroth order moment for kj in kpts.loop(1, mpi=True): - moments[kj, 0] += self.integrals.Mia/self.nkpts + moments[kj, 0] += self.integrals.Mia[kj]/self.nkpts # Get the higher order moments @@ -132,13 +132,13 @@ def build_corrected_dd_moments(self): (self.mo_occ_w[kj], self.mo_occ_w[kj]), ) moments[kj, i] += moments[kj, i - 1] * d.ravel() - tmp += util.einsum("Pa,aQ->PQ", moments[kj, i - 1], self.integrals.Mia.T.conj()) + tmp += util.einsum("Pa,aQ->PQ", moments[kj, i - 1], self.integrals.Mia[kj].T.conj()) tmp = mpi_helper.allreduce(tmp) tmp *= 2.0 tmp /= self.nkpts for kj in kpts.loop(1, mpi=True): - moments[kj, i] += util.einsum("PQ,Qa->Pa", tmp, self.integrals.Mia) + moments[kj, i] += util.einsum("PQ,Qa->Pa", tmp, self.integrals.Mia[kj]) return moments @@ -162,8 +162,11 @@ def build_eta(self, moments_dd): for q in kpts.loop(1): eta_aux = 0 for kj in kpts.loop(1, mpi=True): - kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) + if q==0 and self.fsc is not None: + eta_aux += np.dot(moments_dd[q, kj, n], self.integrals.Mia[kj].T.conj()) + else: + kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) + eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) eta_aux = mpi_helper.allreduce(eta_aux) eta_aux *= 2.0 @@ -181,25 +184,28 @@ def build_eta(self, moments_dd): if q==0 and self.fsc is not None: eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux[1:,1:]/self.nkpts) - if "H" in list(self.fsc): + if "H" in self.fsc: if mpi_helper.rank == 0: if self.gw.diagonal_se: eta[kp, q][x, n][x] += (2/np.pi) * q0 * eta_aux[0,0] else: eta[kp, q][x, n][x, x] += (2/np.pi) * q0 * eta_aux[0,0] - if "W" in list(self.fsc): + if "W" in self.fsc: wing_tmp = util.einsum(f"P,P{pchar}{qchar}->{pqchar}", eta_aux[0,1:], - self.integrals.Lpq[kp, kx]) + self.integrals.Lpx[kp, kx]) wing_tmp = ((q0 ** 2) * ( (cell_vol / (4 * np.pi ** 3)) ** (1 / 2))) * wing_tmp.real if self.gw.diagonal_se: # TODO: Check this eta[kp, q][x, n][x] += 2*wing_tmp[x] else: - eta[kp, q][x, n][x, :] += wing_tmp.T[x, :] - eta[kp, q][x, n][:, x] += wing_tmp[:, x] + eta[kp, q][x, n][x, :] -= wing_tmp.T[x, :] + eta[kp, q][x, n][:, x] -= wing_tmp[:, x] else: eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux/self.nkpts) + return eta + + def kernel(self, exact=False): """ Run the polarizability calculation to compute moments of the From d7ebf2d88c6bd9d986829456d74f80a7c7ae4030 Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 30 Sep 2024 16:16:26 +0100 Subject: [PATCH 31/33] Example and tests for HWB corrections --- examples/20-finite_size_corrections.py | 46 ++++++++++++++++++++++++++ tests/test_kgw.py | 23 ++++++++++--- 2 files changed, 64 insertions(+), 5 deletions(-) create mode 100644 examples/20-finite_size_corrections.py diff --git a/examples/20-finite_size_corrections.py b/examples/20-finite_size_corrections.py new file mode 100644 index 00000000..98a14c38 --- /dev/null +++ b/examples/20-finite_size_corrections.py @@ -0,0 +1,46 @@ +""" +Examples of finite size corrections for `momentGW` calculations for periodic solids. +""" + +import numpy as np +from pyscf.pbc import dft, gto + +from momentGW import KGW, KUGW + +# Define a unit cell +cell = gto.Cell() +cell.atom = "He 0 0 0; He 1 1 1" +cell.basis = "6-31g" +cell.a = np.eye(3) * 3 +cell.verbose = 5 +cell.build() +kpts = cell.make_kpts([3, 1, 1]) + +# Run a DFT calculation +mf = dft.KRKS(cell, kpts) +mf = mf.density_fit() +mf.xc = "b3lyp" +mf.kernel() + +# For any `pbc` calculation a finite size correction can be added for +# the divergence associated with G = 0, q=0 for GDF Coulomb integrals. +# This is done by setting the `fsc` attribute to any combination of +# "H", "W", and/or "B" for the Head, Wings and Body portion of the +# correction. The Ewald correction is added in all cases. The default +# is `None` which disables the correction. + +# RHF reference +gw = KGW(mf) +gw.polarizability = "dRPA" +gw.fsc = "HWB" +gw.kernel(nmom_max=3) + +gw = KGW(mf) +gw.polarizability = "dTDA" +gw.fsc = "HW" +gw.kernel(nmom_max=3) + +gw = KGW(mf) +gw.polarizability = "dRPA" +gw.fsc = "H" +gw.kernel(nmom_max=3) \ No newline at end of file diff --git a/tests/test_kgw.py b/tests/test_kgw.py index 72b70656..eb567ee0 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -127,14 +127,27 @@ def test_dtda_HW_regression(self): nmom_max = 5 kgw = KGW(self.mf) kgw.polarizability = "dtda" - kgw.fc = True + kgw.fsc = "HW" conv, gf, se, _ = kgw.kernel(nmom_max) gf_occ = gf[0].occupied().physical(weight=1e-1) gf_vir = gf[0].virtual().physical(weight=1e-1) - self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7516810595, 6) - self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8925838943, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.0900606812, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8839090190, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7513153970368053, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8922171344774273, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.08971582402363, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8836405468852633, 6) + + def test_drpa_HW_regression(self): + nmom_max = 5 + kgw = KGW(self.mf) + kgw.polarizability = "drpa" + kgw.fsc = "HWB" + conv, gf, se, _ = kgw.kernel(nmom_max) + gf_occ = gf[0].occupied().physical(weight=1e-1) + gf_vir = gf[0].virtual().physical(weight=1e-1) + self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7597568334841686, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.9002286289843537, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.09266415629468, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8899834941075897, 6) def test_dtda_vs_supercell_fock_loop(self): nmom_max = 5 From 7f361eee746c29d8730fc45d9e500c8a9593ec9c Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 30 Sep 2024 16:31:26 +0100 Subject: [PATCH 32/33] Small fix and updated test values --- momentGW/pbc/rpa.py | 2 +- tests/test_kgw.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/momentGW/pbc/rpa.py b/momentGW/pbc/rpa.py index 34ff44be..e8086a02 100644 --- a/momentGW/pbc/rpa.py +++ b/momentGW/pbc/rpa.py @@ -71,7 +71,7 @@ def _build_diag_eri(self): for q in self.kpts.loop(1): for ki in self.kpts.loop(1, mpi=True): - if q == 0 and "B" in self.fsc: + if q == 0 and self.fsc is not None and "B" in self.fsc: diag_eri[q, ki] = ( np.sum(np.abs(self.integrals.Mia[ki]) ** 2, axis=0) / self.nkpts ) diff --git a/tests/test_kgw.py b/tests/test_kgw.py index eb567ee0..0202f92a 100644 --- a/tests/test_kgw.py +++ b/tests/test_kgw.py @@ -131,10 +131,10 @@ def test_dtda_HW_regression(self): conv, gf, se, _ = kgw.kernel(nmom_max) gf_occ = gf[0].occupied().physical(weight=1e-1) gf_vir = gf[0].virtual().physical(weight=1e-1) - self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7513153970368053, 6) - self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8922171344774273, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.08971582402363, 6) - self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.8836405468852633, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-1]), -0.7515760991553542, 6) + self.assertAlmostEqual(np.max(gf_occ.energies[-2]), -0.8924789480797911, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[0]), 1.089967613154295, 6) + self.assertAlmostEqual(np.max(gf_vir.energies[1]), 1.883886215615475, 6) def test_drpa_HW_regression(self): nmom_max = 5 From b2c0fe7cefee9bc8af1b68ac8e9783d1fd1889ae Mon Sep 17 00:00:00 2001 From: mkakcl Date: Mon, 30 Sep 2024 17:38:49 +0100 Subject: [PATCH 33/33] Some changes related to the old PR and linting. --- examples/20-finite_size_corrections.py | 2 +- momentGW/pbc/gw.py | 6 +- momentGW/pbc/ints.py | 10 +- momentGW/pbc/rpa.py | 30 ++--- momentGW/pbc/tda.py | 162 +++++++++++++------------ 5 files changed, 111 insertions(+), 99 deletions(-) diff --git a/examples/20-finite_size_corrections.py b/examples/20-finite_size_corrections.py index 98a14c38..523faeea 100644 --- a/examples/20-finite_size_corrections.py +++ b/examples/20-finite_size_corrections.py @@ -43,4 +43,4 @@ gw = KGW(mf) gw.polarizability = "dRPA" gw.fsc = "H" -gw.kernel(nmom_max=3) \ No newline at end of file +gw.kernel(nmom_max=3) diff --git a/momentGW/pbc/gw.py b/momentGW/pbc/gw.py index 650371d1..c4809171 100644 --- a/momentGW/pbc/gw.py +++ b/momentGW/pbc/gw.py @@ -125,12 +125,12 @@ def build_se_static(self, integrals, **kwargs): `self.diagonal_se`, non-diagonal elements are set to zero. """ if self.fsc is not None: - if len(list(self.fsc))>3: + if len(list(self.fsc)) > 3: raise ValueError( "Finite size corrections require as an input a combination of H, W and B " "for the different finite size corrections (H - Head, W - Wing, B - Body)") - for i,letter in enumerate(list(self.fsc)): - if letter not in ["H","W","B"]: + for i, letter in enumerate(list(self.fsc)): + if letter not in ["H", "W", "B"]: raise ValueError( "Finite size corrections require as an input a combination of H, W and B " "for the different finite size corrections (H - Head, W - Wing, B - Body)") diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py index 3fc1cd08..34953184 100644 --- a/momentGW/pbc/ints.py +++ b/momentGW/pbc/ints.py @@ -6,9 +6,8 @@ import numpy as np from pyscf import lib from pyscf.ao2mo import _ao2mo -from pyscf.pbc import tools +from pyscf.pbc import dft, tools from scipy.linalg import cholesky -from pyscf.pbc import dft from momentGW import logging, mpi_helper, util from momentGW.ints import Integrals, require_compression_metric @@ -209,8 +208,8 @@ def transform(self, do_Lpq=None, do_Lpx=True, do_Lia=True): # Building plane waves for finite size corrections if self.fsc is not None: q_abs = self.kpts.cell.get_abs_kpts(np.array([1e-3, 0, 0]).reshape(1, 3)) - HW_const = np.sqrt(4.0 * np.pi) / np.linalg.norm(q_abs[0]) - pw = HW_const * self.build_pert_term(q_abs[0]) + hwb_const = np.sqrt(4.0 * np.pi) / np.linalg.norm(q_abs[0]) + pw = hwb_const * self.build_pert_term(q_abs[0]) def _ao2mo_e2(Lpq, mo_coeff, orb_slice, out=None): mo_coeff = np.asarray(mo_coeff, order="F") @@ -815,7 +814,8 @@ def Lai(self): @property def Mia(self): - """Get the finite size corrected uncompressed ``(1+ aux, MO, MO)`` integrals.""" + """Get the finite size corrected uncompressed ``(1+ aux, MO, MO)`` + integrals.""" return self._blocks["Mia"] @property diff --git a/momentGW/pbc/rpa.py b/momentGW/pbc/rpa.py index e8086a02..7e3deb68 100644 --- a/momentGW/pbc/rpa.py +++ b/momentGW/pbc/rpa.py @@ -73,7 +73,7 @@ def _build_diag_eri(self): for ki in self.kpts.loop(1, mpi=True): if q == 0 and self.fsc is not None and "B" in self.fsc: diag_eri[q, ki] = ( - np.sum(np.abs(self.integrals.Mia[ki]) ** 2, axis=0) / self.nkpts + np.sum(np.abs(self.integrals.Mia[ki]) ** 2, axis=0) / self.nkpts ) else: kb = self.kpts.member(self.kpts.wrap_around(self.kpts[q] + self.kpts[ki])) @@ -209,7 +209,9 @@ def build_dd_moments(self, integral=None): if "B" in self.fsc: moments[0, kj, i] = corrected_moments[kj, i] else: - moments[0, kj, i] = np.concatenate((np.array([corrected_moments[kj, i][0,:]]), moments[0,kj, i])) + moments[0, kj, i] = np.concatenate( + (np.array([corrected_moments[kj, i][0, :]]), moments[0, kj, i]) + ) return moments else: @@ -286,7 +288,7 @@ def build_corrected_dd_moments(self, integral_fsc=None): Liadinv = self._build_Liadinv(self.integrals.Lia, d, corrected=True)[0] # Get the zeroth order moment - tmp = np.zeros((self.naux[0]+1, self.naux[0]+1), dtype=complex) + tmp = np.zeros((self.naux[0] + 1, self.naux[0] + 1), dtype=complex) inter = 0.0 for ka in kpts.loop(1, mpi=True): tmp += np.dot(Liadinv[ka], Mia[ka].T.conj()) @@ -421,7 +423,7 @@ def eval_offset_integral(self, quad, d, Lia=None): if Lia is None: Lia = self.integrals.Lia Liad = self._build_Liad(Lia, d) - integrals = 2 * Liad / (self.nkpts ** 2) + integrals = 2 * Liad / (self.nkpts**2) kpts = self.kpts @@ -439,7 +441,7 @@ def eval_offset_integral(self, quad, d, Lia=None): for ka in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[ka])) rhs = self.integrals.Lia[ka, kb] * np.exp(-point * d[q, kb]) - rhs /= self.nkpts ** 2 + rhs /= self.nkpts**2 res = np.dot(lhs, rhs) integrals[q, kb] += res * weight * 4 @@ -464,7 +466,7 @@ def eval_offset_integral_fsc(self, quad, d): # Get the integral intermediates Mia = self.integrals.Mia Liad = self._build_Liad(self.integrals.Lia, d, corrected=True)[0] - integrals = 2 * Liad/ (self.nkpts**2) + integrals = 2 * Liad / (self.nkpts**2) kpts = self.kpts @@ -479,8 +481,8 @@ def eval_offset_integral_fsc(self, quad, d): lhs *= 2 for ka in kpts.loop(1, mpi=True): rhs = Mia[ka] * np.exp(-point * d[0, ka]) - rhs /= self.nkpts ** 2 - res = np.dot(lhs,rhs) + rhs /= self.nkpts**2 + res = np.dot(lhs, rhs) integrals[ka] += res * weight * 4 return integrals @@ -631,7 +633,7 @@ def eval_main_integral(self, quad, d, Lia=None): qz = 0.0 for ki in kpts.loop(1, mpi=True): kj = kpts.member(kpts.wrap_around(kpts[q] + kpts[ki])) - f[kj] = 1.0 / (d[q, kj] ** 2 + point ** 2) + f[kj] = 1.0 / (d[q, kj] ** 2 + point**2) pre = (Lia[ki, kj] * f[kj]) * (4 / self.nkpts) qz += np.dot(pre, Liad[q, kj].T.conj()) qz = mpi_helper.allreduce(qz) @@ -641,8 +643,8 @@ def eval_main_integral(self, quad, d, Lia=None): for ka in kpts.loop(1, mpi=True): kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[ka])) - contrib[q, kb] = 2 * np.dot(inner, Lia[ka, kb]) / (self.nkpts ** 2) - value = weight * (contrib[q, kb] * f[kb] * (point ** 2 / np.pi)) + contrib[q, kb] = 2 * np.dot(inner, Lia[ka, kb]) / (self.nkpts**2) + value = weight * (contrib[q, kb] * f[kb] * (point**2 / np.pi)) integral[0, q, kb] += value if i % 2 == 0 and self.report_quadrature_error: @@ -686,15 +688,15 @@ def eval_main_integral_fsc(self, quad, d): f = np.zeros((self.nkpts), dtype=object) qz = 0.0 for ki in kpts.loop(1, mpi=True): - f[ki] = 1.0 / (d[0, ki] ** 2 + point ** 2) + f[ki] = 1.0 / (d[0, ki] ** 2 + point**2) pre = (Mia[ki] * f[ki]) * (4 / self.nkpts) qz += np.dot(pre, Liad[ki].T.conj()) qz = mpi_helper.allreduce(qz) tmp = np.linalg.inv(np.eye(self.naux[0] + 1) + qz) - np.eye(self.naux[0] + 1) inner = np.dot(qz, tmp) for ki in kpts.loop(1, mpi=True): - contrib[ki] = 2 * np.dot(inner, Mia[ki]) / (self.nkpts ** 2) - value = weight * (contrib[ki] * f[ki] * (point ** 2 / np.pi)) + contrib[ki] = 2 * np.dot(inner, Mia[ki]) / (self.nkpts**2) + value = weight * (contrib[ki] * f[ki] * (point**2 / np.pi)) integral[0, ki] += value if i % 2 == 0 and self.report_quadrature_error: diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py index bc2dd810..e8eff868 100644 --- a/momentGW/pbc/tda.py +++ b/momentGW/pbc/tda.py @@ -56,21 +56,34 @@ def build_dd_moments(self): moments : numpy.ndarray Moments of the density-density response at each k-point. """ + + moments = self.build_uncorrected_dd_moments() + if self.fsc is not None: corrected_moments = self.build_corrected_dd_moments() - moments = self.build_uncorrected_dd_moments() for i in range(self.nmom_max + 1): for kj in self.kpts.loop(1, mpi=True): if "B" in self.fsc: moments[0, kj, i] = corrected_moments[kj, i] else: - moments[0, kj, i] = np.concatenate((np.array([corrected_moments[kj, i][0,:]]), moments[0,kj, i])) + moments[0, kj, i] = np.concatenate( + (np.array([corrected_moments[kj, i][0, :]]), moments[0, kj, i]) + ) return moments else: - return self.build_uncorrected_dd_moments() + return moments def build_uncorrected_dd_moments(self): + """Build the moments of the density-density response for + all k-point pairs. + + Returns + ------- + moments : numpy.ndarray + Moments of the density-density response at each k-point pair. + """ + # Initialise the moments kpts = self.kpts moments = np.zeros((self.nkpts, self.nkpts, self.nmom_max + 1), dtype=object) @@ -108,24 +121,25 @@ def build_uncorrected_dd_moments(self): return moments def build_corrected_dd_moments(self): - """Build the head for the moments of the density-density response. + """Build the finite size corrections for the moments of the + density-density response. Returns ------- - head_moments : numpy.ndarray - Head moments of the density-density response at each k-point. + moments : numpy.ndarray + Corrected moments of the density-density response at each + k-point for q=0. """ kpts = self.kpts moments = np.zeros((self.nkpts, self.nmom_max + 1), dtype=object) # Get the zeroth order moment for kj in kpts.loop(1, mpi=True): - moments[kj, 0] += self.integrals.Mia[kj]/self.nkpts - + moments[kj, 0] += self.integrals.Mia[kj] / self.nkpts # Get the higher order moments for i in range(1, self.nmom_max + 1): - tmp = np.zeros((self.naux[0]+1, self.naux[0]+1), dtype=complex) + tmp = np.zeros((self.naux[0] + 1, self.naux[0] + 1), dtype=complex) for kj in kpts.loop(1, mpi=True): d = util.build_1h1p_energies( (self.mo_energy_w[kj], self.mo_energy_w[kj]), @@ -142,70 +156,6 @@ def build_corrected_dd_moments(self): return moments - def build_eta(self, moments_dd): - kpts = self.kpts - - if self.fsc is not None: - cell_vol = kpts.cell.vol - q0 = (6 * np.pi ** 2 / (cell_vol * self.nkpts)) ** (1 / 3) - - # Setup dependent on diagonal SE - if self.gw.diagonal_se: - pqchar = pchar = qchar = "p" - eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo) - else: - pqchar, pchar, qchar = "pq", "p", "q" - eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) - eta = np.zeros((self.nkpts, self.nkpts), dtype=object) - - for n in range(self.nmom_max + 1): - for q in kpts.loop(1): - eta_aux = 0 - for kj in kpts.loop(1, mpi=True): - if q==0 and self.fsc is not None: - eta_aux += np.dot(moments_dd[q, kj, n], self.integrals.Mia[kj].T.conj()) - else: - kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) - eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) - - eta_aux = mpi_helper.allreduce(eta_aux) - eta_aux *= 2.0 - #eta_aux /= self.nkpts - - for kp in kpts.loop(1, mpi=True): - kx = kpts.member(kpts.wrap_around(kpts[kp] - kpts[q])) - - if not isinstance(eta[kp, q], np.ndarray): - eta[kp, q] = np.zeros(eta_shape(kx), dtype=eta_aux.dtype) - - for x in range(self.mo_energy_g[kx].size): - Lp = self.integrals.Lpx[kp, kx][:, :, x] - subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" - if q==0 and self.fsc is not None: - eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux[1:,1:]/self.nkpts) - - if "H" in self.fsc: - if mpi_helper.rank == 0: - if self.gw.diagonal_se: - eta[kp, q][x, n][x] += (2/np.pi) * q0 * eta_aux[0,0] - else: - eta[kp, q][x, n][x, x] += (2/np.pi) * q0 * eta_aux[0,0] - if "W" in self.fsc: - wing_tmp = util.einsum(f"P,P{pchar}{qchar}->{pqchar}", eta_aux[0,1:], - self.integrals.Lpx[kp, kx]) - wing_tmp = ((q0 ** 2) * ( - (cell_vol / (4 * np.pi ** 3)) ** (1 / 2))) * wing_tmp.real - if self.gw.diagonal_se: # TODO: Check this - eta[kp, q][x, n][x] += 2*wing_tmp[x] - else: - eta[kp, q][x, n][x, :] -= wing_tmp.T[x, :] - eta[kp, q][x, n][:, x] -= wing_tmp[:, x] - else: - eta[kp, q][x, n] += util.einsum(subscript, Lp, Lp.conj(), eta_aux/self.nkpts) - - return eta - - def kernel(self, exact=False): """ Run the polarizability calculation to compute moments of the @@ -330,15 +280,75 @@ def build_se_moments(self, moments_dd): moments_vir : numpy.ndarray Moments of the virtual self-energy at each k-point. """ - eta = self.build_eta(moments_dd) + kpts = self.kpts + + if self.fsc is not None: + q0 = (6 * np.pi**2 / (kpts.cell.vol * self.nkpts)) ** (1 / 3) + + # Setup dependent on diagonal SE + if self.gw.diagonal_se: + pqchar = pchar = qchar = "p" + eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo) + else: + pqchar, pchar, qchar = "pq", "p", "q" + eta_shape = lambda k: (self.mo_energy_g[k].size, self.nmom_max + 1, self.nmo, self.nmo) + eta = np.zeros((self.nkpts, self.nkpts), dtype=object) + + for n in range(self.nmom_max + 1): + for q in kpts.loop(1): + eta_aux = 0 + for kj in kpts.loop(1, mpi=True): + if q == 0 and self.fsc is not None: + eta_aux += np.dot(moments_dd[q, kj, n], self.integrals.Mia[kj].T.conj()) + else: + kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj])) + eta_aux += np.dot(moments_dd[q, kb, n], self.integrals.Lia[kj, kb].T.conj()) + + eta_aux = mpi_helper.allreduce(eta_aux) + eta_aux *= 2.0 + + for kp in kpts.loop(1, mpi=True): + kx = kpts.member(kpts.wrap_around(kpts[kp] - kpts[q])) + + if not isinstance(eta[kp, q], np.ndarray): + eta[kp, q] = np.zeros(eta_shape(kx), dtype=eta_aux.dtype) + + for x in range(self.mo_energy_g[kx].size): + Lp = self.integrals.Lpx[kp, kx][:, :, x] + subscript = f"P{pchar},Q{qchar},PQ->{pqchar}" + if q == 0 and self.fsc is not None: + eta[kp, q][x, n] += util.einsum( + subscript, Lp, Lp.conj(), eta_aux[1:, 1:] / self.nkpts + ) + if "H" in self.fsc and mpi_helper.rank == 0: + if self.gw.diagonal_se: + eta[kp, q][x, n][x] += (2 / np.pi) * q0 * eta_aux[0, 0] + else: + eta[kp, q][x, n][x, x] += (2 / np.pi) * q0 * eta_aux[0, 0] + if "W" in self.fsc: + wing_tmp = util.einsum( + f"P,P{pchar}{qchar}->{pqchar}", + eta_aux[0, 1:], + self.integrals.Lpx[kp, kx], + ) + wing_tmp = ( + (q0**2) * ((kpts.cell.vol / (4 * np.pi**3)) ** (1 / 2)) + ) * wing_tmp.real + if self.gw.diagonal_se: + eta[kp, q][x, n][x] += 2 * wing_tmp[x] + else: + eta[kp, q][x, n][x, :] -= wing_tmp.T[x, :] + eta[kp, q][x, n][:, x] -= wing_tmp[:, x] + else: + eta[kp, q][x, n] += util.einsum( + subscript, Lp, Lp.conj(), eta_aux / self.nkpts + ) # Construct the self-energy moments moments_occ, moments_vir = self.convolve(eta) return moments_occ, moments_vir - - @property def naux(self): """Number of auxiliaries."""