Merge pull request #257 from phoebe-team/relaxonsTutorial

Coupled Relaxons tutorial

Author: jcoulter12

doc/sphinx/source/faq.rst

@@ -6,20 +6,23 @@ Frequently asked questions
 
 **Why does my iterative or variational BTE solution not converge, or my relaxons solution have negative eigenvalues?** 
 
-If the iterative solver, or especially the variational solver, does not converge, this could be a sign that your scattering matrix is poorly conditioned. You can check this by running the relaxons solver and noting if there are negative eigenvalues (which indicate the matrix is not positive semi-definite as it should be). This could indicate a number of problems:
+If the iterative solver, or especially the variational solver, does not converge, this could be a sign that your scattering matrix is poorly conditioned. You can check this by running the relaxons solver and noting if there are negative eigenvalues (which indicate the matrix is not positive semi-definite as it should be). 
+
+This could indicate a number of problems:
 
 	* For electron and phonon calcultions, you may have chosen a smearing value which violates the conservation of energy too strongly
 	* For an el-ph calculation, your Wannierization may not be very good, producing noisy el-ph matrix elements, energies, etc. For phonon only calculations, you might have poor force constants. 
 	* Symmetrizing the scattering matrix might help - but if you have many negative eigenvalues, likely this is just going to partially cover up a potentially serious problem. 
 
+For more discussion on this topic, see the :ref:`negativeEigenvalues` section of the relaxons tutorial.
 
 **Can I use the rotationally invariant ASR?**
 
 	* One can take advantage of an improved ASR using the matdyn tool of QE. We recommend this for 2D materials. Here, one should use the ASR=‘all’ tag in matdyn. This should output the force constants using rotational invariance + vanishing stress conditions. Then, when Phoebe is run, use sumRuleFC2 = "no", as applying a further ASR will just add additional noise. 
 	https://www.quantum-espresso.org/Doc/INPUT_MATDYN.html#idm17
 
 
-** I have encountered the error `Error in routine phoebe (1): atom mapping failed`` when running the `phoebe-quantum-espresso` version of `ph.x`. 
+**I have encountered the error `Error in routine phoebe (1): atom mapping failed`` when running the `phoebe-quantum-espresso` version of `ph.x`.**
 
 This is a somewhat rare issue that emerges either because of a general issue with the crystal structure or because of a difference in how symmetries are used in the pw.x and ph.x calculations (often, it seems to appear when fractional translations are found in the sym ops list).
 Typically this can be resolved by a standardized representation of your crystal structure -- generated for example by running `AFLOW <http://www.aflowlib.org/aflow-online/>`_ on your structure, using the button for "structure conversion" and "standard primitive" (or any other equivalent utility.)

doc/sphinx/source/inputFiles.rst

@@ -153,7 +153,7 @@ Phonon BTE Solver
 
 * :ref:`checkNegativeRelaxons`
 
-* :ref:`enforcePositiveSemiDefinite`
+* :ref:`enforceDetailedBalance`
 
 .. raw:: html
 
@@ -256,7 +256,7 @@ Electron BTE Solver
 
 * :ref:`checkNegativeRelaxons`
 
-* :ref:`enforcePositiveSemiDefinite`
+* :ref:`enforceDetailedBalance`
 
 .. raw:: html
 
@@ -332,6 +332,112 @@ EPA Transport
   temperatures = [300.]
   dopings = [1.0e21]
 
