diff --git a/docs/source/images/epsilon_vs_q_hBN.jpg b/docs/source/images/epsilon_vs_q_hBN.jpg
new file mode 100644
index 0000000..505723b
Binary files /dev/null and b/docs/source/images/epsilon_vs_q_hBN.jpg differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 76b8a9a..50b026a 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -44,6 +44,7 @@ Documentation Contents
    outputs/conductivity
    outputs/oscillator_strengths
    outputs/spin
+   outputs/invpesilon
 
 .. toctree::
    :maxdepth: 1
diff --git a/docs/source/input_files.rst b/docs/source/input_files.rst
index 30124f9..eadb3bd 100644
--- a/docs/source/input_files.rst
+++ b/docs/source/input_files.rst
@@ -25,21 +25,21 @@ to determine the dimensionality of the system. The expected format is one vector
 
 **# Filling:** Total number of electrons in the unit cell. Required to identify the Fermi level, which is the reference point in the construction of the excitons. Must be an integer number.
 
-**# BravaisVectors:** List of Bravais vectors :math:`\mathbf{R}` that participate in the construction of the Bloch Hamiltonian. Expected one per line, in format ``x y z``.
+**# BravaisVectors:** List of Bravais vectors :math:`\bm{R}` that participate in the construction of the Bloch Hamiltonian. Expected one per line, in format ``x y z``.
 
-**# FockMatrices:** Matrices :math:`H(\mathbf{R})` that construct the Bloch Hamiltonian :math:`H(\mathbf{k})` . The matrices must
+**# FockMatrices:** Matrices :math:`H(\bm{R})` that construct the Bloch Hamiltonian :math:`H(\bm{k})` . The matrices must
 be fully defined, i.e., they cannot be triangular, since the code does not use hermiticity to generate the Bloch Hamiltonian. The Fock matrices given must follow the ordering given in the block `BravaisVectors`. The matrices can be real or complex, and each one must be separated from the next using the delimiter `&`. In case the matrices are complex, the real and imaginary parts must be separated by a space, and the complex part must carry the imaginary umber symbol (e.g. $1.5 −2.1j$ ). Both $i$ and $j$ can be used.
 
 Optional Blocks
 ---------------
 
-**# [OverlapMatrices]:** In case that the orbitals used are not orthonormal, one can optionally provide the overlap matrices :math:`S(\mathbf{R})`. The overlap in $k$ space is given by:
+**# [OverlapMatrices]:** In case that the orbitals used are not orthonormal, one can optionally provide the overlap matrices :math:`S(\bm{R})`. The overlap in $k$ space is given by:
 
 .. math::
 
    S(\bm{k}) = \sum_{\bm{R}}S(\bm{R})e^{i\bm{k}\cdot\bm{R}}
 
-This is necessary to be able to reproduce the bands, which come from solving the generalized eigenvalue problem :math:`H(\mathbf{k})S(\mathbf{k})\Psi = ES(\mathbf{k})\Psi`. This will be specially necessary if the system was determined using DFT, since in tight-binding we usually assume orthonormality. This block follows the same rules as FockMatrices: each matrix :math:`S(\mathbf{R})` must be separated with the delimiter `&`, and they must follow the order given in `BravaisVectors`.
+This is necessary to be able to reproduce the bands, which come from solving the generalized eigenvalue problem :math:`H(\bm{k})S(\bm{k})\Psi = ES(\bm{k})\Psi`. This will be specially necessary if the system was determined using DFT, since in tight-binding we usually assume orthonormality. This block follows the same rules as FockMatrices: each matrix :math:`S(\bm{R})` must be separated with the delimiter `&`, and they must follow the order given in `BravaisVectors`.
 
 DFT Hamiltonians (CRYSTAL and Wannier90)
 ========================================
@@ -101,9 +101,9 @@ Optional Blocks
 
 **# [ShiftMesh]:** Center submesh at ``kx ky kz`` provided.
 
-**# [TotalMomentum]:** Exciton total center-of-mass momentum :math:`\mathbf{Q}`, expects a vector ``qx qy qz``. Defaults to zero.
+**# [TotalMomentum]:** Exciton total center-of-mass momentum :math:`\bm{Q}`, expects a vector ``qx qy qz``. Defaults to zero.
 
-**# [Reciprocal]:** calculates interaction matrix elements in reciprocal space; It takes an integer argument to specify the number of reciprocal cells to sum over, ``nG``.
+**# [Gcutoff]:** calculates interaction matrix elements in reciprocal space; It takes a real argument to specify the cutoff for the reciprocal lattice vectors. Can be smaller than the one specified in the screening file.
 
 **# [Potential]:** Specify the potential function used in the direct term of the kernel of the BSE. `keldysh` or `coulomb` (defaults to `keldysh`)
 
