diff --git a/examples/20-finite_size_corrections.py b/examples/20-finite_size_corrections.py
new file mode 100644
index 00000000..523faeea
--- /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)
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/base.py b/momentGW/pbc/base.py
index cf2311c0..146631c4 100644
--- a/momentGW/pbc/base.py
+++ b/momentGW/pbc/base.py
@@ -53,14 +53,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
 
@@ -71,7 +71,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 f0991ea6..c4809171 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`.
     """
@@ -69,9 +69,44 @@ 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:
+            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)
+
+        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 +115,8 @@ def build_se_static(self, integrals):
         ----------
         integrals : KIntegrals
             Integrals object.
+        **kwargs : dict, optional
+            Additional keyword arguments.
 
         Returns
         -------
@@ -87,7 +124,20 @@ 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.
         """
-        return super().build_se_static(integrals)
+        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 = {**kwargs, "force_build": True}
+        else:
+            kwargs = {**kwargs, "force_build": False}
+        return super().build_se_static(integrals, **kwargs)
 
     def build_se_moments(self, nmom_max, integrals, **kwargs):
         """Build the moments of the self-energy.
@@ -112,10 +162,10 @@ def build_se_moments(self, nmom_max, integrals, **kwargs):
         """
 
         if self.polarizability.lower() == "dtda":
-            tda = dTDA(self, nmom_max, integrals, **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)
@@ -154,6 +204,8 @@ def ao2mo(self, transform=True):
                 compression=self.compression,
                 compression_tol=self.compression_tol,
                 store_full=self.fock_loop,
+                mo_energy_w=self.mo_energy,
+                fsc=self.fsc,
                 input_path=self.thc_opts["file_path"],
             )
 
@@ -334,7 +386,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 +419,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 bab2c5e4..86b6d55d 100644
--- a/momentGW/pbc/ints.py
+++ b/momentGW/pbc/ints.py
@@ -4,12 +4,11 @@
 
 import functools
 from collections import defaultdict
-
 import h5py
 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 momentGW import logging, mpi_helper, util
@@ -48,6 +47,8 @@ def __init__(
         compression="ia",
         compression_tol=1e-10,
         store_full=False,
+        mo_energy_w=None,
+        fsc=None,
         input_path=None,
     ):
         Integrals.__init__(
@@ -62,10 +63,13 @@ def __init__(
 
         # Options
         self.input_path = input_path
+        self.fsc = fsc
 
         # Attributes
         self.kpts = kpts
         self._naux = None
+        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):
@@ -209,6 +213,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))
+            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")
             if np.iscomplexobj(Lpq):
@@ -222,6 +232,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):
@@ -323,6 +334,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
@@ -331,6 +345,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):
         """
@@ -390,7 +406,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)
 
@@ -400,7 +416,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(naux[q], -1)
+                tmp = tmp.reshape(self.naux[invq], -1)
                 Lai_k += tmp
 
                 Lai[ki, kj] = Lai_k
@@ -544,6 +560,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
         -------
@@ -655,15 +674,57 @@ 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)
 
         # Initialise the Ewald matrix
-        ew = 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
 
+    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.
 
@@ -724,13 +785,19 @@ def get_fock(self, dm, h1e, **kwargs):
         """
         return super().get_fock(dm, h1e, **kwargs)
 
-    @functools.cached_property
+    def reciprocal_lattice(self):
+        """
+        Return the reciprocal lattice vectors.
+        """
+        return self.with_df.cell.reciprocal_vectors()
+
+    @property
     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
@@ -738,6 +805,12 @@ 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 6ce6f704..84c4c9f2 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 and "B" in self.fsc:
+                    diag_eri[q, ki] = (
+                        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]))
+                    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, None
 
     @logging.with_timer("Density-density moments")
     @logging.with_status("Constructing density-density moments")
@@ -184,7 +200,24 @@ 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(
+                            (np.array([corrected_moments[kj, i][0, :]]), moments[0, kj, i])
+                        )
+
+            return moments
+        else:
+            return moments
+
+    def build_uncorrected_dd_moments(self, integral=None):
         if integral is None:
             integral = self.integrate()
 
@@ -206,8 +239,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
@@ -226,10 +260,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
@@ -239,6 +273,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[ka].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[ka].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
@@ -356,7 +439,6 @@ 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])
@@ -366,6 +448,46 @@ def eval_offset_integral(self, quad, d, Lia=None):
 
         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
@@ -533,3 +655,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
diff --git a/momentGW/pbc/tda.py b/momentGW/pbc/tda.py
index 4983fe00..4a96d7d7 100644
--- a/momentGW/pbc/tda.py
+++ b/momentGW/pbc/tda.py
@@ -36,6 +36,18 @@ class dTDA(MoldTDA):
         Default value is `None`.
     """
 
+    def __init__(
+        self,
+        gw,
+        nmom_max,
+        integrals,
+        mo_energy=None,
+        mo_occ=None,
+        fsc=False,
+    ):
+        super().__init__(gw, nmom_max, integrals, mo_energy=mo_energy, mo_occ=mo_occ)
+        self.fsc = fsc
+
     @logging.with_timer("Density-density moments")
     @logging.with_status("Constructing density-density moments")
     def build_dd_moments(self):
@@ -47,6 +59,33 @@ def build_dd_moments(self):
             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()
+            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])
+                        )
+
+            return moments
+        else:
+            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
         naux = self.naux
@@ -61,6 +100,7 @@ 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)
                 for kj in kpts.loop(1, mpi=True):
                     kb = kpts.member(kpts.wrap_around(kpts[q] + kpts[kj]))
 
@@ -70,11 +110,7 @@ def build_dd_moments(self):
                     )
                     moments[q, kb, i] += moments[q, kb, i - 1] * d.ravel()[None]
 
-                tmp = np.zeros((naux[q], 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())
+                    tmp += np.dot(moments[q, kb, i - 1], self.integrals.Lia[kj, kb].T.conj())
 
                 tmp = mpi_helper.allreduce(tmp)
                 tmp *= 2.0
@@ -87,6 +123,42 @@ def build_dd_moments(self):
 
         return moments
 
+    def build_corrected_dd_moments(self):
+        """Build the finite size corrections for the moments of the
+        density-density response.
+
+        Returns
+        -------
+        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
+
+        # 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)
+            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]),
+                )
+                moments[kj, i] += moments[kj, i - 1] * d.ravel()
+                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[kj])
+
+        return moments
+
     def kernel(self, exact=False):
         """
         Run the polarizability calculation to compute moments of the
@@ -211,9 +283,11 @@ def build_se_moments(self, moments_dd):
         moments_vir : numpy.ndarray
             Moments of the virtual self-energy at each k-point.
         """
-
         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"
@@ -223,17 +297,18 @@ 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)
 
-        # 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 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]))
@@ -244,14 +319,45 @@ def build_se_moments(self, moments_dd):
                     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 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
 
-    @functools.cached_property
+    @property
+    def naux(self):
+        """Number of auxiliaries."""
+        return self.integrals.naux
+
+    @property
     def nov(self):
         """Get the number of ov states in W."""
         return np.multiply.outer(
@@ -268,3 +374,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/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.
diff --git a/momentGW/util.py b/momentGW/util.py
index 7c0ffaf4..421236b5 100644
--- a/momentGW/util.py
+++ b/momentGW/util.py
@@ -456,7 +456,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
 
diff --git a/tests/test_kgw.py b/tests/test_kgw.py
index 9a1dbaca..683106ca 100644
--- a/tests/test_kgw.py
+++ b/tests/test_kgw.py
@@ -108,6 +108,47 @@ 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)
+        kgw.polarizability = "dtda"
+        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.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
+        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
 
@@ -168,21 +209,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)
-
 
 if __name__ == "__main__":
     print("Running tests for KGW")