+-------------------------------------
+
+Coupled Electron-Phonon BTE Solver
+-----------------------------------
+
+**Functionality:** Build and solve the coupled Boltzmann Transport Equation (BTE). Use and outputs are described in :ref:`tutorialCBTE`.
+
+.. raw:: html
+
+  <h3>Input variables</h3>
+
+
+:ref:`appName` = "coupledTransport"
+
+* :ref:`phFC2FileName`
+
+* :ref:`phFC3FileName`
+
+* :ref:`phonopyDispFileName`
+
+* :ref:`phonopyBORNFileName`
+
+* :ref:`sumRuleFC2`
+
+* :ref:`electronH0Name`
+
+* :ref:`wsVecFileName`
+
+* :ref:`elphFileName`
+
+* :ref:`kMesh`
+
+* :ref:`qMesh`
+
+* :ref:`temperatures`
+
+* :ref:`dopings`
+
+* :ref:`chemicalPotentials`
+
+* :ref:`smearingMethod`
+
+* :ref:`smearingWidth`
+
+* :ref:`phSmearingWidth`
+
+* :ref:`elSmearingWidth`
+
+* :ref:`dimensionality`
+
+* :ref:`thickness`
+
+* :ref:`windowType`
+
+* :ref:`windowEnergyLimit`
+
+* :ref:`windowPopulationLimit`
+
+* :ref:`solverBTE`
+
+* :ref:`scatteringMatrixInMemory`
+
+* :ref:`symmetrizeMatrix`
+
+* :ref:`fermiLevel`
+
+* :ref:`numOccupiedStates`
+
+* :ref:`numRelaxonsEigenvalues`
+
+* :ref:`checkNegativeRelaxons`
+
+* :ref:`enforceDetailedBalance`
+
+.. raw:: html
+
+  <h3>Sample input file</h3>
+
+::
+
+  appName = "coupledTransport"
+  
+  sumRuleFC2 = "crystal"
+  phFC2FileName = "silicon.fc"
+  electronH0Name = "si_tb.dat"
+  elphFileName = "silicon.phoebe.elph.hdf5"
+  phFC3FileName = "FORCE_CONSTANTS_3RD"
+
+  kMesh = [25, 25, 25]
+  qMesh = [5, 5, 5]
+  temperatures = [200.]
+  dopings = [1.e21]
+
+  smearingMethod = "gaussian"
+  elSmearingWidth = 0.005 eV
+  phSmearingWidth = 0.002 eV  
+  windowType = "population"
+  windowPopulationLimit = 1e-3
+  numOccupiedStates = 4
+
+  useSymmetries = false
+  enforceDetailedBalance = true 
+  symmetrizeMatrix = true
+  scatteringMatrixInMemory = true
+  solverBTE = ["relaxons"]
+
 -----------------------------------
 
 Phonon Lifetimes on a Path
@@ -953,6 +1059,29 @@ smearingWidth
 * **Required:** yes (when :ref:`smearingMethod` = "gaussian")
 
 
+.. _elSmearingWidth:
+
+elSmearingWidth
+^^^^^^^^^^^^^
+
+* **Description:** This parameter allows different smearing values to be used if :ref:`smearingMethod` = "gaussian" and one is using both el and ph scattering types, as in the coupled BTE, where this parameter represents the full-width half-maximum of the Gaussian used to approximate the Dirac-delta conserving energy. Example: elSmearingWidth = 0.5 eV
+
+* **Format:** *double+units*
+
+* **Required:** no
+
+
+.. _phSmearingWidth:
+
+phSmearingWidth
+^^^^^^^^^^^^^
+
+* **Description:** This parameter allows different smearing values to be used if :ref:`smearingMethod` = "gaussian" and one is using both el and ph scattering types, as in the coupled BTE, where this parameter represents the full-width half-maximum of the Gaussian used to approximate the Dirac-delta conserving energy. Example: phSmearingWidth = 0.5 eV
+
+* **Format:** *double+units*
+
+* **Required:** no 
+
 .. _solverBTE:
 
 solverBTE
@@ -1024,12 +1153,12 @@ checkNegativeRelaxons
 
 * **Default:** `true`
 
-.. _enforcePositiveSemiDefinite:
+.. _enforceDetailedBalance:
 
-enforcePositiveSemiDefinite
+enforceDetailedBalance
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-* **Description:** When `enforcePositiveSemiDefinite` is true, we apply a diagonal perturbation to the scattering matrix to make it positive definite. Note that that this is a bit of a band-aid -- you should check that your matrix is not too far from correct before applying this (that it has just a few small negative eigenvalues at most!).
+* **Description:** When `enforceDetailedBalance` is set to true, we recompute the diagonal elements of the scattering matrix by summing up the rows of the matrix (enforcing detailed balance). This is especially relevant to the relaxons solution to the BTE, in order to avoid issues with negative eigenvalues arising from numerical noise. 
 
 * **Format:** *bool*
 

doc/sphinx/source/theory/coupledBTE.rst