@@ -115,6 +115,47 @@ Optional Blocks
 
 **# [Regularization]:** Set the regularization distance used in the real-space method to avoid the electrostatic divergence at $r = 0$ by setting $V (0) = V (a)$, where a is the regularization distance. By default this parameter is set to the unit cell lattice parameter. It is advised to be changed only for supercell calculations.
 
+**# [Percentage]:** Set the regularization distance :math:`q_0` used in the reciprocal-space method to avoid the electrostatic divergence at $q = 0$ by setting $q_0 = \varsigma k_0$, where $\varsigma$ is the percentage of the smallest wavevector in the BZ mesh. By default this parameter is set to `0.5`. If `0.0` is provided, then the singular term in the potential is set to `0.0`.
+
+Screening File Format
+=====================
+
+This file defines how the microscopic dielectric screening is computed and which parameters to use. It uses the same block-based syntax as the exciton file. Currently,  if a screening file is provided, then an exciton file must be provided as well.
+
+Key Blocks
+----------
+
+**# ncells_aux**: Number of unit cells/kpoints, in each direction, of the auxiliary BZ mesh to compute the polarizability matrix element(s).
+
+**# valence.bands:** Number of valence bands included in the calculations.
+
+**# conduction.bands:** Number of conduction bands included in the calculations.
+
+**# spin:** Indicates whether if the spin degree of freedom has been considered in the system model or not (`true` or `false`).
+
+**# gcutoff:** Defines the cutoff ``Gcutoff`` for the reciprocal lattice vectors to be included in the calculation of the dielectric matrix. Defaults to `2.5`. Can be different (higher) than the one specified in the exciton file.
+
+**# function:** Specifies which functionality of the screening the user intends. It can be `dielectric`, `polarizability`, `inversedielectric` or `exciton`.
+
+.. hlist::
+   :columns: 1
+
+   * **dielectric** Computes the dielectric function at a specified momentum and for a chosen matrix element. The values of the polarizability as a function of included valence and conduction bands are computed and printed in a file named `polarizability_convergence.dat`. Each line of the file has the form ``[# valence bands] [# conduction bands] [Re{χ}] [Im{χ}]``.
+   * **polarizability** Computes the polarizability matrix element :math:`\chi_{\bm{G}\bm{G}'}` for the specified pair :math:`\bm{G}`, :math:`\bm{G}'`  in the BZ mesh. The values are printed to the `polarizability_mesh.dat` file, and each line of the file has the structure ``[kx] [ky] [kz] [Re{χ}] [Im{χ}]``.  Currently works only for 2D materials, whence `kz = 0` at all times.
+   * **inversedielectric** Computes the inverse of the dielectric matrix :math:`\epsilon^{-1}(\bm{q})` at the specified momentum. The matrix elements are printed to the file `<label>_invepsilon.dat`. The order of the `G` vectors is the same as the one they appear by when printed to `stdout` when this option is used. Each row in the file has the real and imaginary parts of each matrix element printed separately separated by a space (and conserving the order of the `G` vectors).
+   * **exciton**  Computes the inverse of the dielectric matrix in the BZ mesh, and proceeds with the calculation of the exciton. A file containing all the points in the BZ mesh is generated. The file has a name of the form `kgrid_<ncells>.dat`, where `ncells` is the number of cells in each direction, specified in the exciton file. The object :math:`\epsilon^{-1}(\bm{q})` is printed to the file `<label>_invepsilon.dat`, where `<label>` is the value of the label parameter specified in the exciton file. All the individual matrices :math:`\epsilon^{-1}(\bm{q})` are printed by the same order of momentum points as that of the file containing the `k` points.
+
+**# vectors:** Indicates which pair of reciprocal lattice vectors :math:`\bm{G}`, :math:`\bm{G}'` the polarizability is computed for. Two indices have to be provided as: `<index1> <index2>`. This entry is valid for both functions `dielectric` and `polarizability`. Defaults to `0 0`.
+
+**# momentum:** Specifies which momentum the polarizability, or inverse of the dielectric matrix, is to be computed at. A vector in the form `qx qy qz` has to be given. Defaults to `0.2 0 0`.
+
+Optional Blocks
+---------------
+
+**# [isotropic]:** Indicates whether the system is isotropic (`true` or `false`). Defaults to `false`.
+
+**# [thickness]:** Indicates the thickness of the 2D material. If provided, it enables the Q2D calculation for the **inversedielectric** or **exciton** functions. Defaults to `0.0`.
+
 Absorption File: `kubo_w.in`
 ============================
 
diff --git a/docs/source/methods/BSE.rst b/docs/source/methods/BSE.rst
index e390296..d8d2ab8 100644
--- a/docs/source/methods/BSE.rst
+++ b/docs/source/methods/BSE.rst
@@ -37,6 +37,25 @@ where:
 
 This is the **Tamm-Dancoff approximation (TDA)** form of the BSE.
 
+Screened RPA Coulomb potential
+==============================
+
+If the user chose the `rpa` option for the interaction potential, then in the strict 2D approximation the screened Coulomb potential is also a matrix at each point in the Brillouin zone, whose elements are given by
+
+.. math::
+   W_{\bm{G}\bm{G}'}(\bm{q}) = \sqrt{v_c(\bm{q}+\bm{G})}\, \epsilon^{-1}_{\bm{G}\bm{G}'}(\bm{q}) \sqrt{v_c(\bm{q}+\bm{G}')} \,,
+
+where :math:`\epsilon^{-1}_{\bm{G}\bm{G}'}(\bm{q})` is the inverse RPA dielectric function computed within the strictly 2D approach. If the user provides the thickness :math:`d_{\perp}` of the 2D material in the screening input file, then the quasi-2D mode (Q2D) is enabled, and the screened potential is given by
+
+.. math::
+   \bar{W}_{\bm{G}\bm{G}'}(\bm{q}) = \sqrt{\bar{v}_c(\bm{q}+\bm{G})} \,\bar{\epsilon}^{-1}_{\bm{G}\bm{G}'}(\bm{q}) \sqrt{\bar{v}_c(\bm{q}+\bm{G}')} \,,
+
+where :math:`\bar{\epsilon}^{-1}_{\bm{G}\bm{G}'}(\bm{q})` is the inverse RPA dielectric function computed within the Q2D approach, and :math:`\bar{v}_c` is given by
+
+.. math::
+   \bar{v}_c(\bm{q}) = \int_{-d_{\perp}/2}^{d_{\perp}/2}  \int_{-d_{\perp}/2}^{d_{\perp}/2}  v_c(\bm{q}, z-z') \, \mathrm{d} z \, \mathrm{d} z'\,,
+
+where :math:`v_c(\bm{q},z-z')` is the Coulomb potential in the mixed :math:`(\bm{q},z)`-representation. Formally, it is given by the in-plane Fourier transform of the unscreened real-space-resolved Coulomb potential :math:`V(\bm{r}-\bm{r}')`, while keeping the $z,z'$ variables intact.
 
 .. Interaction Matrix Elements
 .. =============================
diff --git a/docs/source/methods/screening.rst b/docs/source/methods/screening.rst
index a5db98e..dfc85cf 100644
--- a/docs/source/methods/screening.rst
+++ b/docs/source/methods/screening.rst
@@ -9,6 +9,14 @@ Xatu supports two real-space interaction potentials used in the Bethe-Salpeter E
 
 These govern the electron–hole interaction and are used to build the interaction kernel.
 
+Xatu also supports three interaction potentials in reciprocal space to compute the interaction matrix elements:
+
+1. **Bare Coulomb potential**
+2. **Rytova–Keldysh potential**
+3. **Numerical RPA screened potential**
+
+These potentials model the electron–hole interaction in momentum space through the 2D Fourier transform of the real-space potentials: Coulomb, Rytova–Keldysh and the screened potential with discrete translation symmetry.
+
 .. contents::
    :local:
    :depth: 2
@@ -56,4 +64,54 @@ Anisotropic Screening
 
 Xatu supports anisotropic screening in the Rytova–Keldysh model by allowing directional dependence in the screening length. This is implemented by constructing an effective vector :math:`\mathbf{r}_0 = (r_{0}^{x}, r_{0}^{y}, r_{0}^{z})` , and rescaling the coordinates accordingly.
 
-This allows the screening environment to be tuned independently along in-plane and out-of-plane directions -- a generalization that extends beyond conventional isotropic models.
\ No newline at end of file
+This allows the screening environment to be tuned independently along in-plane and out-of-plane directions -- a generalization that extends beyond conventional isotropic models.
+
+Numerical RPA Screened Potential
+=================================
+
+Xatu can compute the microscopic RPA dielectric function 
+
+.. math::
+   \epsilon_{\bm{G}\bm{G}'}(\bm{q}) = \delta_{\bm{G}\bm{G}'} - \sqrt{v_c(\bm{q}+\bm{G})} \chi^0_{\bm{G}\bm{G}'}(\bm{q}) \sqrt{v_c(\bm{q}+\bm{G}')}
+
+in its symmetric form, where :math:`v_c(\bm{q})` is the 2D Fourier transform of the Coulomb potential given by
+
+.. math::
+   v_c(\bm{q}) = \frac{e}{2 \varepsilon_0 |\bm{q}|}\,,
+
+
+and :math:`\chi^0` is the independent-particle polarizability (or the irreducible polarizability within RPA), which for gapped time-reversal symmetric systems is given by
+
+.. math::
+   \chi^0_{\bm{G}\bm{G}'}(\bm{q}) = \frac{2}{A}  \sum_{vc,\bm{k} \sigma} \frac{\langle c,\bm{k}| e^{-i(\bm{q}+\bm{G})\cdot\bm{r}}|v,\bm{k}+\bm{q}\rangle  \langle v,\bm{k}+\bm{q}| e^{i(\bm{q}+\bm{G}')\cdot\bm{r}}|c,\bm{k}\rangle^*}{\epsilon_{v,\bm{k} + \bm{q}} - \epsilon_{c,\bm{k}} }\,.
+
+
+The matrix elements in the sum are computed for Bloch states written within the linear combination of atomic orbitals (LCAO) approximation and in the point-like orbital approximation.
+
+If the thickness :math:`d_{\perp}` is provided in the screening input file, then the dielectric function to be computed is the symmetric quasi-2D one :math:`\bar{\epsilon}_{\bm{G} \bm{G}'} (\bm{q})`.
+
+To illustrate the new functionalities, we show below the macroscopic dielectric function of monolayer hBN computed with our implementation, and explain to the user how they can obtain it themselves.
+
+.. note::
+
+   📄 **Reference Paper:**   
+   `Microscopic screening theory for excitons in two-dimensional materials: A bridge between effective models and ab initio descriptions, arxiv 2603.10966 <https://arxiv.org/abs/2603.10966>`_
+
+.. image:: ../images/epsilon_vs_q_hBN.jpg
+   :width: 75%
+   :align: center
+
+The continuous black curve above was obtained in Phys. Rev. B 92, 245123 (2015), while the dashed black one was obtained using the QEH package for Python from Nano Lett. 2015, 15, 7, 4616-462.
+
+With the script called ``write_screening.cpp`` in the main folder and using the CRYSTAL model file `hBN_base_HSE0.outp` for monolayer hBN with the HSE06 functional, included among the models provided with Xatu, we obtain the red and blue curves with the new screening functionalities.
+The script can be compiled like any other in the main folder. To run it, execute the following command in the terminal:
+
+.. code-block:: bash
+
+   ./write_screening ../examples/material_models/DFT/hBN_base_HSE06.outp ../examples/excitonconfig/hBN_test.txt ../examples/screeningconfig/hBN_DFT_screening.txt <name_of_q_points_file>.dat
+
+similarly to the command-line usage of the `xatu` binary  (but without any flags), and providing the name of the input file containing the list of q-points at the end of the command. Such a file is already provided in the data folder.
+
+The method which reads the input file containing the list of q-points and computes the dielectric matrix accordingly is called ``ExcitonTB::compute_2D_DielectricMatrix(file_name)``, which can be identified in the script.
+To enable the Q2D calculation, a finite value for the thickness parameter must be provided in the screening input file, as mentioned above.
+After concluding the computation, the file containing the inverse dielectric matrix at each q-point is created, analogously to when running the `xatu` binary to compute the exciton using the RPA screened Coulomb potential (see the previous section).
\ No newline at end of file
diff --git a/docs/source/outputs/invpesilon.rst b/docs/source/outputs/invpesilon.rst
new file mode 100644
index 0000000..0774963
--- /dev/null
+++ b/docs/source/outputs/invpesilon.rst
@@ -0,0 +1,71 @@
+================================
+Inverse Dielectric Matrix Output
+================================
+
+Xatu computes the inverse RPA dielectric matrix if the ``-z`` flag enabling the screening functionalities is provided and the user passes either one of the functions `inversedielectric` or `exciton` in the screening input file.
+If the function `exciton` is used, then a file with the format ``kgrid_*.dat`` containing all the momentum vectors in the BZ mesh is created, alongside the file containing the inverse dielectric matrix.
+The file containing all the generated matrix elements has the format ``*_invpesilon.dat``.
+
+Structure of ``*_invepsilon.dat``
+=================================
+
+If the function `inversedielectric` is used, then only a single inverse dielectric matrix is printed with the structure:
+
+ .. code-block:: text
+
+      Re(G0G0)  Im(G0G0)  Re(G0G1)  Im(G0G1)  ...  Re(G0Gn)  Im(G0Gn)
+      .
+      .
+      .
+      Re(G1G0)  Im(G1G0)  Re(G1G1)  Im(G1G1)  ...  Re(G1Gn)  Im(G1Gn)
+
+where *n-1* is the total number of :math:`\bm{G}`-vectors (including the null one). They are sorted by the same order as they are printed when running Xatu.
+
+If the function `exciton` is used, then the inverse dielectric matrix in the entire BZ mesh is printed into the ``*_invpesilon.dat`` file with a table format, with the structure:
+
+.. code-block:: text
+
+      Re(q0G0G0)  Im(q0G0G0)  Re(q0G0G1)  Im(q0G0G1)  ...  Re(q0G0Gn)  Im(q0G0Gn)
+      .
+      .
+      .
+      Re(q0GnG0)  Im(q0GnG0)  Re(q0GnG1)  Im(q0GnG1)  ...  Re(q0GnGn)  Im(q0GnGn)
+      Re(q1G0G0)  Im(q1G0G0)  Re(q1G0G1)  Im(q1G0G1)  ...  Re(q1G0Gn)  Im(q1G0Gn)
+      .
+      .
+      .
+      Re(q1GnG0)  Im(q1GnG0)  Re(q1GnG1)  Im(q1GnG1)  ...  Re(q1GnGn)  Im(q1GnGn)
+      .
+      .
+      .
+      Re(qNG0G0)  Im(qNG0G0)  Re(qNG0G1)  Im(qNG0G1)  ...  Re(qNG0Gn)  Im(qNG0Gn)
+      .
+      .
+      .
+      Re(qNGnG0) Im(qNGnG0)  Re(qNGnG1)  Im(qNGnG1)  ...  Re(qNGnGn)  Im(qNGnGn)
+
+where *N-1* is the total number of :math:`\bm{q}`-points in the BZ mesh.
+The table has *2(n-1)* columns and *n-1* rows for each :math:`\bm{q}`-point.
+The :math:`\bm{G}`-vectors are sorted by the same order as they are printed when running Xatu, and the :math:`\bm{q}`-points are sorted in the same order as they appear in the ``kgrid_*.dat`` file.
+
+The method ``ExcitonTB::readInverseDielectricMatrix(std::string filename)`` provided in the Xatu library can be used to read the inverse dielectric matrix from a file with the name ``filename`` containing a previously computed one in the same format.
+In this way, the user can repeat an exciton calculation with different parameters (e.g., different BSE solver or regularization scheme) without having to recompute the inverse dielectric matrix.
+For an example of its use in a script, see the ``read_screening.cpp`` file in the main folder of the project's repository.
+**Pro tip**: Using the Xatu as an API tool, the user can also use any language of their preference to read the inverse dielectric matrix, invert it to obtain the dielectric matrix and reduce it, and invert it once again if the user wishes to compute the exciton calculation with a smaller ``Gcutoff`` without having to recompute the inverse dielectric matrix.
+
+**Pro tip**: Using the Xatu as an API tool, the user can start with a previously computed inverse dielectric matrix by reading it from a file and augment it with the method ``ExcitonTB::augment_2D_DielectricMatrix(double Gcutoff)`` upon passing a larger ``Gcutoff`` value. In this way, only the missing matrix elements will be computed, without the need to recompute the entire inverse dielectric matrix from scratch. For this, the method ``ExcitonTB::augment_2D_DielectricMatrix(double Gcutoff)`` has to be called after a succesful call of ``ExcitonTB::readInverseDielectricMatrix(std::string filename)``.
+
+
+Structure of ``kgrid_*.dat``
+============================
+
+The file contains a table with one row per k-point and one column per each momentum component:
+
+.. code-block:: text
+
+   qx0 qy0 qz0
+   qx1 qy1 qz1
+   .
+   .
+   .
+   qxN qyN qzN
\ No newline at end of file
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index d4e8c1e..59a9051 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -83,6 +83,9 @@ Diagonalize the Bloch Hamiltonian at k-points specified in a file. Does not comp
 ``-f, --format model | hdf5``  
 Specify format of the system file. Defaults to `model`. Note: HDF5 support requires compilation with `HDF5=1`.
 
+``-z, --screening screeningfile``  
+Enable screening functionalities. An exciton file must be provided alongside the screening file.
+
 Examples
 ========