@@ -11,15 +11,15 @@ Coupled BTE
 | *Coupled electron-phonon hydrodynamics and viscous thermoelectric equations.*
 | `arXiv:2503.07560 <https://arxiv.org/abs/2503.07560>`_
 
-| Here, we outline the high-level perspective of the theory behind this development, but encourage you to study the manuscript before using it. This represents a major development in preparation of the second release of Phoebe.
+| Here, we outline the high-level perspective of the theory behind this development, but encourage you to study the manuscript before using it.
 
 Introduction
 -------------------------------
 When we consider the response of electron or phonon populations to a small electric field or thermal gradient, we typically can focus on calculating the out-of-equilibrium population of only the electron or phonon states, assuming the other particles remain at their equilibrium distribution.
 However, in some cases, the scattering between electrons and phonon means that if the el or ph population goes out-of-equilibrium, it can cause a shift in the other's distribution as well.
 Essentially, this is the **electron-phonon drag effect**, the phenomenon in which out-of-equilibrium electrons induce out-of-equilibrium phonons via electron-phonon scattering, and vice-versa.
 
-In order to predict the phonon drag effect on transport, we have to construct the coupled electron-phonon Boltzmann Transport Equation (epBTE). In Phoebe, we implement the epBTE using a full scattering matrix approach, so that the coupled scattering matrix consists of
+In order to predict the phonon drag effect on transport, we have to construct the coupled electron-phonon Boltzmann Transport Equation (CBTE). In Phoebe, we implement the epBTE using a full scattering matrix approach, so that the coupled scattering matrix consists of
 the standard electron and phonon scattering matrixes, along with off diagonal drag terms, which couple the independent electron and phonon subspaces. Written out,
 
 .. math::

doc/sphinx/source/theory/electronBTE.rst

@@ -21,7 +21,10 @@ The linearized electronic BTE for a system exposed to an applied electric field
 
 where the first term describes the diffusion due to an externally applied electric field :math:`\boldsymbol{E}`, the second the diffusion due to a temperature gradient, and the third term is the linearized scattering operator.
 
-The electron scattering matrix :math:`\Omega_{\lambda,\lambda'}` can be computed as,
+Electron scattering rates
+-----------------------------
+
+The electron scattering matrix :math:`\Omega_{\lambda,\lambda'}` can be computed using electron-phonon scattering as,
 
 .. math::
    \Omega_{\boldsymbol{k}b,\boldsymbol{k}'b'} =&
@@ -59,7 +62,10 @@ This scattering matrix requires us to know the phonon and electron energies, as
 
 The Kronecker delta function on momentum is enforces exactly, while the Dirac delta for conservation of energy is instead approximated with methods as described in the section :ref:`delta_fns`.
 
-Since the scattering matrix :math:`\Omega_{\lambda,\lambda'}` is not inherently symmetric, we perform the transformation, 
+Symmetrization of the scattering matrix
+----------------------------------------
+
+Since the scattering matrix :math:`\Omega_{\lambda,\lambda'}` is not inherently symmetric (which we will need if we want to have a real symmetric matrix, as in the realxons solution), we perform the transformation, 
 
 .. math::
    \tilde{\Omega}_{\lambda \lambda'}

doc/sphinx/source/theory/wannier.rst

@@ -72,11 +72,11 @@ What we need to interpolate instead is:
 First, recall the relation between the wavefunctions in Wannier and Bloch representations,
 
 .. math::
-   \big|m\boldsymbol{R}_e\big> = \sum_{b\boldsymbol{k}} e^{-i\boldsymbol{k}\cdot\boldsymbol{R}_e} U_{mb,\boldsymbol{k}} \big|b\boldsymbol{k}\big>
+   \big|m\boldsymbol{R}_e\big> = \frac{1}{N_k} \sum_{b\boldsymbol{k}} e^{-i\boldsymbol{k}\cdot\boldsymbol{R}_e} U_{mb,\boldsymbol{k}} \big|b\boldsymbol{k}\big>
 
 
 .. math::
-   \big|b\boldsymbol{k}\big> = \frac{1}{N_e} \sum_{m\boldsymbol{R}_e} e^{i\boldsymbol{k}\cdot\boldsymbol{R}_e} U_{bm,\boldsymbol{k}}^\dagger \big|m\boldsymbol{R}_e\big>
+   \big|b\boldsymbol{k}\big> = \sum_{m\boldsymbol{R}_e} e^{i\boldsymbol{k}\cdot\boldsymbol{R}_e} U_{bm,\boldsymbol{k}}^\dagger \big|m\boldsymbol{R}_e\big>
 
 where :math:`N_e` is the number of supercells.
 
@@ -93,8 +93,6 @@ To transform the potential from the reciprocal to the real space representation,
    \frac{1}{N_q}
    \sum_{\boldsymbol{q}\nu} e^{-i\boldsymbol{q}\cdot\boldsymbol{R}_p} [u_{\boldsymbol{q}\kappa}^{\nu}]^{-1} \partial_{\boldsymbol{q}\nu} V(\boldsymbol{r})
 
-
-
 So, we first transform to Wannier space by:
 
 .. math::

doc/sphinx/source/tutorials/bands.rst

@@ -3,7 +3,7 @@
 Basic Band Structure and DoS
 ==============================
 
-Synopsis
+Overview
 --------
 
 In addition to transport and lifetime calculations, Phoebe can also perform basic band structure and density of states calculations. These can be useful when debugging a calculation or checking the quality of a Wannier or Fourier interpolation. The quality of interpolation is often dependent on the number of k or q points used in DFT, so sometimes it can be helpful to run this calculation to ensure that your DFT calculation used a dense enough mesh.

doc/sphinx/source/tutorials/coupledBTE.rst

@@ -0,0 +1,106 @@
+.. _tutorialCBTE:
+
+Coupled BTE Tutorial
+=========================
+
+Overview
+--------
+
+In order to predict thermoelectric transport properties including the phonon drag effect, we have to construct the coupled electron-phonon Boltzmann Transport Equation (epBTE). We do this using a full scattering matrix approach, so that the coupled scattering matrix consists of the standard electron and phonon scattering matrices, along with off diagonal drag terms, which couple the previously independent electron and phonon subspaces. Then, we solve the CBTE using the relaxons solution. 
+
+Before proceeding to the CBTE tutorial, we strongly recommend you first read the :ref:`theoryCBTE` section of the theory documentation, and should also review the :ref:`relaxons` to see how the standard electron and phonon relaxons calculations are used. 
+
+.. note::
+  We note that this is an extreme computational effort -- both in memory and compute requirements, as many scattering rates are required. 
+    
+.. important:: 
+  At this time, we consider it a beta feature -- if you encounter difficulty, please post this on the `discussions page <https://github.com/phoebe-team/phoebe/discussions>`_ of the Phoebe GitHub repository. 
+
+Step 0: Preparing the input files 
+----------------------------------
+
+As in the case of the :ref:`phononElectronTransport`, you will need electron and phonon input files from both the :ref:`elWanTransport` and :ref:`phononTransport` calculations (making sure to use the same crystal structure in both el and ph file generation). 
+These files can be taken directly from the output of the el and ph transport tutorials, which generate files for silicon, as in the example `here <https://github.com/phoebe-team/phoebe/tree/develop/example/Silicon-coupled>`_:
+
++---------------------------+-----------------------------+  
+|**For electrons:**         |        **For phonons:**     |
++===========================+=============================+
+| From Quantum ESPRESSO:    | From phono3py:              | 
+|  - ``*.phoebe.elph.hdf5`` |   - ``fc2.hdf5``            |
+|  - ``*.fc``               |   - ``fc3.hdf5``            |
+|  - ``*_tb.dat``           |   - ``phono3py_disp.yaml``  |
+|                           | Or from shengBTE:           |
+|                           |   - ``*.fc``                |
+|                           |   - ``FORCE_CONSTANTS_3RD`` |
++---------------------------+-----------------------------+  
+
+With these in hand, we can proceed to perform a coupled BTE calculation. For this tutorial, we will perform a calculation of doped silicon, so ``* = si`` or ``silicon``.
+
+Step 1: Running the Coupled BTE calculation 
+-------------------------------------------
+
+This calculation will construct a scattering matrix using electron-phonon, phonon-phonon, phonon-isotope, phonon-electron, and drag term scattering rates. 
+We can set up a coupled BTE calculation using the following example input file, as currently exists in ``phoebe/examples/Silicon-coupled/coupledTransport.in``:: 
+    
+  appName = "coupledTransport"
+  
+  sumRuleFC2 = "crystal"
+  phFC2FileName = "silicon.fc"
+  electronH0Name = "si_tb.dat"
+  elphFileName = "silicon.phoebe.elph.hdf5"
+  phFC3FileName = "FORCE_CONSTANTS_3RD"
+
+  kMesh = [25, 25, 25]
+  qMesh = [5, 5, 5]
+  temperatures = [200.]
+  dopings = [1.e21]
+
+  smearingMethod = "gaussian"
+  elSmearingWidth = 0.005 eV
+  phSmearingWidth = 0.002 eV  
+  windowType = "population"
+  windowPopulationLimit = 1e-3
+  numOccupiedStates = 4
+
+  useSymmetries = false
+  enforceDetailedBalance = true 
+  symmetrizeMatrix = true
+  scatteringMatrixInMemory = true
+  solverBTE = ["relaxons"]
+  
+| Most parameters used here are applicable to all relaxon calculations and are described in the relaxons tutorial. 
+| **Here we address a few points specific to the coupled BTE calculation:**
+
+  - The population window limit and k/q-grids here are *very* coarse. One should increase them and the grids used in the calculation until it is converged. 
+  - We have chosen a ``qMesh`` which is commensurate with our ``kMesh``. This is required in the coupled BTE calculation so that the same electron and phonon states are used in the electron-phonon, phonon-electron, and drag contributions to the scattering matrix. 
+  - As with all relaxons calculations, we here choose Gaussian smearing. However now, we have ``elSmearingWidth`` and ``phSmearingWidth`` listed separately. This is required because electron and phonon energy scales are dramatically different, resulting in different requirements for mesh samplings and as a result differences in smearing values. ``phSmearingWidth`` will apply to phonon-phonon and phonon-isotope scattering, and the ``elSmearingWidth`` applies to electron-phonon, phonon-electron, and drag terms. 
+  - Because the coupled BTE calculation is very sensitive to interpolation error with respect to the quality of the electron-phonon matrix elements, it's very likely that we will need to apply ``enforceDetailedBalance = true`` to enforce that the diagonal of the scattering matrix can be reconstructed from the off-diagonal scattering rates (detailed balance).
+  - Note, as with earlier electron and phonon only relaxons solutions, here the use of symmetries is still a research problem, so we have to have ``useSymmetries = false``.
+
+.. note::
+  A very important part of this calculation is that it's extremely easy to come up with negative eigenvalues. 
+  Please check carefully the appendix section of the relevant paper by Coulter et al. for discussion about the positive semi-definite nature of the coupled scattering matrix.
+
+This should be run just as in the other tutorial::
+
+  export OMP_NUM_THREADS=4
+  mpirun -np 4 /path/to/phoebe/build/phoebe -in coupledTransport.in > coupledTransport.out
+  
+We note that this one can take a bit of time, so we recommend using a cluster or workstation for the demonstration. 
+The parallelism of this calculation follows the same principles as that of the standard relaxons solution as discussed in :ref:`relaxonsParallelism`.
+
+Output and Post-Processing 
+--------------------------
+
+As in the standard relaxons solution to the BTE, we have outputs containing bulk transport coefficients, including:
+  
+  * ``relaxons_coupled_viscosity.json``
+  * ``relaxons_coupled_transport_coefficients.json``
+  * ``relaxons_coupled_relaxation_times.json``
+  
+These files contain the viscosity, electrical and thermal conductivity, and relaxation times calculated from the coupled relaxons solution in ``json`` format. 
+
+As before, we will also have the files ``relaxons_el_eigenvectors.hdf5`` or ``relaxons_ph_eigenvectors.hdf5``, which contain the el and ph contributions to each coupled relaxon eigenvector of the scattering matrix. 
+These again can be plotted on the Wigner-Seitz cell using the script provided in ``phoebe/plotScripts/relaxons_eigenvectors.py``.
+
+Finally, the contents of ``relaxons_coupled_real_space_coefficients.json`` can be used to parameterize the Viscous Thermoelectric Equations as in Coulter, Rajkov, and Simoncelli (2025). `arXiv:2503.07560 <https://arxiv.org/abs/2503.07560>`_