Lumen wiki - User contributions [en]

Lumen release notes

2026-07-07T13:48:19Z

Wikiadmin: /* Lumen Releases */

This page contains the list of the releases of the lumen code.

Starting from version 2.0, the release notes are available in a dedicated page on the wiki, while the source code is available on gitlab.

== Lumen Releases ==

* Lumen 2.1: [[Lumen 2.1|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.1.0 Version 2.1.0] (06/07/2026)

* Lumen 2.0: [[Lumen 2.0|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.1 Version 2.0.1] (29/01/2026); [https://zenodo.org/records/20121198 Zenodo 20121198]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Version 2.0.0] (04/12/2025);

== Previous versions ==

* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Lumen 2.1

2026-07-06T08:14:26Z

Wikiadmin: /* Fixes */

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3. A Partial sync with branch 5.4 of the yambo code was done.

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. Tutorial: [[SHG from Floquet formulation]]

* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver. Tutorials: [[LiF#Detection of excitons symmetries (requires lumen 2.1 or better)| Excitons symmetries in LiF]]; [[How to analyse excitons#Detection of excitons symmetries (requires lumen 2.1 or better)| Excitons symmetries in 2D-hBN]]

* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the [https://www.lumen-code.org/wiki/index.php?title=Install_Lumen_with_external_chase_libraries CHASE library] as [[BSE solvers overview|BSE solver]]
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.
* Ypp now prints the absolute value of the dipoles

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable
* Fixed missing variable in the ndb.QP_interpolated database

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

Lumen 2.1

2026-06-16T09:51:10Z

2026-06-11T09:51:54Z

Wikiadmin: /* Lumen Releases */

== Lumen Releases ==

* Lumen 2.1 (11/06/2026) [[Lumen 2.1|Release Notes]]
* [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Lumen 2.0] (04/12/2025) [[Lumen 2.0|Release Notes]]
* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Install Lumen on MacOsX (from Yuncheng Mao)

2026-05-29T08:41:42Z

Wikiadmin: Created page with " = Lumen 2.0.1 macOS Compilation Guide = This guide records a working macOS/Homebrew recipe for compiling the Lumen executables with MPI parallel execution enabled. Two build configurations are documented: * '''MPI+SLK+HDF5_IO''': Serial HDF5/NetCDF I/O. Simpler setup using Homebrew libraries directly. MPI execution and ScaLAPACK are still active. * '''MPI+SLK+HDF5_MPI_IO''': Parallel HDF5/NetCDF I/O. Requires building HDF5, NetCDF-C, and NetCDF-Fortran from source with..."

= Lumen 2.0.1 macOS Compilation Guide =
This guide records a working macOS/Homebrew recipe for compiling the Lumen executables with MPI parallel execution enabled. Two build configurations are documented:

* '''MPI+SLK+HDF5_IO''': Serial HDF5/NetCDF I/O. Simpler setup using Homebrew libraries directly. MPI execution and ScaLAPACK are still active.
* '''MPI+SLK+HDF5_MPI_IO''': Parallel HDF5/NetCDF I/O. Requires building HDF5, NetCDF-C, and NetCDF-Fortran from source with MPI support into a local prefix.

The commands below intentionally avoid conda. Run configure, build, and verification commands through the clean environment wrapper shown here so that Homebrew compilers and libraries are used first.

== Verified Platform ==
This recipe was verified on:

* macOS on Apple Silicon (<code>apple@arm</code>)
* Homebrew under <code>/opt/homebrew</code>
* GNU Fortran from Homebrew GCC
* Open MPI from Homebrew

The final verified executable set from either build is:<pre>
bin/a2y
bin/p2y
bin/yambo
bin/yambo_nl
bin/yambo_ph
bin/yambo_rt
bin/yambo_sc
bin/ypp
bin/ypp_nl
bin/ypp_ph
bin/ypp_rt
bin/ypp_sc
</pre><code>p2y</code> is enabled in this configuration. <code>e2y</code> is disabled by <code>config/mk/global/defs.mk</code>.

== Required Homebrew Packages ==
Install the needed tools and libraries:

<source lang="bash"> brew install gcc open-mpi openblas scalapack fftw netcdf netcdf-fortran hdf5 automake libtool wget cmake </source>

The verified versions were:<pre>
automake 1.18.1
cmake 4.x
fftw 3.3.11
gcc 15.2.0_1
hdf5 2.1.1
libtool 2.5.4
netcdf 4.10.0
netcdf-fortran 4.6.2
open-mpi 5.0.9
openblas 0.3.33
scalapack 2.2.3
</pre>Check that the tools resolve to Homebrew, not conda:

<source lang="bash"> which mpicc mpifort mpirun gfortran nc-config nf-config h5fc cmake </source>

Expected paths are under <code>/opt/homebrew/bin</code>.

== Clean Build Environment ==
Use this wrapper for all configure, build, and verification commands:

<source lang="bash"> env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin <command> </source>

'''Why:'''

* <code>env -i</code> removes conda variables and library paths.
* <code>PATH</code> puts Homebrew and system tools in a known order.
* <code>TMPDIR=/private/tmp</code> avoids macOS compiler temp-file failures seen with a stripped environment.

The examples below use the explicit path form that was verified:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin <command> </source>G

== Build Configuration A: MPI+SLK+HDF5_IO (Serial I/O) ==
This is the simpler configuration. Homebrew HDF5 and NetCDF are used directly, but parallel HDF5 I/O is disabled. MPI execution and ScaLAPACK are still active.

Check whether Homebrew NetCDF has parallel I/O support: <source lang="bash"> nc-config --has-parallel4 </source>

If this prints <code>no</code> (the typical Homebrew default), use this configuration.

=== Configure (Serial I/O) ===
Run from the repository root: <source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin ./configure \
CC=mpicc \
FC=gfortran \
F77=gfortran \
MPICC=mpicc \
MPIFC=mpifort \
MPIF77=mpifort \
LDFLAGS=-Wl,-rpath,/opt/homebrew/lib \
--enable-mpi \
--enable-par-linalg \
--disable-hdf5-par-io \
--with-blas-libs="-L/opt/homebrew/opt/openblas/lib -lopenblas" \
--with-lapack-libs="-L/opt/homebrew/opt/openblas/lib -lopenblas" \
--with-scalapack-libs="-L/opt/homebrew/opt/scalapack/lib -lscalapack" \
--with-fft-path=/opt/homebrew/opt/fftw \
--with-netcdf-path=/opt/homebrew/opt/netcdf \
--with-netcdff-path=/opt/homebrew/opt/netcdf-fortran \
--with-hdf5-path=/opt/homebrew/opt/hdf5
</source>

Expected <code>config/report</code> highlights:<pre>
[E] MPI
[E] SCALAPACK
[E] BLAS : OpenBLAS
[E] LAPACK : OpenBLAS
[E] FFT : FFTW v3
[E] NETCDF
[E] NETCDFF
[E] HDF5 : Serial_lib
[-] Parallel I/O
</pre>
----

== Build Configuration B: MPI+SLK+HDF5_MPI_IO (Parallel I/O) ==
This configuration enables parallel HDF5/NetCDF I/O, which allows Lumen to read and write HDF5 files collectively across MPI ranks. It requires building HDF5, NetCDF-C, and NetCDF-Fortran from source because Homebrew's <code>netcdf</code> formula is not compiled with parallel HDF5 support.

The local libraries are installed into <code>./external/parallel-libs/</code> so that Homebrew packages remain untouched.

=== Build Local Parallel HDF5 ===
Download, configure, and install HDF5 with MPI support: <source lang="bash"> PARLIB=$PWD/external/parallel-libs mkdir -p $PARLIB/src cd $PARLIB/src

# HDF5 1.14.5

wget -O hdf5-1.14.5.tar.gz "https://github.com/HDFGroup/hdf5/releases/download/hdf5_1.14.5/hdf5-1.14.5.tar.gz" tar xzf hdf5-1.14.5.tar.gz cd hdf5-1.14.5

/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
./configure --prefix=$PARLIB --enable-fortran --enable-parallel CC=mpicc FC=mpifort F77=mpifort
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
make -j$(sysctl -n hw.ncpu)
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
make install
</source>

Verify parallel support: <source lang="bash"> $PARLIB/bin/h5pcc -showconfig 2>&1 | grep "Parallel HDF5" </source>

Expected: <code>Parallel HDF5: yes</code>

=== Build Local Parallel NetCDF-C ===
<source lang="bash"> cd $PARLIB/src

# NetCDF-C 4.9.2

wget -O netcdf-c-4.9.2.tar.gz "https://github.com/Unidata/netcdf-c/archive/refs/tags/v4.9.2.tar.gz" tar xzf netcdf-c-4.9.2.tar.gz cd netcdf-c-4.9.2

/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake -S . -B build \
-DCMAKE_INSTALL_PREFIX=$PARLIB \
-DCMAKE_C_COMPILER=mpicc \
-DENABLE_PNETCDF=OFF \
-DENABLE_PARALLEL4=ON \
-DENABLE_HDF5=ON \
-DHDF5_DIR=$PARLIB \
-DENABLE_DAP=OFF \
-DENABLE_TESTS=OFF \
-DENABLE_BYTERANGE=OFF \
-DENABLE_DISKLESS=OFF
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake --build build -j$(sysctl -n hw.ncpu)
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake --install build
</source>

Verify parallel support: <source lang="bash"> $PARLIB/bin/nc-config --has-parallel4 </source>

Expected: <code>yes</code>

=== Build Local Parallel NetCDF-Fortran ===
<source lang="bash"> cd $PARLIB/src

# NetCDF-Fortran 4.6.1

wget -O netcdf-fortran-4.6.1.tar.gz "https://github.com/Unidata/netcdf-fortran/archive/refs/tags/v4.6.1.tar.gz" tar xzf netcdf-fortran-4.6.1.tar.gz cd netcdf-fortran-4.6.1

/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake -S . -B build \
-DCMAKE_INSTALL_PREFIX=$PARLIB \
-DCMAKE_C_COMPILER=mpicc \
-DCMAKE_Fortran_COMPILER=mpifort \
-DENABLE_TESTS=OFF \
-DCMAKE_PREFIX_PATH=$PARLIB \
-DnetCDF_INCLUDE_DIR=$PARLIB/include \
-DnetCDF_LIBRARIES=$PARLIB/lib/libnetcdf.dylib
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake --build build -j$(sysctl -n hw.ncpu)
/usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin \
cmake --install build
</source>

The cmake configure output should show <code>Parallel IO: yes</code> and <code>NetCDF4 Parallel IO: yes</code>.

=== Summary of Local Parallel Libraries ===
After the steps above, <code>./external/parallel-libs/</code> contains:
{| class="wikitable"
!Library
!Version
!Parallel
|-
|HDF5
|1.14.5
|yes
|-
|NetCDF-C
|4.9.2
|yes (parallel4)
|-
|NetCDF-Fortran
|4.6.1
|yes
|}
These are self-contained and do not affect Homebrew.

=== Configure (Parallel I/O) ===
Run from the repository root, pointing to the local parallel libraries:

<source lang="bash"> PARLIB=$PWD/external/parallel-libs /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin ./configure \
CC=mpicc \
FC=gfortran \
F77=gfortran \
MPICC=mpicc \
MPIFC=mpifort \
MPIF77=mpifort \
LDFLAGS="-Wl,-rpath,$PARLIB/lib" \
--enable-mpi \
--enable-par-linalg \
--enable-hdf5-par-io \
--with-blas-libs="-L/opt/homebrew/opt/openblas/lib -lopenblas" \
--with-lapack-libs="-L/opt/homebrew/opt/openblas/lib -lopenblas" \
--with-scalapack-libs="-L/opt/homebrew/opt/scalapack/lib -lscalapack" \
--with-fft-path=/opt/homebrew/opt/fftw \
--with-netcdf-path=$PARLIB \
--with-netcdff-path=$PARLIB \
--with-hdf5-path=$PARLIB
</source>

Expected <code>config/report</code> highlights:<pre>
[E] MPI
[E] SCALAPACK
[E] BLAS : OpenBLAS
[E] LAPACK : OpenBLAS
[E] FFT : FFTW v3
[C] NETCDF
[C] NETCDFF
[E] HDF5 : Parallel_lib
[X] Parallel I/O (via HDF5)
</pre>''Note: NetCDF and NetCDF-Fortran are marked <code>[C]</code> (to be compiled) because the build system compiles them internally against the parallel HDF5. The <code>[E]</code> HDF5 points to the local parallel build.''
----

== Build the Internal Libraries ==
After configure (either configuration), build the internal libraries first:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 iotk ldiago devxlib </source>

Build bundled LibXC:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 libxc </source>

Do not force Homebrew LibXC for this Lumen source. Current Homebrew LibXC is too new for the Fortran module symbols expected by this code. The bundled <code>libxc-5.2.3</code> builds and links correctly.

If the internal tarballs are missing and the machine has network access, these targets download the needed archives into <code>lib/archive/</code>:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make ldiago-dl iotk-dl devxlib-dl libxc-dl </source>

Expected archives:<pre>
lib/archive/Ldiago-0.5.1.tar.gz
lib/archive/iotk-y1.2.2.tar.gz
lib/archive/devicexlib-0.8.5.tar.gz
lib/archive/libxc-5.2.3.tar.gz
</pre>

== Build the Executables ==
Build the project executables one target at a time. Do not put multiple project targets in the same make invocation: project flags are selected from <code>MAKECMDGOALS</code>, and the internal library/stamp setup is sensitive to target switching.

Use <code>-j1</code> for project executables. Higher parallelism can race Fortran module availability during target switches; the compiled programs themselves are still MPI-enabled.

During builds on macOS, this warning can appear many times:<pre>
grep: /proc/cpuinfo: No such file or directory
</pre>It is harmless. The build scripts probe Linux CPU information, but the absence of <code>/proc/cpuinfo</code> on macOS does not prevent compilation.

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 yambo </source>

Then build the remaining enabled targets:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 ypp /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 a2y /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 p2y /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 yambo_ph /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 ypp_ph /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 yambo_sc /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 ypp_sc /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 yambo_rt /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 ypp_rt /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 yambo_nl /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -j1 ypp_nl </source>

In this configuration <code>p2y</code> is enabled and <code>e2y</code> is disabled. The expected outputs are:<pre>
bin/a2y
bin/p2y
bin/yambo
bin/yambo_nl
bin/yambo_ph
bin/yambo_rt
bin/yambo_sc
bin/ypp
bin/ypp_nl
bin/ypp_ph
bin/ypp_rt
bin/ypp_sc
</pre>

== Precision Switching for NL Targets ==
The <code>yambo_nl</code> and <code>ypp_nl</code> targets are compiled with <code>-D_DOUBLE</code> (double precision). The build system's object save/restore mechanism handles switching between single and double precision automatically: when you build <code>yambo</code> (single precision) first and then <code>yambo_nl</code> (double precision), the build system saves the single-precision objects and compiles fresh double-precision ones.

If the save/restore mechanism gets confused (e.g., after a failed build or manual cleanup), remove all stamps and library objects for a clean rebuild:

<source lang="bash"> rm -f config/stamps_and_lists/*.stamp find lib/qe_pseudo lib/slatec lib/math77 lib/local -name '*.lock' -delete 2>/dev/null find lib/qe_pseudo lib/slatec lib/math77 lib/local -name 'objects.save' -type d -exec rm -rf {} + 2>/dev/null find lib/qe_pseudo lib/slatec lib/math77 lib/local -name 'DOUBLE_objects.save' -type d -exec rm -rf {} + 2>/dev/null find src -name '*.o' -delete 2>/dev/null find src -name '*.mod' -delete 2>/dev/null </source>

Then rebuild the targets in order. Single-precision targets first, then NL:

<source lang="bash">

# Single precision targets

make -j1 yambo ypp a2y p2y yambo_ph ypp_ph yambo_sc ypp_sc yambo_rt ypp_rt

# Double precision targets

make -j1 yambo_nl ypp_nl </source>

== Stale Stamp Recovery ==
If a previous failed or cleaned build leaves stamps but missing generated libraries/modules, remove only the stale stamp for the failing library and rerun the target. Examples encountered during this build:

<source lang="bash"> rm -f config/stamps_and_lists/lib_Y_tools.a.stamp rm -f config/stamps_and_lists/libqe_pseudo.a.stamp </source>

The symptom for the <code>libqe_pseudo</code> case is usually:<pre>
Fatal Error: Cannot open module file 'qe_pseudo_m.mod' for reading
</pre>This was seen while switching to the <code>ypp</code> and <code>ypp_nl</code> targets after previous project builds. The module had been built in <code>lib/qe_pseudo/*objects.save/</code>, but the stale stamp prevented the build system from recreating <code>include/qe_pseudo_m.mod</code>.

If you prefer not to remove the stamp, move it aside instead: <source lang="bash"> mv config/stamps_and_lists/libqe_pseudo.a.stamp config/stamps_and_lists/libqe_pseudo.a.stamp.stale </source>

Then rerun the failed target with <code>make -j1 <target></code>.

If <code>libiotk.a</code> is missing after a partial clean, rebuild and copy it:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin make -C lib/iotk/iotk/src libiotk.a /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin ranlib lib/iotk/iotk/src/libiotk.a /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin cp lib/iotk/iotk/src/libiotk.a lib/external/gfortran/mpifort/lib/ </source>

Prefer targeted stamp recovery over broad clean commands. The source tarball is not necessarily a git checkout, and <code>make veryclean</code> can assume git metadata.

== Verification ==
Check the files:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin ls -l bin /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin file bin/* </source>

Expected:<pre>
Mach-O 64-bit executable arm64
</pre>Check linked libraries:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin otool -L bin/* </source>

The dynamic libraries should resolve to <code>/opt/homebrew/opt/...</code>, <code>./external/parallel-libs/</code> (for the parallel I/O build), or <code>/usr/lib</code>. There should be no conda, miniforge, or anaconda paths.

This command should print nothing:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin otool -L bin/* | rg -i 'conda|anaconda|miniforge|mambaforge' </source>

Check help output:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin zsh -lc 'for exe in a2y p2y yambo yambo_nl yambo_ph yambo_rt yambo_sc ypp ypp_nl ypp_ph ypp_rt ypp_sc; do ./bin/$exe -h >/tmp/lumen_${exe}_help.out 2>&1; st=$?; echo "$exe $st"; done' </source>

The verified run printed exit status <code>0</code> for every executable.

Check MPI launch:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin mpirun -np 2 ./bin/yambo -h /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin mpirun -np 2 ./bin/yambo_nl -h /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin mpirun -np 2 ./bin/ypp -h /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin mpirun -np 2 ./bin/ypp_nl -h </source>

In restricted sandboxes, Open MPI may fail with socket <code>bind()</code> permission errors. On a normal macOS shell, the two-rank smoke tests should print help from both MPI ranks and exit with status 0.

== Running Later ==
Use the same clean environment for production runs if conda is installed:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin mpirun -np 4 ./bin/yambo_nl </source>

For single-process runs:

<source lang="bash"> /usr/bin/env -i HOME="$HOME" TMPDIR=/private/tmp PATH=/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin ./bin/yambo </source>

Installation

2026-05-29T08:40:16Z

Wikiadmin:

In this page we explain how to install Lumen on different machine with different compiler and configurations.

===General installations===
* [[Install Lumen on Ubuntu/LinuxMint]]
* [[Install Lumen on Ubuntu/LinuxMint with Intel compiler]]
* [[Install Lumen on Ubuntu/LinuxMint with NVfortran compiler]]
* [[Install Lumenon MacOS (ARM)]]
* [[Install Lumen on a Debian machine with some libraries provided via modules]]
* [[Install Lumen on MacOsX (from Yuncheng Mao)]]

===Libraries===
* [[Compile the libraries in Lumen independently]]
* [[Install Lumen with external HDF5 and NetCDF libraries]]
* [[Install Lumen with external chase libraries|Install Lumen with external ChASE libraries]]

===Installations on specific machines===
* [[Install Lumen Claudia cluster in Marseille|Install Lumen on Claudia cluster (Marseille]])
* [[Install Lumen on Mesocentre (Marseille)]]
* [[Install Lumen on IRENE machine (TGCC Joliot Curie)]]
* [[Install Lumen on ISM cluster (Rome)]]
* [[Install Lumen on Leonardo|Install Lumen on Leonardo (Cineca, Bologna)]]
* [[Install Lumen on Tirant|Install Lumen/LetzElPhC/Yambopy on Tirant (Univ. of Valencia Cluster)]]
* [[Install Lumen/LetzElPhC/Yambopy on Lynx(Univ. of Valencia)| Install Lumen on Lynx (Univ. of Valencia Cluster)]]
* [[Install Lumen on Viper (MPCDF, Munich)]]

First steps: a walk through from DFT to optical properties

2026-05-19T13:38:03Z

Wikiadmin: /* Step 0: Download the Files */

In this tutorial you will learn how to calculate optical spectra using Yambo, starting from a DFT calculation and ending with a look at local field effects in the optical response.

We will use a 3D system (bulk hBN) and a 2D system (hBN sheet).

[[File:HBN-bulk-3x3-annotated.png|x200px|Atomic structure of bulk hBN]]
[[File:HBN2.png|x200px|Atomic structure of 2D hBN]]

Before starting, you need to obtain the tarballs for hBN and hBN-2D. See instructions on the [[Tutorials|main tutorials page]].

The steps are the following: 

===Step 0: Download the Files===
You can download inputs file for [https://media.yambo-code.eu/educational/tutorials/files/hBN.tar.gz hBN bulk] and [https://media.yambo-code.eu/educational/tutorials/files/hBN-2D.tar.gz hBN-2D].

===Step 1: DFT calculation of bulk hBN and conversion to Yambo===
Follow the DFT and p2y module on '''[[Bulk material: h-BN|database generation for bulk hBN]]''' and then '''return to this tutorial "First steps..." '''

There is no need to complete the DFT and p2y module on '''[[2D material: h-BN sheet|database generation in 2D hBN]]''', but you may do so if you wish.

===Step 2: Initialization of Yambo databases===
Use the ''SAVE'' folders that are already provided, rather than any ones you may have generated previously.

Follow the module on '''[[Initialization]]''', for both hBN and 2D-hBN, and '''return to this tutorial "First steps..." '''

===Step 3: Yambo's command line interface===
Yambo uses a command-line interface to select tasks, generate input files, and control the runtime behavior.



For Yambo 5.0: Follow the module on '''[[Input file generation and command line options 5.0|Input file generation and command line options (5.0)]]''' for bulk hBN and then return to this tutorial "First steps...".

===Step 4: Optical absorption in hBN===
Follow the module on '''[[Optics at the independent particle level]]''' for bulk hBN and '''return to this tutorial "First steps..."'''

===Step 5: Optical absorption in 2D BN===

To complete the tutorial you should use the preprepared ''SAVE'' folder for 2D BN and follow the module on '''[[Local fields]]'''.

 


{| style="width:100%" border="1"
|style="width:15%; text-align:left"|Prev: [[Tutorials|Tutorials Home]]
|style="width:50%; text-align:center"|Now: [[Tutorials|Tutorials Home]] --> [[First steps: a walk through from DFT to optical properties|First steps]]
|style="width:35%; text-align:right"|Next: If you did everything, choose another tutorial in the menu
|-
|}

Exciton-phonon coupling and luminescence

2026-04-30T08:00:10Z

Wikiadmin: /* Exciton lifetimes */

[[File:Tdgw-phonon-usc-01-1024x829.jpg|thumb|right|400px|(c) Zhenglu Li, University of Southern California]]

In this advanced tutorial, we will calculate exciton-phonon interactions from first principles by interfacing DFPT (for phonon calculations) and BSE (for exciton calculations).

The DFTP calculations are run with Quantum ESPRESSO, while the many-body GW-BSE calculations are run with Yambo. Finally, the exciton-phonon interaction will be obtained by combining and postprocessing the databases computed in the two previous runs. The great advantage of this workflow is that the calculations can be run in the irreducible Brillouin zones both for the electronic momenta ('''k''') and the transfer momenta ('''Q''', '''q''') of excitons and phonons, thus speeding up considerably the jobs while reducing the IO and memory load.

We will first compute the exciton-phonon coupling matrix elements: these are the building blocks needed to construct experimental observables such as phonon-assisted optical spectra (such as luminescence), Raman spectra and exciton lifetimes. We will do this in the case of monolayer MoS2, a 2D system with large spin-orbit interaction.

As an example of application, we will consider the case of phonon-assisted luminescence. We will do this in the case of bulk hBN, a layered indirect insulator with strong electron-phonon coupling.

'''Note:''' this tutorial will be updated when new exc-ph tools become available in Yambopy (including full-python postprocessing, Raman spectra, interpolated lifetimes, etc).

== Requirements ==

This is an advanced topic: we assume that you already know something about the theory<ref name="toyozawa2003" /><ref name="antonius2017" /><ref name="cudazzo2020" /><ref name="paleari2019_PhD" /><ref name="paleari2022" /><ref name="lechifflart2023_PhD" /> and applications<ref name="paleari2019" /><ref name="cannuccia2019" /><ref name="chen2020" /><ref name="zanfrognini2023" /><ref name="lechifflart2023" /><ref name="chan2023" /><ref name="marini2024" /><ref name="murali2025" /> of exciton-phonon physics.

Also, we assume that you already know how to run both a basic '''Yambo''' GW-BSE calculation and a DFPT phonon calculation with '''Quantum ESPRESSO'''.

Besides the QE executables <code>pw.x</code> and <code>ph.x</code>, we also use the yambo phonon-specific executable <code>yambo_ph</code> and the python utility '''Yambopy'''. The auxiliary code '''LetzElPhC''' (executable <code>lelphc</code>) will be used to obtain the electron-phonon matrix elements by reading the same electronic wavefunctions used by Yambo (and stored in the <code>SAVE</code> directory), while also making full use of crystal symmetries. [https://github.com/yambo-code/LetzElPhC LetzElPhC] will be run by Yambopy, but it must nonetheless be installed. Finally, the exciton-phonon properties can be computed either using <code>yambo_ph</code> or using Yambopy itself. Both cases will be covered in this tutorial.

[[File:Workflow scheme.png|800px|center]]

== Step 0: Pseudopotentials, equilibrium structure and convergence ==

In a real calculation, it is important to ensure that both the pseudopotential and the lattice parameters that we are using are compatible and perform well for the electronic excited states and for the lattice vibrations simultaneously. Furthermore, you have to make sure that the wave function cutoff <code>ecutwfc</code> is converged with respect to the DFPT step and not just to the DFT one. This is in addition to the other customary convergence tests for DFT, DFPT, GW and BSE calculations.

This is often the most time-demanding step when starting on a new system.

For the sake of this tutorial, we assume that we have already done all these tests and we are starting the final workflow to get the exciton-phonon properties.

== Step 1: scf calculation ==

First of all, we run a standard scf calculation with <code>pw.x</code> for Yambo. We stick with non-symmorphic symmetries. At the end, we will have the QE <code>save</code> directory.

This is the input <code>mos2.scf</code>:

&control
wf_collect = .true.,
calculation = "scf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Here we are using full relativistic pseudopotentials from the SG-15 database.

We can run it on our machine (for example using 4 MPI tasks) as:

mpirun -np 4 pw.x -inp mos2.scf > scf.out

== Step 2: nscf calculation for Yambo ==

Copy the QE <code>save</code> directory from the scf calculation and run the nscf calculation for any number of empty states, with the correct <code>k</code>-grid we want to use in Yambo. Here we are using a badly underconverged grid of 6x6x1.

This reciprocal-space grid will also match the momentum transfer <code>q</code> grid on which excitons and phonons will be defined!

The electronic wavefunctions computed at this step and stored in the new nscf <code>save</code> directory will be used both by Yambo and by the electron-phonon code: this is important because using different sets of wavefunctions would lead to a phase mismatch issue in the exciton-phonon matrix elements.

The nscf input <code>mos2.nscf</code> is
&control
wf_collect = .true.,
calculation = "nscf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
nbnd = 250
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Again, we run the calculation

mpirun -np 4 pw.x -inp mos2.nscf > nscf.out

== Step 3: dvscf phonon calculation ==

Now we run the phonon calculation.

Copy the <code>save</code> directory from the '''scf''' calculation and run <code>ph.x</code> for a dvscf calculation with a standard <code>q</code>-grid matching the <code>k</code>-grid we wanna use in Yambo.

At the end, we will have the <code>_ph0</code> directory containing the variation of the self-consistent potential, <math>\Delta V_{SCF}(q)</math>, and the <code>*.dyn</code> files with the phonon energies and eigenvectors.

NB: one could further refine the phonon energies by enforcing the acoustic sum rule, including non-analytic long-range contributions, interpolating to finer grids... all of this can be done within Quantum ESPRESSO and will not be covered in this version of the tutorial.

The input is <code>mos2.dvscf</code>:

mos2_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='mos2',
fildvscf = 'mos2-dvscf',
electron_phonon = 'dvscf',
fildyn='mos2.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=1
/

And now we run as

nohup mpirun -np 8 ph.x -inp mos2.dvscf > dvscf.out &

This time we use nohup and more processes because this calculation may take some time. It is a good idea to set <code>recover=.true.</code> as in a real calculation you will easily breach walltime, and in this way you can safely restart.

== Step 4: create Yambo <code>SAVE</code> directory ==

This is just the standard Yambo initialization: run

p2y

and then

yambo

in the '''nscf''' <code>save</code> folder and then move the newly generated <code>SAVE</code> directory to a convenient place.

== Step 5: run a BSE calculation ==

Now we switch from QE to Yambo. Here, we forgo the GW step for simplicity (we can use a scissor operator to open the band gap).

This calculation has a couple of differences with respect to a standard BSE calculation for optical absorption. We can look at the input file <code>bse.in</code>:

# Runlevels
optics # [R OPT] Optics
rim_cut # [R RIM CUT] Coulomb potential
bss # [R BSS] Bethe Salpeter Equation solver
em1s # [R Xs] Static Inverse Dielectric Matrix
bse # [R BSE] Bethe Salpeter Equation.
bsk # [R BSK] Bethe Salpeter Equation kernel
# RIM and cutoff settings
RandQpts=1000000 # [RIM] Number of random q-points in the BZ
RandGvec= 100 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
# Static screening
X_and_IO_CPU= "1 1 4 2 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
Chimod= "hartree" # [X] IP/Hartree/ALDA/LRC/BSfxc
% BndsRnXs
1 | 200 | # [Xs] Polarization function bands
%
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
NGsBlkXs= 8000 mRy # [Xs] Response block size
# BSE
BS_CPU= "4.1.2" # [PARALLEL] CPUs for each role
BS_ROLEs= "k.eh.t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_diago=4 # [PARALLEL] CPUs for matrix diagonalization
BSEmod= "causal" # [BSE] resonant/causal/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 40000 mRy # [BSK] Exchange components
ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 8000 mRy # [BSK] Screened interaction block size
Lkind="full" #[BSE,X] bar(default)/full/tilde
% KfnQP_E
1.000000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BEnRange
0.00000 | 4.00000 | eV # [BSS] Energy range
%
% BDmRange
0.05000 | 0.05000 | eV # [BSS] Damping range
%
BEnSteps= 2000 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
% BSEQptR
 1 | 7 | # [BSK] Transferred momenta range
%
% BSEBands
25 | 28 | # [BSK] Bands range
%
WRbsWF # [BSS] Write to disk excitonic the FWs

This file was generated using the command: <code> yambo -X s -o b -k sex -y d -r</code> 

First of all, we compute the excitons for all the momenta in the irreducible Brillouin zone for our discrete grid via the <code>BSEQptR</code> variable. This will be a '''finite-momentum''' BSE calculation, analogous to the phonon one.

Second, we change the variable <code>Lkind</code> from <code>bar</code> to <code>full</code>. In Yambo, <code>Lkind="bar"</code>, which is the default for optical absorption, means that we are computing the excitonic response function without the long-range component of the exchange interaction. This cannot be used when computing the exciton momentum dependence, where the long-range exchange interaction can play a role, therefore we have to include it with <code>Lkind="full"</code>. This allows for the calculation of the excitonic longitudinal-transverse splitting (in 3D systems) as well.

We can now run the code:

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

At the end of the calculation, we have obtained the <code>ndb.BS_diago_Q*</code> databases inside the directory <code>bse_Lfull</code>. They contain information on the exciton energies and wavefunctions at each momentum. Do not forget to check the report and logs of your calculation in the same directory to make sure that the code is doing what you want.

== Step 6: obtain the electron-phonon matrix elements ==

We have finished the heavy simulations. Now it's time for the postprocessing. The first order of business is the reconstruction of the electron-phonon coupling matrix elements from the dvscf results and the electronic wavefunctions.

In order to do this, we will run the <code>lelphc</code> executable of the '''LetzElPhC''' code. We will run via command line using yambopy, although it will be instructive to have look at the <code>lelphc</code> input files later.

We run in the same directory where the Yambo <code>SAVE</code> is (remember than you can also virtually move it with a symbolic link).

Type:

yambopy l2y

to see the help for the calculation. Have a look at the various flags and their description.
For example, if we want to do a serial run of LetzElPhC for bands from <math>n_i</math> to <math>n_f</math>, we should type:

yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here <math>n_i</math> and <math>n_f</math> are integers representing the initial and final band indices. '''Please note''': in this case, the band index starts from 1 and the the interval of bands read by the code is <math>[n_i,n_f]</math> including the extrema of the interval.

These should coincide with those used for the Bethe-Salpeter kernel, i.e. those specified in the <code>BSEBands</code> variable of the BSE input file (this is not strictly necessary, but certainly efficient since these calculations use a lot of disk space).

For our system, we want to do a parallel calculation with 4 qpools and 2 kpools. In addition, we want to explicitly specify the path of the <code>lelphc</code> executable and avoid automatically deleting the LetzElPhC data. So we type:

yambopy l2y -ph path/of/dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe

At the end, check your directory: you should find the <code>lelphc.in</code> input file that was run, let's inspect it:

# LetzElPhC input for yambo generated by yambopy
nqpool = 2
nkpool = 4
start_bnd = 25
end_bnd = 28
save_dir = ./SAVE
kernel = dfpt
ph_save_dir = dvscf/ph_save
convention = yambo

Notice the variable <code>convention=yambo</code>: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band <math>n</math> and momentum <math>k-q</math> to band <math>m</math> and momentum <math>k</math>. In the "forward" momentum transfer convention (the more standard one), the transitions go from <math>nk</math> to <math>mk+q</math>. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as <math>\langle mk|dV|nk-q\rangle</math>. This will have consequences also in the formulation of the ''exciton''-phonon coupling matrix element.

You will also find a <code>ndb.elph</code> database that contains the computed electron-phonon matrix elements, in the output format of <code>LetzElPhC</code>.

In addition, if you check the <code>SAVE</code> folder:

ls SAVE/ndb.elph_gkkp*

you will see that yambopy has created the Yambo-compatible electron-phonon databases.

If you want to run LetzElPhC directly, without using yambopy, you can refer to its [https://github.com/yambo-code/LetzElPhC/blob/main/docs/tex_doc/main.pdf|user guide].
Keep in mind that in this case, in order to convert the database to the <code>ndb.elph_gkkp*</code> databases of Yambo, you will then need a couple of lines of python using the Yambopy class <code>ConvertElectronPhononDB</code> in <code>yambopy/letzelph_interface/lelph2y.py</code>.

If you don't need the <code>ndb.elph_gkkp*</code> databases -- e.g., because you are doing the exc-ph postprocessing with yambopy -- you can run <code>yambopy l2y ...</code> with the <code>--no_gkkp</code> option.

If you don't need the <code>ndb.elph</code> databases -- e.g., because you are using <code>yambo_ph</code> for the exc-ph part -- you can run <code>yambopy l2y ...</code> with the <code>--no_lelphc_dbs</code> option.

== Step 7: Obtain the exciton-phonon coupling ==

Now, we can finally access our basic building block for exciton-phonon physics. This could be done entirely in python (using '''Yambopy'''), or by running '''Yambo'''.

* For the '''Yambo postprocessing''' case (running yambo inputs with the <code>yambo_ph</code> executable), '''follow the alternative route to Steps 7-8 [[Exciton-phonon coupling and luminescence - Yambo postprocessing|at this link]]'''.
* For the '''Yambopy postprocessing''' case (using flexible python scripting), '''keep following Steps 7-8 on this page'''.

Our objective is obtaining the following quantity:

<math>\mathcal{G}^\mu_{\alpha\lambda}(0,q)=\sum_{vv^\prime c k} A^{\alpha, *}_{cv^\prime} (k, q) g_{vv^\prime}^\mu (k,q) A^{\lambda}_{cv}(k,0) - \sum_{cc^\prime vk} A^{\alpha, *}_{c^\prime v} (k+q, q) g_{c^\prime c}^\mu (k+q,q) A^{\lambda}_{cv}(k, 0)</math>

Here, <math>A^{\alpha}_{cv}(k,q)</math> are the exciton coefficients extracted from the eigenvector of the two-particles Hamiltonian during the BSE calculation in step 5, while <math>g_{nm}^\mu (k,q)</math> are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton <math>\lambda</math> undergoes phonon-mediated scattering to state <math>\alpha</math> via phonon mode <math>\mu</math>. The scattering can happen for the hole (valence, <math>v</math>) or for the electron (conduction, <math>c</math>).

NB:

(1) This is written in the "backward" momentum transfer convention used by Yambo. The momentum dependence is different in the "forward" transfer convention.

(2) For simplicity, this is written for zero initial exciton momentum. This means that one of the two states involved in the phonon-mediated scattering process will be in the optical limit (and possibly an optically generated exciton), while the other state can have any momentum: this momentum will be the same as the phonon one. This matrix element can be used to describe phonon-assisted absorption and emission spectra.

In order to calculate this quantity using python, we need the <code>ndb.elph</code> databases natively generated by the <code>LetzElPhC</code> code, as well as the BSE databases <code>ndb.BS_diago_Q*</code> containing the information on exciton wavefunctions and energies.

Next, we write a python user script importing the yambopy exciton-phonon tools. You can find a version of this script in <code>yambopy/tutorials/exciton-phonon/calculate_excph.py</code>.

import numpy as np
from yambopy import YamboLatticeDB,YamboWFDB,LetzElphElectronPhononDB #Specific Yambopy dbs loading classes
from yambopy.exciton_phonon.excph_matrix_elements import exciton_phonon_matelem #Specific exc-ph functions

path = '1L_MoS2'
bands_range=[24,28] # 2 valence bands, 2 conduction bands
nexc = 12 # number of excitonic states

bsepath = f'{path}/bse-allq_full' # Path to BSE calculation (Lin=Lout)
savepath = f'{path}/SAVE' # Yambo SAVE
ndb_elph = f'{path}/ndb.elph' # LetzElPhC electron-phonon database (any convention)

# Read lattice
lattice = YamboLatticeDB.from_db_file(filename=f'{savepath}/ns.db1')
# Read electron-phonon
elph = LetzElphElectronPhononDB(ndb_elph,read_all=False)
# Read wave functions
wfcs = YamboWFDB(filename='ns.wf',save=savepath,latdb=lattice,bands_range=bands_range)
# Calculate exciton-phonon matrix elements
exph = exciton_phonon_matelem(lattice,elph,wfcs,BSE_dir=bsepath,nexc_in=nexc,nexc_out=nexc,dmat_mode='save',exph_file='MoS2_Ex-ph.npy')

In this script, we can select the number exciton states <math>\lambda</math> with <code>nexc</code>, the single-particle bands range with <code>bands_range</code>. '''Please note''': in this case, the band index starts from 0 (python indexing) and the the interval of bands read by the code is <math>[n_i,n_f)</math>, therefore the last value (28 in the example) is excluded. That is, in the example we are selecting the 25th, 26th, 27th and 28th bands. Those bands have to be present in both the electron-phonon and BSE calculations.

Here we calculate the couplings of the first twelve states at each finite-<math>q</math> point including <math>q=0</math>. We also include all the nine phonon modes of monolayer MoS2. We also need access to the Yambo <code>SAVE</code> directory in order to read the electronic wavefunctions: they are used to compute the electronic representation matrices (<code>dmat</code> in the code). The argument <code>dmat_mode</code> can be set to <code>'load'</code> for subsequent calculations.

When we are satisfied with the input, we run the code:

python calculate_excph.py

If you check the output, you should find the <code>MoS2_Ex-ph.npy</code> binary file in the directory where you ran.

=== Analysis of the couplings ===

It is a good idea to have a look at what we computed up to now in order to make sure nothing has gone wrong.

It is not easy to know what to expect (apart from symmetry and gauge compliance of the matrix elements), but one can work out the exciton-phonon selection rules in advance, check that the magnitude is reasonable, etc.

It is also not easy to meaningfully plot this quantity. We have to make sure that we are not breaking degenerate states, otherwise the plots will not be invariant.

First of all, we have to know our system: in monolayer MoS2, the first four excitons are all doubly degenerate. The first exciton responsible for a bright peak in the absorption spectrum (the '''A''' peak), is the second state, corresponding to state indices <code>(3,4)</code> in fortran indexing or <code>(2,3)</code> in python indexing.

All these information can be obtained by analyzing the BSE results (this stuff is explained in the BSE tutorials) and by knowledge of the system or class of systems from the literature.

Thus, a good quantity to plot may be the norm of the matrix elements, summed over the degenerate subspace of exciton A, for a certain number of scattered final states mediated by certain phonon modes:

<math>F_A(q)= \sqrt{ \sum_{\alpha \in A,\lambda,\mu} |\mathcal{G}_{\alpha\lambda}^\mu (0,q)|^2 }</math>

In order to do this, we create a python script <code>analyse_excph.py</code> in which we first load the excph dabatases.
You can find a version of this script in the yambopy directory, in <code>tutorials/exciton-phonon</code>.
First, we select the exciton and phonon states to be included in <code>F_A</code>, together with the path of databases and plot details:

# Exciton "in" states
exc_in = [2,3] # First bright peak (A: 2,3 -- B: 6,7)
exc_out = [0,1,2,3] # first 4 states (dispersion of dark triplet state and A)
ph_in = 'all'
# Paths of databases
ns_db1 =f'{path}/SAVE/ns.db1'
ns_ypy = 'MoS2_Ex-ph.npy'

...

Then, we load the data:

# Read lattice and k-space info
ylat = YamboLatticeDB.from_db_file(filename=ns_db1)
print(ylat)

# Load exc-ph database
X_py = np.load(ns_ypy)
G_squared = np.abs(X_py)**2.

The quantity <math>F_A(q)</math> is obtained from a dedicated function as:

if exc_in == 'all': exc_in = range(G_squared.shape[2])
if exc_out == 'all': exc_out = range(G_squared.shape[3])
if ph_in == 'all': ph_in = range(G_squared.shape[1])

G_squared = G_squared[:, ph_in, :, :].sum(axis=(1))
G_squared = G_squared[:, exc_in, :].sum(axis=(1))
G_squared = G_squared[:, exc_out].sum(axis=(1))

F_q = np.sqrt( G_squared )*ha2ev # Switch from Ha to eV

And finally, we have to make a plotting function. For this tutorial we will use a custom scatterplot employing some of the plotting tools provided by yambopy (but you can do whatever you want).

plot_2D_excph(qgrid,G2_to_plot,rlat=ylat.rlat,plt_cbar=True,\
marker='H',s=700,cmap='magma')
...

You can get more experience on using Yambopy for these kinds of visualization by following the [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy Yambopy tutorials]. In fact, remember that these scripts and all the other Yambopy tutorial scripts are just suggestions, not source code written in stone: if you know <code>numpy</code> and <code>matplotlib</code> you can do your own analysis and your own plots, you just need to import the required Yambopy modules to load the data.

In our case, the resulting plot is the following.

[[File:1L MoS2 MoS2 Ex-ph.png|500px|center]]

This can be checked against Fig. 2(d) of reference <ref name="chan2023" />, although you have to keep in mind that our results are badly undersampled in terms of the reciprocal-space grid, as can be easily seen, and the quantity plotted is not exactly the same. However, the main features are already there since they are dictated mostly by crystal symmetries.

Now that we have the exciton-phonon matrix elements, we can use them to build several kinds of observables. Below, we give an example related to phonon-assisted luminescence, but we may update this tutorial in the future to include more cases.

== Step 8: Compute phonon-assisted luminescence ==

[[File:Luminescence scheme.png|250px|right]]

We want to compute the experimental optical signature due to the phonon-assisted recombination of an exciton (as sketched in the figure).

The signal from the phonon replicas can be modeled as a second-order scattering process involving one phonon and one photon:

<math>I^{Sat}(\omega)=\frac{1}{ N_q} \frac{1}{3} \sum_{\epsilon s\beta \mu q} \frac{1}{E_{\beta q}-s\Omega_{\mu q}}\left|\sum_\alpha\frac{ D^{\epsilon}_\alpha \mathcal{G}_{\beta \alpha}^{\mu,*}(q)}{E_\alpha -E_{\beta q} +s\Omega_{\mu q}+\mathrm{i}\eta}\right|^2 \frac{N^{exc}_{\beta q}(T_{exc})[\frac{1+s}{2}+n_{\mu q}(T)]}{\omega -[E_{\beta q}-s\Omega_{\mu q}]+\mathrm{i}\eta}</math>

In this equation, the oscillator strength of the peak is given by the exciton-phonon coupling matrix elements <math>\mathcal{G}</math> multiplied by the exciton dipoles <math>D</math> (they are called "residuals" in Yambo). Here <math>E_\lambda</math> and <math>E_{\alpha q}</math> are the energies of the optical and finite-momentum excitons, respectively, while <math>\Omega_{\mu q}</math> are the phonon energies.

Here, <math>n_{\mu q}(T)</math> is the temperature-dependent phonon Bose-Einstein occupation function. As it can be seen, <math>s=1</math> corresponds to processes of phonon ''emission'' (<math>\propto n(T)+1</math>), while <math>s=-1</math> corresponds to processes of phonon ''absorption'' (<math>\propto n(T)</math>). Therefore, <math>I^{Sat}_{PL}(\omega;T)</math> describes ''light'' emission by recombining excitons mediated by either ''phonon'' absorption or emission.

The quantity <math> N_{\alpha q}(T_{exc})</math> is the exciton occupation function. Luminescence is technically an out-of-equilibrium process, but we can assume that for very low density of excitations and in steady-state conditions, the exciton population can be approximately described by an equilibrium distribution evaluated at an effective temperature. Here, we use the Boltzmann distribution. Experimentally, <math>T_{exc}</math> tends to coincide with the lattice temperature <math>T</math> more or less above 100 K, while at very low temperature (< 10 K), <math>T_{exc}</math> may vary between 10-50 K. It goes without saying that this needs to carefully be checked in your realistic calculations.

Finally, the spectrum is averaged over the polarization directions of the emitted photons (<math>\epsilon=x,y,z</math> representing the respective dipole components).

=== Running the jobs ===

In order to study luminescence in a paradigmatic system, we switch to bulk hexagonal boron nitride and we repeat the workflow. As you can easily see, one can think about automatizing the execution of all these calculations via scripting or more advanced tools. However, in the case of very large simulations (memory-limited or disk-space limited) or for systems whose electronic and lattice properties are fragile with respect to tiny calculation details, one must be very careful and run many basic tests.

Fortunately, we are running a fast underconverged example. We use LDA pseudopotentials from the pseudo-dojo library and the following are the calculations steps.

1. Input <code>hbn.scf</code>:

&control
calculation='scf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR',
outdir = './tmp'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.scf > scf.out

2. Input <code>hbn.nscf</code>:

&control
calculation='nscf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR'
outdir = './'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
nbnd = 120
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.nscf > nscf.out

3. Input <code>hbn.dvscf</code>:

hbn_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='hBN',
fildvscf = 'hBN-dvscf',
electron_phonon = 'dvscf',
fildyn='hBN.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=2
/

nohup mpirun -np 8 pw.x -inp hbn.dvscf > dvscf.out &

4. Input <code>bse.in</code> (we include 2 valence and 2 conduction bands):

optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
em1s
DIP_CPU= "1 8 1" # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_and_IO_CPU= "1 1 1 8 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
BS_CPU= "8 1 1" # [PARALLEL] CPUs for each role
BS_ROLEs= "k eh t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
BS_nCPU_LinAlg_DIAGO=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
% QpntsRXs
1 | 14 | # [Xs] Transferred momenta
%
% BndsRnXs
1 | 120 | # [Xs] Polarization function bands
%
NGsBlkXs= 10 Ry # [Xs] Response block size
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Lkind= "Lfull" # [BSE] Lbar (default) / full
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
% DipBands
1 | 120 | # [DIP] Bands range for dipoles
%
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
DipComputed= "R V P" # [DIP] [default R P V; extra P2 Spin Orb]
BSENGexx= 30000 Ry # [BSK] Exchange components
#ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 9000 Ry # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% KfnQP_E
1.25997 | 1.08816 | 1.12683 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BSEQptR
 1 | 14 | # [BSK] Transferred momenta range
%
% BSEBands
 7 | 10 | # [BSK] Bands range
%
% BEnRange
0.50000 | 8.00000 | eV # [BSS] Energy range
%
% BDmRange
0.050000 | 0.050000 | eV # [BSS] Damping range
%
BEnSteps= 1000 # [BSS] Energy steps
% BLongDir
1.000000 | 1.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
WRbsWF # [BSS] Write to disk excitonic the WFs

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

Importantly, since we want to describe the phonon-assisted recombination process of an *optical* exciton (i.e., emitting a transverse photon), this time we also run an additional calculation at `Q=0` omitting the nonanalytic long-range Coulomb exchange. Make a second input <code>bse_Lbar.in</code> with the following changes:

Lkind= "Lbar" # [BSE] Lbar (default) / full
% BSEQptR
 1 | 1 | # [BSK] Transferred momenta range
%

4b. So now we make a second BSE run in a different directory specified by <code>-J</code>. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases <code>ndb.em1s*</code> that we do not want to recompute from scratch.

mpirun -np 8 yambo -F bse_Lbar.in -J bse_Lbar,bse_Lfull -C bse_Lbar

5. Now we run <code>lelphc</code> (with or without yambopy: in the latter case remember the option <code>-D</code>) to get the el-ph matrix elements, particularly the <code>ndb.elph</code> databases:

yambopy l2y -ph path/of/dvscf/hbn.dvscf -b 7 10 -par 4 2 -D

=== Luminescence calculation ===

6. And finally we calculate exciton-phonon matrix elements and the luminescence spectrum in one go using a python script importing the Yambopy exciton-phonon tools: in order to do this, we need to take a look at all the necessary input variables for the formula written above.

We can start from the script <code>luminescence.py</code> available in <code>tutorials/exciton-phonon</code>.

import numpy as np
import matplotlib.pyplot as plt
from yambopy.exciton_phonon.excph_luminescence import exc_ph_luminescence
from yambopy.exciton_phonon.excph_input_data import exc_ph_get_inputs

We import the necessary tools...

path = '3D_hBN'
# Path to BSE calculation (Lout--> response is Lfull)
bsepath = f'{path}/bse_Lfull' # ndb.BS_diago_Q* databases are needed
# Path to BSE calculation for optically active exciton (Lin --> response is Lbar)
bseBARpath = f'{path}/bse_Lbar' # ndb.BS_diago_Q1 database is needed
# Path to electron-phonon calculation
elphpath = path # ndb.elph is needed
# Path to unprojected dipoles matrix elements (optional)
dipolespath = bsepath # ndb.dipoles is needed (optional)
# Path to lattice and k-space info
savepath = f'{path}/SAVE' # ns.db1 database is needed

... and we specify the paths to the necessary databases. Note that, if we want to perform polarization averaging, we have to recompute the excitonic dipoles in python as seen here.

What about <code>bseBARpath</code>? This variable points to the directory where the databases for the optical (zero-momentum) excitons <math>\alpha</math> (which may be computed with <code>Lkind='Lbar'</code>) is located, which can be different from the directory with the full indirect exciton dispersion <math>\beta</math> (usually computed with <code>Lkind='Lfull'</code>, however <code>Lkind='Ltilde'</code> can also be used if one is interested in "irreducible" excitons). This makes it possible to compute the coupling between different exciton kinds.

bands_range=[6,10] # 2 valence, 2 conduction bands
phonons_range=[0,12] # All phonons
nexc_out = 12 # 12 excitonic states at each momentum (Lout)
nexc_in = 12 # 12 excitonic states at Q=0 (Lin)
T_ph = 10 # Lattice temperature
T_exc = 10 # Effective excitonic temperature

emin=4.4 # Energy range and plot details (in eV)
emax=4.7
estep=0.0002
broad = 0.005 # Broadening parameter for peak width (in eV)

Next, we specify the parameters for the calculation. We include valence bands from 7 to 10 (6 to 9 in python index) and the contribution of all 12 phonon modes. We consider 12 excitonic states for the coupling matrix elements.

Here <code>T_ph</code> and <code>T_exc</code> are the lattice and excitonic temperatures, respectively.

# We calculate and load all the inputs:
# * Exciton-phonon matrix elements
# * Excitonic dipole matrix elements
# * Exciton energies
# * Phonon energies
# * We specify bse_path2=bseBARpath meaning we use Lbar calculation for Q=0 excitons
input_data = exc_ph_get_inputs(savepath,elphpath,bsepath,\
bse_path2=bseBARpath,dipoles_path=dipolespath,\
nexc_in=12,nexc_out=12,\
bands_range=[6,10],phonons_range=[0,12])

ph_energies, exc_energies, exc_energies_in, G, exc_dipoles = input_data

The function <code>exc_ph_get_inputs</code> gives us all the input data for luminescence in the correct format: phonon energies, exciton energies, optical exciton energies (if needed), exciton-phonon matrix elements (calculated on the fly and printed to file), unprojected exciton dipoles (optional, calculated on the fly and printed to file).
The exc-ph matrix element calculation is just a wrapper of the same tools that we have tested in the above section on MoS2.

Finally, there is the calculation of the luminescence spectrum via the function <code>exc_ph_luminescence</code>:
# We calculate the luminescence spectrum including the input data from before
w,PL = exc_ph_luminescence(T_ph,ph_energies,exc_energies,exc_dipoles,G,\
exc_energies_in=exc_energies_in,exc_temp=T_exc,\
nexc_out=nexc_out,nexc_in=nexc_in,emin=emin,emax=emax,\
estep=estep,broad=broad)

If you want to print the luminescence data for later plotting, you can also add the following line to the script:
# Save to file
data = np.column_stack((w, PL))
np.savetxt("hBN_luminescence_12x12x1.dat", data, fmt="%.8f")

It is now time to run the whole script:

python luminescence.py

NB: By using Yambopy for Step 6, we have limited the use of the Yambo code to just step 4 in the entire workflow. This option is more flexible, as it allows for a greater degree of control by the user. On the other hand the Yambo postprocessing route features a Yambo-style input that doesn't require python knowledge and the calculation is currently faster in fortran. However, the luminescence expression computed in Yambo is a slightly different than this one: it is more approximated in the description of the satellite oscillator strengths, but it explicitly includes the renormalization of the direct exciton peak. You can check the differences [[Exciton-phonon coupling and luminescence - Yambo postprocessing|here]].

=== Results ===
We can plot the results of the step 6 calculation. If we do it in the same script we can add something like this:
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_xlim(emin,emax)
ax.set_ylim(0,np.max(PL)*1.1)
ax.get_yaxis().set_visible(False)

ax.plot(w, PL, '-',c='red', label="AA' hBN luminescence")

plt.legend()
plt.savefig('hBN_luminescence.png')
plt.show()

[[File:HBN luminescence satellites.png|500px|center]]

Here, the signal corresponds to a finite-momentum exciton that recombines with the help of several different phonon modes, both optical and acoustic. Each phonon mode whose coupling with the exciton is allowed can generate a peak, and the energy shifts of these peaks with respect to the initial exciton energy correspond to the phonon energies. This result is underconverged, but the main features are all there. In the plot, we show a more converged example using a 12x12x4 grid (all the other parameters being equal). These plots can be compared with Fig. 4(a) of reference <ref name="zanfrognini2023" />.

==Exciton lifetimes==

Using the exciton-phonon coupling it is possible also to calculate exciton life-time due to the scattering with phonon. For this kind of calculation input file can be generated with the command ''yambo_ph -excph l'':

excph # [R] Exction-phonon
ExcGkkp # [R][EXCPH] Exciton-Phonon Matrix Elelements
ExcPhLifeT # [R][EXCPH] Exciton-Phonon Life-Times
BoseTemp= 0.000000 eV # Bosonic Temperature
% ELPhExcStates
1 | 4 | # [EXCPH] Incoming (external) exciton states
%
% ELPhExcSum
1 | 8 | # [EXCPH] Outgoing (virtual) exciton states
%
LoutPath= "none" # [EXCPH] Path of the outgoing L
FANdEtresh= 0.100000E-5 eV # [ELPH] Energy treshold for Fan denominator
% ElPhModes
1 | 12 | # [ELPH] Phonon modes included
%
EXCTemp= 0.000000 eV # [EXCPH] Excitonic Temperature (for luminescence spectra)
AlphaQ= 0.000000 # [EXCPH] Excitonic band structure 2D distortion
#ExcPhOffDiago # [EXCPH] Exciton-Phonon off-diagonal self-energy
#NoMatrxEl # [EXCPH] Set all exciton-phonon matrix elements to one (for testing purpose)
LDamping= 0.1 meV # [EXCPH] Damping of exc-ph self-energy

notice that in this case LDamping is the broadening of the exciton-phonon self energy and it should be a very small value of the order of the phonon life-times around the meV.<ref>[https://arxiv.org/abs/1411.5628 First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong an-harmonicity: The example of palladium hydrides], L. Paulatto, et al. Phys. Rev. B '''91''', 054304 (2015) </ref>. Running this input you will get:

<-> Exciton lifetimes:
<-> Exciton [1] Lifetime: 0.697314 meV
<-> Exciton [2] Lifetime: 0.697314 meV
<-> Exciton [3] Lifetime: 20.22155 meV
<-> Exciton [4] Lifetime: 32.95292 meV

the life-time of the first four excitons. Note that this result is not converged respect to the Q-points. In order to speedup convergence that double-grid<ref name='lechifflart2023'/> can be used also in the life-time calculations, see tutorial on electron-phonon coupling for tricks on the convergence of the double-grid.
-->

== References ==

<references>
<ref name="toyozawa2003" >Toyozawa, Yutaka, and Chris Oxlade, ''Optical processes in solids'', [https://m.booksee.org/book/1121964?force_lang=en Cambridge University Press, (2003)]. </ref>
<ref name='lechifflart2023'>P. Lechifflart, F. Paleari, D. Sangalli, C. Attaccalite, ''First-principles study of luminescence in hexagonal boron nitride single layer: Exciton-phonon coupling and the role of substrate'',
[https://doi.org/10.1103/PhysRevMaterials.7.024006 Phys. Rev. M, '''7''' (2), 024006 (2023)]; [https://arxiv.org/abs/2212.10407 arXiv2212.1047]</ref>
<ref name='cannuccia2019'>E. Cannuccia, B. Monserrat and C. Attaccalite, ''Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride'', [https://doi.org/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]; [https://arxiv.org/abs/1807.11797 arXiv1807.11797]</ref>
<ref name='paleari2019'>F. Paleari et al., ''Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride'', [https://doi.org/10.1103/PhysRevLett.122.187401 Phys. Rev. Lett. '''122''', 187401 (2019)]; [https://arxiv.org/abs/1810.08976 arXiv1810.089776] </ref>
<ref name='chen2020'>''Exciton-Phonon Interaction and Relaxation Times from First Principles'',
Hsiao-Yi Chen, Davide Sangalli, and Marco Bernardi, [https://doi.org/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]; [https://arxiv.org/abs/2002.08913 arXiv 2002.08913 (2020)]</ref>
<ref name="lechifflart2023_PhD">P. Lechifflart, ''Exciton-phonon coupling and phonon-assisted luminescence in hexagonal Boron Nitride nanostructures'', [https://hal.science/tel-04266805v1 PhD Thesis, University of Marseille (2023)]; [https://www.yambo-code.eu/wiki/images/5/54/These_final.pdf From the yambo website]</ref>
<ref name='paleari2019_PhD'>F. Paleari, ''First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride'', [https://wwwen.uni.lu/research/fstm/dphyms/people/fulvio_paleari PhD thesis, University of Luxembourg (2019)]</ref>
<ref name='zanfrognini2023'>M. Zanfrognini et al., ''Distinguishing different stackings in layered materials via luminescence spectroscopy'', [https://doi.org/10.1103/PhysRevLett.131.206902 Phys. Rev. Lett. '''131''', 206902 (2023)]; [https://arxiv.org/abs/2305.17554 arXiv 2305.17554] </ref>
<ref name='marini2024'>G. Marini, M. Calandra, P. Cudazzo, ''Optical absorption and photoluminescence of single layer boron nitride from a first principles cumulant approach'', [https://doi.org/10.1021/acs.nanolett.4c00669 Nano Lett., '''24''', 20, 6017 (2024)]; [https://arxiv.org/abs/2402.03826 arXiv 2402.03826 (2024)]</ref>
<ref name='antonius2017'>G. Antonius, S. G. Louie, ''Theory of exciton-phonon coupling'', [https://doi.org/10.1103/PhysRevB.105.085111 Phys. Rev. B, '''105''', 085111 (2022)]; [https://arxiv.org/abs/1705.04245 arXiv1705.04245 (2017)]</ref>
<ref name='paleari2022'> F. Paleari, and A. Marini, ''Exciton-phonon interaction calls for a revision of the “exciton” concept'', [https://doi.org/10.1103/PhysRevB.106.125403 Phys. Rev. B, '''106''', 125403 (2022)]; [https://arxiv.org/abs/2205.02783 arXiv 2205.02783]</ref>
<ref name='cudazzo2020'> P. Cudazzo, ''First-principles description of the exciton-phonon interaction: A cumulant approach'', [https://doi.org/10.1103/PhysRevB.102.045136 Phys. Rev. B, '''102''', 045136 (2020)]; [https://orbilu.uni.lu/bitstream/10993/44769/1/main.pdf Open access pdf from Luxembourg University]</ref>
<ref name='chan2023'>Y-h Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, ''Exciton Lifetime and Optical Line Width Profile via Exciton–Phonon Interactions: Theory and First-Principles Calculations for Monolayer MoS2'', [https://doi.org/10.1021/acs.nanolett.3c00732 Nano Lett., '''23''', 9 (2023)]; [https://arxiv.org/abs/2212.08451 arXiv 2212.08451 (2023)]</ref>
<ref name='murali2025'>M. Nalabothula, S. Reichardt, L. Wirtz, ''Origin of Interlayer Exciton–Phonon Coupling in 2D Heterostructures'', [https://doi.org/10.1021/acs.nanolett.5c00355 Nano Lett., '''25''', 15 (2025)], [https://arxiv.org/abs/2407.16111 arXiv 2407.16111 (2025)]</ref>
</references>

Exciton-phonon coupling and luminescence

2026-04-30T07:58:35Z

Wikiadmin: /* Exciton lifetimes */

[[File:Tdgw-phonon-usc-01-1024x829.jpg|thumb|right|400px|(c) Zhenglu Li, University of Southern California]]

In this advanced tutorial, we will calculate exciton-phonon interactions from first principles by interfacing DFPT (for phonon calculations) and BSE (for exciton calculations).

The DFTP calculations are run with Quantum ESPRESSO, while the many-body GW-BSE calculations are run with Yambo. Finally, the exciton-phonon interaction will be obtained by combining and postprocessing the databases computed in the two previous runs. The great advantage of this workflow is that the calculations can be run in the irreducible Brillouin zones both for the electronic momenta ('''k''') and the transfer momenta ('''Q''', '''q''') of excitons and phonons, thus speeding up considerably the jobs while reducing the IO and memory load.

We will first compute the exciton-phonon coupling matrix elements: these are the building blocks needed to construct experimental observables such as phonon-assisted optical spectra (such as luminescence), Raman spectra and exciton lifetimes. We will do this in the case of monolayer MoS2, a 2D system with large spin-orbit interaction.

As an example of application, we will consider the case of phonon-assisted luminescence. We will do this in the case of bulk hBN, a layered indirect insulator with strong electron-phonon coupling.

'''Note:''' this tutorial will be updated when new exc-ph tools become available in Yambopy (including full-python postprocessing, Raman spectra, interpolated lifetimes, etc).

== Requirements ==

This is an advanced topic: we assume that you already know something about the theory<ref name="toyozawa2003" /><ref name="antonius2017" /><ref name="cudazzo2020" /><ref name="paleari2019_PhD" /><ref name="paleari2022" /><ref name="lechifflart2023_PhD" /> and applications<ref name="paleari2019" /><ref name="cannuccia2019" /><ref name="chen2020" /><ref name="zanfrognini2023" /><ref name="lechifflart2023" /><ref name="chan2023" /><ref name="marini2024" /><ref name="murali2025" /> of exciton-phonon physics.

Also, we assume that you already know how to run both a basic '''Yambo''' GW-BSE calculation and a DFPT phonon calculation with '''Quantum ESPRESSO'''.

Besides the QE executables <code>pw.x</code> and <code>ph.x</code>, we also use the yambo phonon-specific executable <code>yambo_ph</code> and the python utility '''Yambopy'''. The auxiliary code '''LetzElPhC''' (executable <code>lelphc</code>) will be used to obtain the electron-phonon matrix elements by reading the same electronic wavefunctions used by Yambo (and stored in the <code>SAVE</code> directory), while also making full use of crystal symmetries. [https://github.com/yambo-code/LetzElPhC LetzElPhC] will be run by Yambopy, but it must nonetheless be installed. Finally, the exciton-phonon properties can be computed either using <code>yambo_ph</code> or using Yambopy itself. Both cases will be covered in this tutorial.

[[File:Workflow scheme.png|800px|center]]

== Step 0: Pseudopotentials, equilibrium structure and convergence ==

In a real calculation, it is important to ensure that both the pseudopotential and the lattice parameters that we are using are compatible and perform well for the electronic excited states and for the lattice vibrations simultaneously. Furthermore, you have to make sure that the wave function cutoff <code>ecutwfc</code> is converged with respect to the DFPT step and not just to the DFT one. This is in addition to the other customary convergence tests for DFT, DFPT, GW and BSE calculations.

This is often the most time-demanding step when starting on a new system.

For the sake of this tutorial, we assume that we have already done all these tests and we are starting the final workflow to get the exciton-phonon properties.

== Step 1: scf calculation ==

First of all, we run a standard scf calculation with <code>pw.x</code> for Yambo. We stick with non-symmorphic symmetries. At the end, we will have the QE <code>save</code> directory.

This is the input <code>mos2.scf</code>:

&control
wf_collect = .true.,
calculation = "scf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Here we are using full relativistic pseudopotentials from the SG-15 database.

We can run it on our machine (for example using 4 MPI tasks) as:

mpirun -np 4 pw.x -inp mos2.scf > scf.out

== Step 2: nscf calculation for Yambo ==

Copy the QE <code>save</code> directory from the scf calculation and run the nscf calculation for any number of empty states, with the correct <code>k</code>-grid we want to use in Yambo. Here we are using a badly underconverged grid of 6x6x1.

This reciprocal-space grid will also match the momentum transfer <code>q</code> grid on which excitons and phonons will be defined!

The electronic wavefunctions computed at this step and stored in the new nscf <code>save</code> directory will be used both by Yambo and by the electron-phonon code: this is important because using different sets of wavefunctions would lead to a phase mismatch issue in the exciton-phonon matrix elements.

The nscf input <code>mos2.nscf</code> is
&control
wf_collect = .true.,
calculation = "nscf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
nbnd = 250
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Again, we run the calculation

mpirun -np 4 pw.x -inp mos2.nscf > nscf.out

== Step 3: dvscf phonon calculation ==

Now we run the phonon calculation.

Copy the <code>save</code> directory from the '''scf''' calculation and run <code>ph.x</code> for a dvscf calculation with a standard <code>q</code>-grid matching the <code>k</code>-grid we wanna use in Yambo.

At the end, we will have the <code>_ph0</code> directory containing the variation of the self-consistent potential, <math>\Delta V_{SCF}(q)</math>, and the <code>*.dyn</code> files with the phonon energies and eigenvectors.

NB: one could further refine the phonon energies by enforcing the acoustic sum rule, including non-analytic long-range contributions, interpolating to finer grids... all of this can be done within Quantum ESPRESSO and will not be covered in this version of the tutorial.

The input is <code>mos2.dvscf</code>:

mos2_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='mos2',
fildvscf = 'mos2-dvscf',
electron_phonon = 'dvscf',
fildyn='mos2.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=1
/

And now we run as

nohup mpirun -np 8 ph.x -inp mos2.dvscf > dvscf.out &

This time we use nohup and more processes because this calculation may take some time. It is a good idea to set <code>recover=.true.</code> as in a real calculation you will easily breach walltime, and in this way you can safely restart.

== Step 4: create Yambo <code>SAVE</code> directory ==

This is just the standard Yambo initialization: run

p2y

and then

yambo

in the '''nscf''' <code>save</code> folder and then move the newly generated <code>SAVE</code> directory to a convenient place.

== Step 5: run a BSE calculation ==

Now we switch from QE to Yambo. Here, we forgo the GW step for simplicity (we can use a scissor operator to open the band gap).

This calculation has a couple of differences with respect to a standard BSE calculation for optical absorption. We can look at the input file <code>bse.in</code>:

# Runlevels
optics # [R OPT] Optics
rim_cut # [R RIM CUT] Coulomb potential
bss # [R BSS] Bethe Salpeter Equation solver
em1s # [R Xs] Static Inverse Dielectric Matrix
bse # [R BSE] Bethe Salpeter Equation.
bsk # [R BSK] Bethe Salpeter Equation kernel
# RIM and cutoff settings
RandQpts=1000000 # [RIM] Number of random q-points in the BZ
RandGvec= 100 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
# Static screening
X_and_IO_CPU= "1 1 4 2 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
Chimod= "hartree" # [X] IP/Hartree/ALDA/LRC/BSfxc
% BndsRnXs
1 | 200 | # [Xs] Polarization function bands
%
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
NGsBlkXs= 8000 mRy # [Xs] Response block size
# BSE
BS_CPU= "4.1.2" # [PARALLEL] CPUs for each role
BS_ROLEs= "k.eh.t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_diago=4 # [PARALLEL] CPUs for matrix diagonalization
BSEmod= "causal" # [BSE] resonant/causal/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 40000 mRy # [BSK] Exchange components
ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 8000 mRy # [BSK] Screened interaction block size
Lkind="full" #[BSE,X] bar(default)/full/tilde
% KfnQP_E
1.000000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BEnRange
0.00000 | 4.00000 | eV # [BSS] Energy range
%
% BDmRange
0.05000 | 0.05000 | eV # [BSS] Damping range
%
BEnSteps= 2000 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
% BSEQptR
 1 | 7 | # [BSK] Transferred momenta range
%
% BSEBands
25 | 28 | # [BSK] Bands range
%
WRbsWF # [BSS] Write to disk excitonic the FWs

This file was generated using the command: <code> yambo -X s -o b -k sex -y d -r</code> 

First of all, we compute the excitons for all the momenta in the irreducible Brillouin zone for our discrete grid via the <code>BSEQptR</code> variable. This will be a '''finite-momentum''' BSE calculation, analogous to the phonon one.

Second, we change the variable <code>Lkind</code> from <code>bar</code> to <code>full</code>. In Yambo, <code>Lkind="bar"</code>, which is the default for optical absorption, means that we are computing the excitonic response function without the long-range component of the exchange interaction. This cannot be used when computing the exciton momentum dependence, where the long-range exchange interaction can play a role, therefore we have to include it with <code>Lkind="full"</code>. This allows for the calculation of the excitonic longitudinal-transverse splitting (in 3D systems) as well.

We can now run the code:

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

At the end of the calculation, we have obtained the <code>ndb.BS_diago_Q*</code> databases inside the directory <code>bse_Lfull</code>. They contain information on the exciton energies and wavefunctions at each momentum. Do not forget to check the report and logs of your calculation in the same directory to make sure that the code is doing what you want.

== Step 6: obtain the electron-phonon matrix elements ==

We have finished the heavy simulations. Now it's time for the postprocessing. The first order of business is the reconstruction of the electron-phonon coupling matrix elements from the dvscf results and the electronic wavefunctions.

In order to do this, we will run the <code>lelphc</code> executable of the '''LetzElPhC''' code. We will run via command line using yambopy, although it will be instructive to have look at the <code>lelphc</code> input files later.

We run in the same directory where the Yambo <code>SAVE</code> is (remember than you can also virtually move it with a symbolic link).

Type:

yambopy l2y

to see the help for the calculation. Have a look at the various flags and their description.
For example, if we want to do a serial run of LetzElPhC for bands from <math>n_i</math> to <math>n_f</math>, we should type:

yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here <math>n_i</math> and <math>n_f</math> are integers representing the initial and final band indices. '''Please note''': in this case, the band index starts from 1 and the the interval of bands read by the code is <math>[n_i,n_f]</math> including the extrema of the interval.

These should coincide with those used for the Bethe-Salpeter kernel, i.e. those specified in the <code>BSEBands</code> variable of the BSE input file (this is not strictly necessary, but certainly efficient since these calculations use a lot of disk space).

For our system, we want to do a parallel calculation with 4 qpools and 2 kpools. In addition, we want to explicitly specify the path of the <code>lelphc</code> executable and avoid automatically deleting the LetzElPhC data. So we type:

yambopy l2y -ph path/of/dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe

At the end, check your directory: you should find the <code>lelphc.in</code> input file that was run, let's inspect it:

# LetzElPhC input for yambo generated by yambopy
nqpool = 2
nkpool = 4
start_bnd = 25
end_bnd = 28
save_dir = ./SAVE
kernel = dfpt
ph_save_dir = dvscf/ph_save
convention = yambo

Notice the variable <code>convention=yambo</code>: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band <math>n</math> and momentum <math>k-q</math> to band <math>m</math> and momentum <math>k</math>. In the "forward" momentum transfer convention (the more standard one), the transitions go from <math>nk</math> to <math>mk+q</math>. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as <math>\langle mk|dV|nk-q\rangle</math>. This will have consequences also in the formulation of the ''exciton''-phonon coupling matrix element.

You will also find a <code>ndb.elph</code> database that contains the computed electron-phonon matrix elements, in the output format of <code>LetzElPhC</code>.

In addition, if you check the <code>SAVE</code> folder:

ls SAVE/ndb.elph_gkkp*

you will see that yambopy has created the Yambo-compatible electron-phonon databases.

If you want to run LetzElPhC directly, without using yambopy, you can refer to its [https://github.com/yambo-code/LetzElPhC/blob/main/docs/tex_doc/main.pdf|user guide].
Keep in mind that in this case, in order to convert the database to the <code>ndb.elph_gkkp*</code> databases of Yambo, you will then need a couple of lines of python using the Yambopy class <code>ConvertElectronPhononDB</code> in <code>yambopy/letzelph_interface/lelph2y.py</code>.

If you don't need the <code>ndb.elph_gkkp*</code> databases -- e.g., because you are doing the exc-ph postprocessing with yambopy -- you can run <code>yambopy l2y ...</code> with the <code>--no_gkkp</code> option.

If you don't need the <code>ndb.elph</code> databases -- e.g., because you are using <code>yambo_ph</code> for the exc-ph part -- you can run <code>yambopy l2y ...</code> with the <code>--no_lelphc_dbs</code> option.

== Step 7: Obtain the exciton-phonon coupling ==

Now, we can finally access our basic building block for exciton-phonon physics. This could be done entirely in python (using '''Yambopy'''), or by running '''Yambo'''.

* For the '''Yambo postprocessing''' case (running yambo inputs with the <code>yambo_ph</code> executable), '''follow the alternative route to Steps 7-8 [[Exciton-phonon coupling and luminescence - Yambo postprocessing|at this link]]'''.
* For the '''Yambopy postprocessing''' case (using flexible python scripting), '''keep following Steps 7-8 on this page'''.

Our objective is obtaining the following quantity:

<math>\mathcal{G}^\mu_{\alpha\lambda}(0,q)=\sum_{vv^\prime c k} A^{\alpha, *}_{cv^\prime} (k, q) g_{vv^\prime}^\mu (k,q) A^{\lambda}_{cv}(k,0) - \sum_{cc^\prime vk} A^{\alpha, *}_{c^\prime v} (k+q, q) g_{c^\prime c}^\mu (k+q,q) A^{\lambda}_{cv}(k, 0)</math>

Here, <math>A^{\alpha}_{cv}(k,q)</math> are the exciton coefficients extracted from the eigenvector of the two-particles Hamiltonian during the BSE calculation in step 5, while <math>g_{nm}^\mu (k,q)</math> are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton <math>\lambda</math> undergoes phonon-mediated scattering to state <math>\alpha</math> via phonon mode <math>\mu</math>. The scattering can happen for the hole (valence, <math>v</math>) or for the electron (conduction, <math>c</math>).

NB:

(1) This is written in the "backward" momentum transfer convention used by Yambo. The momentum dependence is different in the "forward" transfer convention.

(2) For simplicity, this is written for zero initial exciton momentum. This means that one of the two states involved in the phonon-mediated scattering process will be in the optical limit (and possibly an optically generated exciton), while the other state can have any momentum: this momentum will be the same as the phonon one. This matrix element can be used to describe phonon-assisted absorption and emission spectra.

In order to calculate this quantity using python, we need the <code>ndb.elph</code> databases natively generated by the <code>LetzElPhC</code> code, as well as the BSE databases <code>ndb.BS_diago_Q*</code> containing the information on exciton wavefunctions and energies.

Next, we write a python user script importing the yambopy exciton-phonon tools. You can find a version of this script in <code>yambopy/tutorials/exciton-phonon/calculate_excph.py</code>.

import numpy as np
from yambopy import YamboLatticeDB,YamboWFDB,LetzElphElectronPhononDB #Specific Yambopy dbs loading classes
from yambopy.exciton_phonon.excph_matrix_elements import exciton_phonon_matelem #Specific exc-ph functions

path = '1L_MoS2'
bands_range=[24,28] # 2 valence bands, 2 conduction bands
nexc = 12 # number of excitonic states

bsepath = f'{path}/bse-allq_full' # Path to BSE calculation (Lin=Lout)
savepath = f'{path}/SAVE' # Yambo SAVE
ndb_elph = f'{path}/ndb.elph' # LetzElPhC electron-phonon database (any convention)

# Read lattice
lattice = YamboLatticeDB.from_db_file(filename=f'{savepath}/ns.db1')
# Read electron-phonon
elph = LetzElphElectronPhononDB(ndb_elph,read_all=False)
# Read wave functions
wfcs = YamboWFDB(filename='ns.wf',save=savepath,latdb=lattice,bands_range=bands_range)
# Calculate exciton-phonon matrix elements
exph = exciton_phonon_matelem(lattice,elph,wfcs,BSE_dir=bsepath,nexc_in=nexc,nexc_out=nexc,dmat_mode='save',exph_file='MoS2_Ex-ph.npy')

In this script, we can select the number exciton states <math>\lambda</math> with <code>nexc</code>, the single-particle bands range with <code>bands_range</code>. '''Please note''': in this case, the band index starts from 0 (python indexing) and the the interval of bands read by the code is <math>[n_i,n_f)</math>, therefore the last value (28 in the example) is excluded. That is, in the example we are selecting the 25th, 26th, 27th and 28th bands. Those bands have to be present in both the electron-phonon and BSE calculations.

Here we calculate the couplings of the first twelve states at each finite-<math>q</math> point including <math>q=0</math>. We also include all the nine phonon modes of monolayer MoS2. We also need access to the Yambo <code>SAVE</code> directory in order to read the electronic wavefunctions: they are used to compute the electronic representation matrices (<code>dmat</code> in the code). The argument <code>dmat_mode</code> can be set to <code>'load'</code> for subsequent calculations.

When we are satisfied with the input, we run the code:

python calculate_excph.py

If you check the output, you should find the <code>MoS2_Ex-ph.npy</code> binary file in the directory where you ran.

=== Analysis of the couplings ===

It is a good idea to have a look at what we computed up to now in order to make sure nothing has gone wrong.

It is not easy to know what to expect (apart from symmetry and gauge compliance of the matrix elements), but one can work out the exciton-phonon selection rules in advance, check that the magnitude is reasonable, etc.

It is also not easy to meaningfully plot this quantity. We have to make sure that we are not breaking degenerate states, otherwise the plots will not be invariant.

First of all, we have to know our system: in monolayer MoS2, the first four excitons are all doubly degenerate. The first exciton responsible for a bright peak in the absorption spectrum (the '''A''' peak), is the second state, corresponding to state indices <code>(3,4)</code> in fortran indexing or <code>(2,3)</code> in python indexing.

All these information can be obtained by analyzing the BSE results (this stuff is explained in the BSE tutorials) and by knowledge of the system or class of systems from the literature.

Thus, a good quantity to plot may be the norm of the matrix elements, summed over the degenerate subspace of exciton A, for a certain number of scattered final states mediated by certain phonon modes:

<math>F_A(q)= \sqrt{ \sum_{\alpha \in A,\lambda,\mu} |\mathcal{G}_{\alpha\lambda}^\mu (0,q)|^2 }</math>

In order to do this, we create a python script <code>analyse_excph.py</code> in which we first load the excph dabatases.
You can find a version of this script in the yambopy directory, in <code>tutorials/exciton-phonon</code>.
First, we select the exciton and phonon states to be included in <code>F_A</code>, together with the path of databases and plot details:

# Exciton "in" states
exc_in = [2,3] # First bright peak (A: 2,3 -- B: 6,7)
exc_out = [0,1,2,3] # first 4 states (dispersion of dark triplet state and A)
ph_in = 'all'
# Paths of databases
ns_db1 =f'{path}/SAVE/ns.db1'
ns_ypy = 'MoS2_Ex-ph.npy'

...

Then, we load the data:

# Read lattice and k-space info
ylat = YamboLatticeDB.from_db_file(filename=ns_db1)
print(ylat)

# Load exc-ph database
X_py = np.load(ns_ypy)
G_squared = np.abs(X_py)**2.

The quantity <math>F_A(q)</math> is obtained from a dedicated function as:

if exc_in == 'all': exc_in = range(G_squared.shape[2])
if exc_out == 'all': exc_out = range(G_squared.shape[3])
if ph_in == 'all': ph_in = range(G_squared.shape[1])

G_squared = G_squared[:, ph_in, :, :].sum(axis=(1))
G_squared = G_squared[:, exc_in, :].sum(axis=(1))
G_squared = G_squared[:, exc_out].sum(axis=(1))

F_q = np.sqrt( G_squared )*ha2ev # Switch from Ha to eV

And finally, we have to make a plotting function. For this tutorial we will use a custom scatterplot employing some of the plotting tools provided by yambopy (but you can do whatever you want).

plot_2D_excph(qgrid,G2_to_plot,rlat=ylat.rlat,plt_cbar=True,\
marker='H',s=700,cmap='magma')
...

You can get more experience on using Yambopy for these kinds of visualization by following the [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy Yambopy tutorials]. In fact, remember that these scripts and all the other Yambopy tutorial scripts are just suggestions, not source code written in stone: if you know <code>numpy</code> and <code>matplotlib</code> you can do your own analysis and your own plots, you just need to import the required Yambopy modules to load the data.

In our case, the resulting plot is the following.

[[File:1L MoS2 MoS2 Ex-ph.png|500px|center]]

This can be checked against Fig. 2(d) of reference <ref name="chan2023" />, although you have to keep in mind that our results are badly undersampled in terms of the reciprocal-space grid, as can be easily seen, and the quantity plotted is not exactly the same. However, the main features are already there since they are dictated mostly by crystal symmetries.

Now that we have the exciton-phonon matrix elements, we can use them to build several kinds of observables. Below, we give an example related to phonon-assisted luminescence, but we may update this tutorial in the future to include more cases.

== Step 8: Compute phonon-assisted luminescence ==

[[File:Luminescence scheme.png|250px|right]]

We want to compute the experimental optical signature due to the phonon-assisted recombination of an exciton (as sketched in the figure).

The signal from the phonon replicas can be modeled as a second-order scattering process involving one phonon and one photon:

<math>I^{Sat}(\omega)=\frac{1}{ N_q} \frac{1}{3} \sum_{\epsilon s\beta \mu q} \frac{1}{E_{\beta q}-s\Omega_{\mu q}}\left|\sum_\alpha\frac{ D^{\epsilon}_\alpha \mathcal{G}_{\beta \alpha}^{\mu,*}(q)}{E_\alpha -E_{\beta q} +s\Omega_{\mu q}+\mathrm{i}\eta}\right|^2 \frac{N^{exc}_{\beta q}(T_{exc})[\frac{1+s}{2}+n_{\mu q}(T)]}{\omega -[E_{\beta q}-s\Omega_{\mu q}]+\mathrm{i}\eta}</math>

In this equation, the oscillator strength of the peak is given by the exciton-phonon coupling matrix elements <math>\mathcal{G}</math> multiplied by the exciton dipoles <math>D</math> (they are called "residuals" in Yambo). Here <math>E_\lambda</math> and <math>E_{\alpha q}</math> are the energies of the optical and finite-momentum excitons, respectively, while <math>\Omega_{\mu q}</math> are the phonon energies.

Here, <math>n_{\mu q}(T)</math> is the temperature-dependent phonon Bose-Einstein occupation function. As it can be seen, <math>s=1</math> corresponds to processes of phonon ''emission'' (<math>\propto n(T)+1</math>), while <math>s=-1</math> corresponds to processes of phonon ''absorption'' (<math>\propto n(T)</math>). Therefore, <math>I^{Sat}_{PL}(\omega;T)</math> describes ''light'' emission by recombining excitons mediated by either ''phonon'' absorption or emission.

The quantity <math> N_{\alpha q}(T_{exc})</math> is the exciton occupation function. Luminescence is technically an out-of-equilibrium process, but we can assume that for very low density of excitations and in steady-state conditions, the exciton population can be approximately described by an equilibrium distribution evaluated at an effective temperature. Here, we use the Boltzmann distribution. Experimentally, <math>T_{exc}</math> tends to coincide with the lattice temperature <math>T</math> more or less above 100 K, while at very low temperature (< 10 K), <math>T_{exc}</math> may vary between 10-50 K. It goes without saying that this needs to carefully be checked in your realistic calculations.

Finally, the spectrum is averaged over the polarization directions of the emitted photons (<math>\epsilon=x,y,z</math> representing the respective dipole components).

=== Running the jobs ===

In order to study luminescence in a paradigmatic system, we switch to bulk hexagonal boron nitride and we repeat the workflow. As you can easily see, one can think about automatizing the execution of all these calculations via scripting or more advanced tools. However, in the case of very large simulations (memory-limited or disk-space limited) or for systems whose electronic and lattice properties are fragile with respect to tiny calculation details, one must be very careful and run many basic tests.

Fortunately, we are running a fast underconverged example. We use LDA pseudopotentials from the pseudo-dojo library and the following are the calculations steps.

1. Input <code>hbn.scf</code>:

&control
calculation='scf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR',
outdir = './tmp'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.scf > scf.out

2. Input <code>hbn.nscf</code>:

&control
calculation='nscf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR'
outdir = './'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
nbnd = 120
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.nscf > nscf.out

3. Input <code>hbn.dvscf</code>:

hbn_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='hBN',
fildvscf = 'hBN-dvscf',
electron_phonon = 'dvscf',
fildyn='hBN.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=2
/

nohup mpirun -np 8 pw.x -inp hbn.dvscf > dvscf.out &

4. Input <code>bse.in</code> (we include 2 valence and 2 conduction bands):

optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
em1s
DIP_CPU= "1 8 1" # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_and_IO_CPU= "1 1 1 8 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
BS_CPU= "8 1 1" # [PARALLEL] CPUs for each role
BS_ROLEs= "k eh t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
BS_nCPU_LinAlg_DIAGO=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
% QpntsRXs
1 | 14 | # [Xs] Transferred momenta
%
% BndsRnXs
1 | 120 | # [Xs] Polarization function bands
%
NGsBlkXs= 10 Ry # [Xs] Response block size
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Lkind= "Lfull" # [BSE] Lbar (default) / full
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
% DipBands
1 | 120 | # [DIP] Bands range for dipoles
%
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
DipComputed= "R V P" # [DIP] [default R P V; extra P2 Spin Orb]
BSENGexx= 30000 Ry # [BSK] Exchange components
#ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 9000 Ry # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% KfnQP_E
1.25997 | 1.08816 | 1.12683 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BSEQptR
 1 | 14 | # [BSK] Transferred momenta range
%
% BSEBands
 7 | 10 | # [BSK] Bands range
%
% BEnRange
0.50000 | 8.00000 | eV # [BSS] Energy range
%
% BDmRange
0.050000 | 0.050000 | eV # [BSS] Damping range
%
BEnSteps= 1000 # [BSS] Energy steps
% BLongDir
1.000000 | 1.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
WRbsWF # [BSS] Write to disk excitonic the WFs

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

Importantly, since we want to describe the phonon-assisted recombination process of an *optical* exciton (i.e., emitting a transverse photon), this time we also run an additional calculation at `Q=0` omitting the nonanalytic long-range Coulomb exchange. Make a second input <code>bse_Lbar.in</code> with the following changes:

Lkind= "Lbar" # [BSE] Lbar (default) / full
% BSEQptR
 1 | 1 | # [BSK] Transferred momenta range
%

4b. So now we make a second BSE run in a different directory specified by <code>-J</code>. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases <code>ndb.em1s*</code> that we do not want to recompute from scratch.

mpirun -np 8 yambo -F bse_Lbar.in -J bse_Lbar,bse_Lfull -C bse_Lbar

5. Now we run <code>lelphc</code> (with or without yambopy: in the latter case remember the option <code>-D</code>) to get the el-ph matrix elements, particularly the <code>ndb.elph</code> databases:

yambopy l2y -ph path/of/dvscf/hbn.dvscf -b 7 10 -par 4 2 -D

=== Luminescence calculation ===

6. And finally we calculate exciton-phonon matrix elements and the luminescence spectrum in one go using a python script importing the Yambopy exciton-phonon tools: in order to do this, we need to take a look at all the necessary input variables for the formula written above.

We can start from the script <code>luminescence.py</code> available in <code>tutorials/exciton-phonon</code>.

import numpy as np
import matplotlib.pyplot as plt
from yambopy.exciton_phonon.excph_luminescence import exc_ph_luminescence
from yambopy.exciton_phonon.excph_input_data import exc_ph_get_inputs

We import the necessary tools...

path = '3D_hBN'
# Path to BSE calculation (Lout--> response is Lfull)
bsepath = f'{path}/bse_Lfull' # ndb.BS_diago_Q* databases are needed
# Path to BSE calculation for optically active exciton (Lin --> response is Lbar)
bseBARpath = f'{path}/bse_Lbar' # ndb.BS_diago_Q1 database is needed
# Path to electron-phonon calculation
elphpath = path # ndb.elph is needed
# Path to unprojected dipoles matrix elements (optional)
dipolespath = bsepath # ndb.dipoles is needed (optional)
# Path to lattice and k-space info
savepath = f'{path}/SAVE' # ns.db1 database is needed

... and we specify the paths to the necessary databases. Note that, if we want to perform polarization averaging, we have to recompute the excitonic dipoles in python as seen here.

What about <code>bseBARpath</code>? This variable points to the directory where the databases for the optical (zero-momentum) excitons <math>\alpha</math> (which may be computed with <code>Lkind='Lbar'</code>) is located, which can be different from the directory with the full indirect exciton dispersion <math>\beta</math> (usually computed with <code>Lkind='Lfull'</code>, however <code>Lkind='Ltilde'</code> can also be used if one is interested in "irreducible" excitons). This makes it possible to compute the coupling between different exciton kinds.

bands_range=[6,10] # 2 valence, 2 conduction bands
phonons_range=[0,12] # All phonons
nexc_out = 12 # 12 excitonic states at each momentum (Lout)
nexc_in = 12 # 12 excitonic states at Q=0 (Lin)
T_ph = 10 # Lattice temperature
T_exc = 10 # Effective excitonic temperature

emin=4.4 # Energy range and plot details (in eV)
emax=4.7
estep=0.0002
broad = 0.005 # Broadening parameter for peak width (in eV)

Next, we specify the parameters for the calculation. We include valence bands from 7 to 10 (6 to 9 in python index) and the contribution of all 12 phonon modes. We consider 12 excitonic states for the coupling matrix elements.

Here <code>T_ph</code> and <code>T_exc</code> are the lattice and excitonic temperatures, respectively.

# We calculate and load all the inputs:
# * Exciton-phonon matrix elements
# * Excitonic dipole matrix elements
# * Exciton energies
# * Phonon energies
# * We specify bse_path2=bseBARpath meaning we use Lbar calculation for Q=0 excitons
input_data = exc_ph_get_inputs(savepath,elphpath,bsepath,\
bse_path2=bseBARpath,dipoles_path=dipolespath,\
nexc_in=12,nexc_out=12,\
bands_range=[6,10],phonons_range=[0,12])

ph_energies, exc_energies, exc_energies_in, G, exc_dipoles = input_data

The function <code>exc_ph_get_inputs</code> gives us all the input data for luminescence in the correct format: phonon energies, exciton energies, optical exciton energies (if needed), exciton-phonon matrix elements (calculated on the fly and printed to file), unprojected exciton dipoles (optional, calculated on the fly and printed to file).
The exc-ph matrix element calculation is just a wrapper of the same tools that we have tested in the above section on MoS2.

Finally, there is the calculation of the luminescence spectrum via the function <code>exc_ph_luminescence</code>:
# We calculate the luminescence spectrum including the input data from before
w,PL = exc_ph_luminescence(T_ph,ph_energies,exc_energies,exc_dipoles,G,\
exc_energies_in=exc_energies_in,exc_temp=T_exc,\
nexc_out=nexc_out,nexc_in=nexc_in,emin=emin,emax=emax,\
estep=estep,broad=broad)

If you want to print the luminescence data for later plotting, you can also add the following line to the script:
# Save to file
data = np.column_stack((w, PL))
np.savetxt("hBN_luminescence_12x12x1.dat", data, fmt="%.8f")

It is now time to run the whole script:

python luminescence.py

NB: By using Yambopy for Step 6, we have limited the use of the Yambo code to just step 4 in the entire workflow. This option is more flexible, as it allows for a greater degree of control by the user. On the other hand the Yambo postprocessing route features a Yambo-style input that doesn't require python knowledge and the calculation is currently faster in fortran. However, the luminescence expression computed in Yambo is a slightly different than this one: it is more approximated in the description of the satellite oscillator strengths, but it explicitly includes the renormalization of the direct exciton peak. You can check the differences [[Exciton-phonon coupling and luminescence - Yambo postprocessing|here]].

=== Results ===
We can plot the results of the step 6 calculation. If we do it in the same script we can add something like this:
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_xlim(emin,emax)
ax.set_ylim(0,np.max(PL)*1.1)
ax.get_yaxis().set_visible(False)

ax.plot(w, PL, '-',c='red', label="AA' hBN luminescence")

plt.legend()
plt.savefig('hBN_luminescence.png')
plt.show()

[[File:HBN luminescence satellites.png|500px|center]]

Here, the signal corresponds to a finite-momentum exciton that recombines with the help of several different phonon modes, both optical and acoustic. Each phonon mode whose coupling with the exciton is allowed can generate a peak, and the energy shifts of these peaks with respect to the initial exciton energy correspond to the phonon energies. This result is underconverged, but the main features are all there. In the plot, we show a more converged example using a 12x12x4 grid (all the other parameters being equal). These plots can be compared with Fig. 4(a) of reference <ref name="zanfrognini2023" />.

==Exciton lifetimes==

Using the exciton-phonon coupling it is possible also to calculate exciton life-time due to the scattering with phonon. For this kind of calculation input file can be generated with the command ''yambo_ph -excph l'':

excph # [R] Exction-phonon
ExcGkkp # [R][EXCPH] Exciton-Phonon Matrix Elelements
ExcPhLifeT # [R][EXCPH] Exciton-Phonon Life-Times
BoseTemp= 0.000000 eV # Bosonic Temperature
% ELPhExcStates
1 | 4 | # [EXCPH] Incoming (external) exciton states
%
% ELPhExcSum
1 | 8 | # [EXCPH] Outgoing (virtual) exciton states
%
LoutPath= "none" # [EXCPH] Path of the outgoing L
FANdEtresh= 0.100000E-5 eV # [ELPH] Energy treshold for Fan denominator
% ElPhModes
1 | 12 | # [ELPH] Phonon modes included
%
EXCTemp= 0.000000 eV # [EXCPH] Excitonic Temperature (for luminescence spectra)
AlphaQ= 0.000000 # [EXCPH] Excitonic band structure 2D distortion
#ExcPhOffDiago # [EXCPH] Exciton-Phonon off-diagonal self-energy
#NoMatrxEl # [EXCPH] Set all exciton-phonon matrix elements to one (for testing purpose)
LDamping= 0.1 meV # [EXCPH] Damping of exc-ph self-energy

notice that in this case LDamping is the broadening of the exciton-phonon self energy and it should be a very small value of the order of the phonon life-times around the meV.<ref>[https://arxiv.org/abs/1411.5628 First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong an-harmonicity: The example of palladium hydrides], L. Paulatto, et al. Phys. Rev. B '''91''', 054304 (2015) </ref>. Running this input you will get:

<-> Exciton lifetimes:
<-> Exciton [1] Lifetime: 0.697314 meV
<-> Exciton [2] Lifetime: 0.697314 meV
<-> Exciton [3] Lifetime: 20.22155 meV
<-> Exciton [4] Lifetime: 32.95292 meV

the life-time of the first four excitons. Note that this result is not converged respect to the Q-points. In order to speedup convergence that double-grid<ref name='pierre' /> can be used also in the life-time calculations, see tutorial on electron-phonon coupling for tricks on the convergence of the double-grid.
-->

== References ==

<references>
<ref name="toyozawa2003" >Toyozawa, Yutaka, and Chris Oxlade, ''Optical processes in solids'', [https://m.booksee.org/book/1121964?force_lang=en Cambridge University Press, (2003)]. </ref>
<ref name='lechifflart2023'>P. Lechifflart, F. Paleari, D. Sangalli, C. Attaccalite, ''First-principles study of luminescence in hexagonal boron nitride single layer: Exciton-phonon coupling and the role of substrate'',
[https://doi.org/10.1103/PhysRevMaterials.7.024006 Phys. Rev. M, '''7''' (2), 024006 (2023)]; [https://arxiv.org/abs/2212.10407 arXiv2212.1047]</ref>
<ref name='cannuccia2019'>E. Cannuccia, B. Monserrat and C. Attaccalite, ''Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride'', [https://doi.org/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]; [https://arxiv.org/abs/1807.11797 arXiv1807.11797]</ref>
<ref name='paleari2019'>F. Paleari et al., ''Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride'', [https://doi.org/10.1103/PhysRevLett.122.187401 Phys. Rev. Lett. '''122''', 187401 (2019)]; [https://arxiv.org/abs/1810.08976 arXiv1810.089776] </ref>
<ref name='chen2020'>''Exciton-Phonon Interaction and Relaxation Times from First Principles'',
Hsiao-Yi Chen, Davide Sangalli, and Marco Bernardi, [https://doi.org/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]; [https://arxiv.org/abs/2002.08913 arXiv 2002.08913 (2020)]</ref>
<ref name="lechifflart2023_PhD">P. Lechifflart, ''Exciton-phonon coupling and phonon-assisted luminescence in hexagonal Boron Nitride nanostructures'', [https://hal.science/tel-04266805v1 PhD Thesis, University of Marseille (2023)]; [https://www.yambo-code.eu/wiki/images/5/54/These_final.pdf From the yambo website]</ref>
<ref name='paleari2019_PhD'>F. Paleari, ''First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride'', [https://wwwen.uni.lu/research/fstm/dphyms/people/fulvio_paleari PhD thesis, University of Luxembourg (2019)]</ref>
<ref name='zanfrognini2023'>M. Zanfrognini et al., ''Distinguishing different stackings in layered materials via luminescence spectroscopy'', [https://doi.org/10.1103/PhysRevLett.131.206902 Phys. Rev. Lett. '''131''', 206902 (2023)]; [https://arxiv.org/abs/2305.17554 arXiv 2305.17554] </ref>
<ref name='marini2024'>G. Marini, M. Calandra, P. Cudazzo, ''Optical absorption and photoluminescence of single layer boron nitride from a first principles cumulant approach'', [https://doi.org/10.1021/acs.nanolett.4c00669 Nano Lett., '''24''', 20, 6017 (2024)]; [https://arxiv.org/abs/2402.03826 arXiv 2402.03826 (2024)]</ref>
<ref name='antonius2017'>G. Antonius, S. G. Louie, ''Theory of exciton-phonon coupling'', [https://doi.org/10.1103/PhysRevB.105.085111 Phys. Rev. B, '''105''', 085111 (2022)]; [https://arxiv.org/abs/1705.04245 arXiv1705.04245 (2017)]</ref>
<ref name='paleari2022'> F. Paleari, and A. Marini, ''Exciton-phonon interaction calls for a revision of the “exciton” concept'', [https://doi.org/10.1103/PhysRevB.106.125403 Phys. Rev. B, '''106''', 125403 (2022)]; [https://arxiv.org/abs/2205.02783 arXiv 2205.02783]</ref>
<ref name='cudazzo2020'> P. Cudazzo, ''First-principles description of the exciton-phonon interaction: A cumulant approach'', [https://doi.org/10.1103/PhysRevB.102.045136 Phys. Rev. B, '''102''', 045136 (2020)]; [https://orbilu.uni.lu/bitstream/10993/44769/1/main.pdf Open access pdf from Luxembourg University]</ref>
<ref name='chan2023'>Y-h Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, ''Exciton Lifetime and Optical Line Width Profile via Exciton–Phonon Interactions: Theory and First-Principles Calculations for Monolayer MoS2'', [https://doi.org/10.1021/acs.nanolett.3c00732 Nano Lett., '''23''', 9 (2023)]; [https://arxiv.org/abs/2212.08451 arXiv 2212.08451 (2023)]</ref>
<ref name='murali2025'>M. Nalabothula, S. Reichardt, L. Wirtz, ''Origin of Interlayer Exciton–Phonon Coupling in 2D Heterostructures'', [https://doi.org/10.1021/acs.nanolett.5c00355 Nano Lett., '''25''', 15 (2025)], [https://arxiv.org/abs/2407.16111 arXiv 2407.16111 (2025)]</ref>
</references>

Exciton-phonon coupling and luminescence

2026-04-30T07:58:17Z

Wikiadmin:

[[File:Tdgw-phonon-usc-01-1024x829.jpg|thumb|right|400px|(c) Zhenglu Li, University of Southern California]]

In this advanced tutorial, we will calculate exciton-phonon interactions from first principles by interfacing DFPT (for phonon calculations) and BSE (for exciton calculations).

The DFTP calculations are run with Quantum ESPRESSO, while the many-body GW-BSE calculations are run with Yambo. Finally, the exciton-phonon interaction will be obtained by combining and postprocessing the databases computed in the two previous runs. The great advantage of this workflow is that the calculations can be run in the irreducible Brillouin zones both for the electronic momenta ('''k''') and the transfer momenta ('''Q''', '''q''') of excitons and phonons, thus speeding up considerably the jobs while reducing the IO and memory load.

We will first compute the exciton-phonon coupling matrix elements: these are the building blocks needed to construct experimental observables such as phonon-assisted optical spectra (such as luminescence), Raman spectra and exciton lifetimes. We will do this in the case of monolayer MoS2, a 2D system with large spin-orbit interaction.

As an example of application, we will consider the case of phonon-assisted luminescence. We will do this in the case of bulk hBN, a layered indirect insulator with strong electron-phonon coupling.

'''Note:''' this tutorial will be updated when new exc-ph tools become available in Yambopy (including full-python postprocessing, Raman spectra, interpolated lifetimes, etc).

== Requirements ==

This is an advanced topic: we assume that you already know something about the theory<ref name="toyozawa2003" /><ref name="antonius2017" /><ref name="cudazzo2020" /><ref name="paleari2019_PhD" /><ref name="paleari2022" /><ref name="lechifflart2023_PhD" /> and applications<ref name="paleari2019" /><ref name="cannuccia2019" /><ref name="chen2020" /><ref name="zanfrognini2023" /><ref name="lechifflart2023" /><ref name="chan2023" /><ref name="marini2024" /><ref name="murali2025" /> of exciton-phonon physics.

Also, we assume that you already know how to run both a basic '''Yambo''' GW-BSE calculation and a DFPT phonon calculation with '''Quantum ESPRESSO'''.

Besides the QE executables <code>pw.x</code> and <code>ph.x</code>, we also use the yambo phonon-specific executable <code>yambo_ph</code> and the python utility '''Yambopy'''. The auxiliary code '''LetzElPhC''' (executable <code>lelphc</code>) will be used to obtain the electron-phonon matrix elements by reading the same electronic wavefunctions used by Yambo (and stored in the <code>SAVE</code> directory), while also making full use of crystal symmetries. [https://github.com/yambo-code/LetzElPhC LetzElPhC] will be run by Yambopy, but it must nonetheless be installed. Finally, the exciton-phonon properties can be computed either using <code>yambo_ph</code> or using Yambopy itself. Both cases will be covered in this tutorial.

[[File:Workflow scheme.png|800px|center]]

== Step 0: Pseudopotentials, equilibrium structure and convergence ==

In a real calculation, it is important to ensure that both the pseudopotential and the lattice parameters that we are using are compatible and perform well for the electronic excited states and for the lattice vibrations simultaneously. Furthermore, you have to make sure that the wave function cutoff <code>ecutwfc</code> is converged with respect to the DFPT step and not just to the DFT one. This is in addition to the other customary convergence tests for DFT, DFPT, GW and BSE calculations.

This is often the most time-demanding step when starting on a new system.

For the sake of this tutorial, we assume that we have already done all these tests and we are starting the final workflow to get the exciton-phonon properties.

== Step 1: scf calculation ==

First of all, we run a standard scf calculation with <code>pw.x</code> for Yambo. We stick with non-symmorphic symmetries. At the end, we will have the QE <code>save</code> directory.

This is the input <code>mos2.scf</code>:

&control
wf_collect = .true.,
calculation = "scf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Here we are using full relativistic pseudopotentials from the SG-15 database.

We can run it on our machine (for example using 4 MPI tasks) as:

mpirun -np 4 pw.x -inp mos2.scf > scf.out

== Step 2: nscf calculation for Yambo ==

Copy the QE <code>save</code> directory from the scf calculation and run the nscf calculation for any number of empty states, with the correct <code>k</code>-grid we want to use in Yambo. Here we are using a badly underconverged grid of 6x6x1.

This reciprocal-space grid will also match the momentum transfer <code>q</code> grid on which excitons and phonons will be defined!

The electronic wavefunctions computed at this step and stored in the new nscf <code>save</code> directory will be used both by Yambo and by the electron-phonon code: this is important because using different sets of wavefunctions would lead to a phase mismatch issue in the exciton-phonon matrix elements.

The nscf input <code>mos2.nscf</code> is
&control
wf_collect = .true.,
calculation = "nscf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
nbnd = 250
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Again, we run the calculation

mpirun -np 4 pw.x -inp mos2.nscf > nscf.out

== Step 3: dvscf phonon calculation ==

Now we run the phonon calculation.

Copy the <code>save</code> directory from the '''scf''' calculation and run <code>ph.x</code> for a dvscf calculation with a standard <code>q</code>-grid matching the <code>k</code>-grid we wanna use in Yambo.

At the end, we will have the <code>_ph0</code> directory containing the variation of the self-consistent potential, <math>\Delta V_{SCF}(q)</math>, and the <code>*.dyn</code> files with the phonon energies and eigenvectors.

NB: one could further refine the phonon energies by enforcing the acoustic sum rule, including non-analytic long-range contributions, interpolating to finer grids... all of this can be done within Quantum ESPRESSO and will not be covered in this version of the tutorial.

The input is <code>mos2.dvscf</code>:

mos2_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='mos2',
fildvscf = 'mos2-dvscf',
electron_phonon = 'dvscf',
fildyn='mos2.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=1
/

And now we run as

nohup mpirun -np 8 ph.x -inp mos2.dvscf > dvscf.out &

This time we use nohup and more processes because this calculation may take some time. It is a good idea to set <code>recover=.true.</code> as in a real calculation you will easily breach walltime, and in this way you can safely restart.

== Step 4: create Yambo <code>SAVE</code> directory ==

This is just the standard Yambo initialization: run

p2y

and then

yambo

in the '''nscf''' <code>save</code> folder and then move the newly generated <code>SAVE</code> directory to a convenient place.

== Step 5: run a BSE calculation ==

Now we switch from QE to Yambo. Here, we forgo the GW step for simplicity (we can use a scissor operator to open the band gap).

This calculation has a couple of differences with respect to a standard BSE calculation for optical absorption. We can look at the input file <code>bse.in</code>:

# Runlevels
optics # [R OPT] Optics
rim_cut # [R RIM CUT] Coulomb potential
bss # [R BSS] Bethe Salpeter Equation solver
em1s # [R Xs] Static Inverse Dielectric Matrix
bse # [R BSE] Bethe Salpeter Equation.
bsk # [R BSK] Bethe Salpeter Equation kernel
# RIM and cutoff settings
RandQpts=1000000 # [RIM] Number of random q-points in the BZ
RandGvec= 100 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
# Static screening
X_and_IO_CPU= "1 1 4 2 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
Chimod= "hartree" # [X] IP/Hartree/ALDA/LRC/BSfxc
% BndsRnXs
1 | 200 | # [Xs] Polarization function bands
%
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
NGsBlkXs= 8000 mRy # [Xs] Response block size
# BSE
BS_CPU= "4.1.2" # [PARALLEL] CPUs for each role
BS_ROLEs= "k.eh.t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_diago=4 # [PARALLEL] CPUs for matrix diagonalization
BSEmod= "causal" # [BSE] resonant/causal/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 40000 mRy # [BSK] Exchange components
ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 8000 mRy # [BSK] Screened interaction block size
Lkind="full" #[BSE,X] bar(default)/full/tilde
% KfnQP_E
1.000000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BEnRange
0.00000 | 4.00000 | eV # [BSS] Energy range
%
% BDmRange
0.05000 | 0.05000 | eV # [BSS] Damping range
%
BEnSteps= 2000 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
% BSEQptR
 1 | 7 | # [BSK] Transferred momenta range
%
% BSEBands
25 | 28 | # [BSK] Bands range
%
WRbsWF # [BSS] Write to disk excitonic the FWs

This file was generated using the command: <code> yambo -X s -o b -k sex -y d -r</code> 

First of all, we compute the excitons for all the momenta in the irreducible Brillouin zone for our discrete grid via the <code>BSEQptR</code> variable. This will be a '''finite-momentum''' BSE calculation, analogous to the phonon one.

Second, we change the variable <code>Lkind</code> from <code>bar</code> to <code>full</code>. In Yambo, <code>Lkind="bar"</code>, which is the default for optical absorption, means that we are computing the excitonic response function without the long-range component of the exchange interaction. This cannot be used when computing the exciton momentum dependence, where the long-range exchange interaction can play a role, therefore we have to include it with <code>Lkind="full"</code>. This allows for the calculation of the excitonic longitudinal-transverse splitting (in 3D systems) as well.

We can now run the code:

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

At the end of the calculation, we have obtained the <code>ndb.BS_diago_Q*</code> databases inside the directory <code>bse_Lfull</code>. They contain information on the exciton energies and wavefunctions at each momentum. Do not forget to check the report and logs of your calculation in the same directory to make sure that the code is doing what you want.

== Step 6: obtain the electron-phonon matrix elements ==

We have finished the heavy simulations. Now it's time for the postprocessing. The first order of business is the reconstruction of the electron-phonon coupling matrix elements from the dvscf results and the electronic wavefunctions.

In order to do this, we will run the <code>lelphc</code> executable of the '''LetzElPhC''' code. We will run via command line using yambopy, although it will be instructive to have look at the <code>lelphc</code> input files later.

We run in the same directory where the Yambo <code>SAVE</code> is (remember than you can also virtually move it with a symbolic link).

Type:

yambopy l2y

to see the help for the calculation. Have a look at the various flags and their description.
For example, if we want to do a serial run of LetzElPhC for bands from <math>n_i</math> to <math>n_f</math>, we should type:

yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here <math>n_i</math> and <math>n_f</math> are integers representing the initial and final band indices. '''Please note''': in this case, the band index starts from 1 and the the interval of bands read by the code is <math>[n_i,n_f]</math> including the extrema of the interval.

These should coincide with those used for the Bethe-Salpeter kernel, i.e. those specified in the <code>BSEBands</code> variable of the BSE input file (this is not strictly necessary, but certainly efficient since these calculations use a lot of disk space).

For our system, we want to do a parallel calculation with 4 qpools and 2 kpools. In addition, we want to explicitly specify the path of the <code>lelphc</code> executable and avoid automatically deleting the LetzElPhC data. So we type:

yambopy l2y -ph path/of/dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe

At the end, check your directory: you should find the <code>lelphc.in</code> input file that was run, let's inspect it:

# LetzElPhC input for yambo generated by yambopy
nqpool = 2
nkpool = 4
start_bnd = 25
end_bnd = 28
save_dir = ./SAVE
kernel = dfpt
ph_save_dir = dvscf/ph_save
convention = yambo

Notice the variable <code>convention=yambo</code>: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band <math>n</math> and momentum <math>k-q</math> to band <math>m</math> and momentum <math>k</math>. In the "forward" momentum transfer convention (the more standard one), the transitions go from <math>nk</math> to <math>mk+q</math>. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as <math>\langle mk|dV|nk-q\rangle</math>. This will have consequences also in the formulation of the ''exciton''-phonon coupling matrix element.

You will also find a <code>ndb.elph</code> database that contains the computed electron-phonon matrix elements, in the output format of <code>LetzElPhC</code>.

In addition, if you check the <code>SAVE</code> folder:

ls SAVE/ndb.elph_gkkp*

you will see that yambopy has created the Yambo-compatible electron-phonon databases.

If you want to run LetzElPhC directly, without using yambopy, you can refer to its [https://github.com/yambo-code/LetzElPhC/blob/main/docs/tex_doc/main.pdf|user guide].
Keep in mind that in this case, in order to convert the database to the <code>ndb.elph_gkkp*</code> databases of Yambo, you will then need a couple of lines of python using the Yambopy class <code>ConvertElectronPhononDB</code> in <code>yambopy/letzelph_interface/lelph2y.py</code>.

If you don't need the <code>ndb.elph_gkkp*</code> databases -- e.g., because you are doing the exc-ph postprocessing with yambopy -- you can run <code>yambopy l2y ...</code> with the <code>--no_gkkp</code> option.

If you don't need the <code>ndb.elph</code> databases -- e.g., because you are using <code>yambo_ph</code> for the exc-ph part -- you can run <code>yambopy l2y ...</code> with the <code>--no_lelphc_dbs</code> option.

== Step 7: Obtain the exciton-phonon coupling ==

Now, we can finally access our basic building block for exciton-phonon physics. This could be done entirely in python (using '''Yambopy'''), or by running '''Yambo'''.

* For the '''Yambo postprocessing''' case (running yambo inputs with the <code>yambo_ph</code> executable), '''follow the alternative route to Steps 7-8 [[Exciton-phonon coupling and luminescence - Yambo postprocessing|at this link]]'''.
* For the '''Yambopy postprocessing''' case (using flexible python scripting), '''keep following Steps 7-8 on this page'''.

Our objective is obtaining the following quantity:

<math>\mathcal{G}^\mu_{\alpha\lambda}(0,q)=\sum_{vv^\prime c k} A^{\alpha, *}_{cv^\prime} (k, q) g_{vv^\prime}^\mu (k,q) A^{\lambda}_{cv}(k,0) - \sum_{cc^\prime vk} A^{\alpha, *}_{c^\prime v} (k+q, q) g_{c^\prime c}^\mu (k+q,q) A^{\lambda}_{cv}(k, 0)</math>

Here, <math>A^{\alpha}_{cv}(k,q)</math> are the exciton coefficients extracted from the eigenvector of the two-particles Hamiltonian during the BSE calculation in step 5, while <math>g_{nm}^\mu (k,q)</math> are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton <math>\lambda</math> undergoes phonon-mediated scattering to state <math>\alpha</math> via phonon mode <math>\mu</math>. The scattering can happen for the hole (valence, <math>v</math>) or for the electron (conduction, <math>c</math>).

NB:

(1) This is written in the "backward" momentum transfer convention used by Yambo. The momentum dependence is different in the "forward" transfer convention.

(2) For simplicity, this is written for zero initial exciton momentum. This means that one of the two states involved in the phonon-mediated scattering process will be in the optical limit (and possibly an optically generated exciton), while the other state can have any momentum: this momentum will be the same as the phonon one. This matrix element can be used to describe phonon-assisted absorption and emission spectra.

In order to calculate this quantity using python, we need the <code>ndb.elph</code> databases natively generated by the <code>LetzElPhC</code> code, as well as the BSE databases <code>ndb.BS_diago_Q*</code> containing the information on exciton wavefunctions and energies.

Next, we write a python user script importing the yambopy exciton-phonon tools. You can find a version of this script in <code>yambopy/tutorials/exciton-phonon/calculate_excph.py</code>.

import numpy as np
from yambopy import YamboLatticeDB,YamboWFDB,LetzElphElectronPhononDB #Specific Yambopy dbs loading classes
from yambopy.exciton_phonon.excph_matrix_elements import exciton_phonon_matelem #Specific exc-ph functions

path = '1L_MoS2'
bands_range=[24,28] # 2 valence bands, 2 conduction bands
nexc = 12 # number of excitonic states

bsepath = f'{path}/bse-allq_full' # Path to BSE calculation (Lin=Lout)
savepath = f'{path}/SAVE' # Yambo SAVE
ndb_elph = f'{path}/ndb.elph' # LetzElPhC electron-phonon database (any convention)

# Read lattice
lattice = YamboLatticeDB.from_db_file(filename=f'{savepath}/ns.db1')
# Read electron-phonon
elph = LetzElphElectronPhononDB(ndb_elph,read_all=False)
# Read wave functions
wfcs = YamboWFDB(filename='ns.wf',save=savepath,latdb=lattice,bands_range=bands_range)
# Calculate exciton-phonon matrix elements
exph = exciton_phonon_matelem(lattice,elph,wfcs,BSE_dir=bsepath,nexc_in=nexc,nexc_out=nexc,dmat_mode='save',exph_file='MoS2_Ex-ph.npy')

In this script, we can select the number exciton states <math>\lambda</math> with <code>nexc</code>, the single-particle bands range with <code>bands_range</code>. '''Please note''': in this case, the band index starts from 0 (python indexing) and the the interval of bands read by the code is <math>[n_i,n_f)</math>, therefore the last value (28 in the example) is excluded. That is, in the example we are selecting the 25th, 26th, 27th and 28th bands. Those bands have to be present in both the electron-phonon and BSE calculations.

Here we calculate the couplings of the first twelve states at each finite-<math>q</math> point including <math>q=0</math>. We also include all the nine phonon modes of monolayer MoS2. We also need access to the Yambo <code>SAVE</code> directory in order to read the electronic wavefunctions: they are used to compute the electronic representation matrices (<code>dmat</code> in the code). The argument <code>dmat_mode</code> can be set to <code>'load'</code> for subsequent calculations.

When we are satisfied with the input, we run the code:

python calculate_excph.py

If you check the output, you should find the <code>MoS2_Ex-ph.npy</code> binary file in the directory where you ran.

=== Analysis of the couplings ===

It is a good idea to have a look at what we computed up to now in order to make sure nothing has gone wrong.

It is not easy to know what to expect (apart from symmetry and gauge compliance of the matrix elements), but one can work out the exciton-phonon selection rules in advance, check that the magnitude is reasonable, etc.

It is also not easy to meaningfully plot this quantity. We have to make sure that we are not breaking degenerate states, otherwise the plots will not be invariant.

First of all, we have to know our system: in monolayer MoS2, the first four excitons are all doubly degenerate. The first exciton responsible for a bright peak in the absorption spectrum (the '''A''' peak), is the second state, corresponding to state indices <code>(3,4)</code> in fortran indexing or <code>(2,3)</code> in python indexing.

All these information can be obtained by analyzing the BSE results (this stuff is explained in the BSE tutorials) and by knowledge of the system or class of systems from the literature.

Thus, a good quantity to plot may be the norm of the matrix elements, summed over the degenerate subspace of exciton A, for a certain number of scattered final states mediated by certain phonon modes:

<math>F_A(q)= \sqrt{ \sum_{\alpha \in A,\lambda,\mu} |\mathcal{G}_{\alpha\lambda}^\mu (0,q)|^2 }</math>

In order to do this, we create a python script <code>analyse_excph.py</code> in which we first load the excph dabatases.
You can find a version of this script in the yambopy directory, in <code>tutorials/exciton-phonon</code>.
First, we select the exciton and phonon states to be included in <code>F_A</code>, together with the path of databases and plot details:

# Exciton "in" states
exc_in = [2,3] # First bright peak (A: 2,3 -- B: 6,7)
exc_out = [0,1,2,3] # first 4 states (dispersion of dark triplet state and A)
ph_in = 'all'
# Paths of databases
ns_db1 =f'{path}/SAVE/ns.db1'
ns_ypy = 'MoS2_Ex-ph.npy'

...

Then, we load the data:

# Read lattice and k-space info
ylat = YamboLatticeDB.from_db_file(filename=ns_db1)
print(ylat)

# Load exc-ph database
X_py = np.load(ns_ypy)
G_squared = np.abs(X_py)**2.

The quantity <math>F_A(q)</math> is obtained from a dedicated function as:

if exc_in == 'all': exc_in = range(G_squared.shape[2])
if exc_out == 'all': exc_out = range(G_squared.shape[3])
if ph_in == 'all': ph_in = range(G_squared.shape[1])

G_squared = G_squared[:, ph_in, :, :].sum(axis=(1))
G_squared = G_squared[:, exc_in, :].sum(axis=(1))
G_squared = G_squared[:, exc_out].sum(axis=(1))

F_q = np.sqrt( G_squared )*ha2ev # Switch from Ha to eV

And finally, we have to make a plotting function. For this tutorial we will use a custom scatterplot employing some of the plotting tools provided by yambopy (but you can do whatever you want).

plot_2D_excph(qgrid,G2_to_plot,rlat=ylat.rlat,plt_cbar=True,\
marker='H',s=700,cmap='magma')
...

You can get more experience on using Yambopy for these kinds of visualization by following the [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy Yambopy tutorials]. In fact, remember that these scripts and all the other Yambopy tutorial scripts are just suggestions, not source code written in stone: if you know <code>numpy</code> and <code>matplotlib</code> you can do your own analysis and your own plots, you just need to import the required Yambopy modules to load the data.

In our case, the resulting plot is the following.

[[File:1L MoS2 MoS2 Ex-ph.png|500px|center]]

This can be checked against Fig. 2(d) of reference <ref name="chan2023" />, although you have to keep in mind that our results are badly undersampled in terms of the reciprocal-space grid, as can be easily seen, and the quantity plotted is not exactly the same. However, the main features are already there since they are dictated mostly by crystal symmetries.

Now that we have the exciton-phonon matrix elements, we can use them to build several kinds of observables. Below, we give an example related to phonon-assisted luminescence, but we may update this tutorial in the future to include more cases.

== Step 8: Compute phonon-assisted luminescence ==

[[File:Luminescence scheme.png|250px|right]]

We want to compute the experimental optical signature due to the phonon-assisted recombination of an exciton (as sketched in the figure).

The signal from the phonon replicas can be modeled as a second-order scattering process involving one phonon and one photon:

<math>I^{Sat}(\omega)=\frac{1}{ N_q} \frac{1}{3} \sum_{\epsilon s\beta \mu q} \frac{1}{E_{\beta q}-s\Omega_{\mu q}}\left|\sum_\alpha\frac{ D^{\epsilon}_\alpha \mathcal{G}_{\beta \alpha}^{\mu,*}(q)}{E_\alpha -E_{\beta q} +s\Omega_{\mu q}+\mathrm{i}\eta}\right|^2 \frac{N^{exc}_{\beta q}(T_{exc})[\frac{1+s}{2}+n_{\mu q}(T)]}{\omega -[E_{\beta q}-s\Omega_{\mu q}]+\mathrm{i}\eta}</math>

In this equation, the oscillator strength of the peak is given by the exciton-phonon coupling matrix elements <math>\mathcal{G}</math> multiplied by the exciton dipoles <math>D</math> (they are called "residuals" in Yambo). Here <math>E_\lambda</math> and <math>E_{\alpha q}</math> are the energies of the optical and finite-momentum excitons, respectively, while <math>\Omega_{\mu q}</math> are the phonon energies.

Here, <math>n_{\mu q}(T)</math> is the temperature-dependent phonon Bose-Einstein occupation function. As it can be seen, <math>s=1</math> corresponds to processes of phonon ''emission'' (<math>\propto n(T)+1</math>), while <math>s=-1</math> corresponds to processes of phonon ''absorption'' (<math>\propto n(T)</math>). Therefore, <math>I^{Sat}_{PL}(\omega;T)</math> describes ''light'' emission by recombining excitons mediated by either ''phonon'' absorption or emission.

The quantity <math> N_{\alpha q}(T_{exc})</math> is the exciton occupation function. Luminescence is technically an out-of-equilibrium process, but we can assume that for very low density of excitations and in steady-state conditions, the exciton population can be approximately described by an equilibrium distribution evaluated at an effective temperature. Here, we use the Boltzmann distribution. Experimentally, <math>T_{exc}</math> tends to coincide with the lattice temperature <math>T</math> more or less above 100 K, while at very low temperature (< 10 K), <math>T_{exc}</math> may vary between 10-50 K. It goes without saying that this needs to carefully be checked in your realistic calculations.

Finally, the spectrum is averaged over the polarization directions of the emitted photons (<math>\epsilon=x,y,z</math> representing the respective dipole components).

=== Running the jobs ===

In order to study luminescence in a paradigmatic system, we switch to bulk hexagonal boron nitride and we repeat the workflow. As you can easily see, one can think about automatizing the execution of all these calculations via scripting or more advanced tools. However, in the case of very large simulations (memory-limited or disk-space limited) or for systems whose electronic and lattice properties are fragile with respect to tiny calculation details, one must be very careful and run many basic tests.

Fortunately, we are running a fast underconverged example. We use LDA pseudopotentials from the pseudo-dojo library and the following are the calculations steps.

1. Input <code>hbn.scf</code>:

&control
calculation='scf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR',
outdir = './tmp'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.scf > scf.out

2. Input <code>hbn.nscf</code>:

&control
calculation='nscf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR'
outdir = './'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
nbnd = 120
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.nscf > nscf.out

3. Input <code>hbn.dvscf</code>:

hbn_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='hBN',
fildvscf = 'hBN-dvscf',
electron_phonon = 'dvscf',
fildyn='hBN.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=2
/

nohup mpirun -np 8 pw.x -inp hbn.dvscf > dvscf.out &

4. Input <code>bse.in</code> (we include 2 valence and 2 conduction bands):

optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
em1s
DIP_CPU= "1 8 1" # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_and_IO_CPU= "1 1 1 8 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
BS_CPU= "8 1 1" # [PARALLEL] CPUs for each role
BS_ROLEs= "k eh t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
BS_nCPU_LinAlg_DIAGO=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
% QpntsRXs
1 | 14 | # [Xs] Transferred momenta
%
% BndsRnXs
1 | 120 | # [Xs] Polarization function bands
%
NGsBlkXs= 10 Ry # [Xs] Response block size
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Lkind= "Lfull" # [BSE] Lbar (default) / full
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
% DipBands
1 | 120 | # [DIP] Bands range for dipoles
%
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
DipComputed= "R V P" # [DIP] [default R P V; extra P2 Spin Orb]
BSENGexx= 30000 Ry # [BSK] Exchange components
#ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 9000 Ry # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% KfnQP_E
1.25997 | 1.08816 | 1.12683 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BSEQptR
 1 | 14 | # [BSK] Transferred momenta range
%
% BSEBands
 7 | 10 | # [BSK] Bands range
%
% BEnRange
0.50000 | 8.00000 | eV # [BSS] Energy range
%
% BDmRange
0.050000 | 0.050000 | eV # [BSS] Damping range
%
BEnSteps= 1000 # [BSS] Energy steps
% BLongDir
1.000000 | 1.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
WRbsWF # [BSS] Write to disk excitonic the WFs

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

Importantly, since we want to describe the phonon-assisted recombination process of an *optical* exciton (i.e., emitting a transverse photon), this time we also run an additional calculation at `Q=0` omitting the nonanalytic long-range Coulomb exchange. Make a second input <code>bse_Lbar.in</code> with the following changes:

Lkind= "Lbar" # [BSE] Lbar (default) / full
% BSEQptR
 1 | 1 | # [BSK] Transferred momenta range
%

4b. So now we make a second BSE run in a different directory specified by <code>-J</code>. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases <code>ndb.em1s*</code> that we do not want to recompute from scratch.

mpirun -np 8 yambo -F bse_Lbar.in -J bse_Lbar,bse_Lfull -C bse_Lbar

5. Now we run <code>lelphc</code> (with or without yambopy: in the latter case remember the option <code>-D</code>) to get the el-ph matrix elements, particularly the <code>ndb.elph</code> databases:

yambopy l2y -ph path/of/dvscf/hbn.dvscf -b 7 10 -par 4 2 -D

=== Luminescence calculation ===

6. And finally we calculate exciton-phonon matrix elements and the luminescence spectrum in one go using a python script importing the Yambopy exciton-phonon tools: in order to do this, we need to take a look at all the necessary input variables for the formula written above.

We can start from the script <code>luminescence.py</code> available in <code>tutorials/exciton-phonon</code>.

import numpy as np
import matplotlib.pyplot as plt
from yambopy.exciton_phonon.excph_luminescence import exc_ph_luminescence
from yambopy.exciton_phonon.excph_input_data import exc_ph_get_inputs

We import the necessary tools...

path = '3D_hBN'
# Path to BSE calculation (Lout--> response is Lfull)
bsepath = f'{path}/bse_Lfull' # ndb.BS_diago_Q* databases are needed
# Path to BSE calculation for optically active exciton (Lin --> response is Lbar)
bseBARpath = f'{path}/bse_Lbar' # ndb.BS_diago_Q1 database is needed
# Path to electron-phonon calculation
elphpath = path # ndb.elph is needed
# Path to unprojected dipoles matrix elements (optional)
dipolespath = bsepath # ndb.dipoles is needed (optional)
# Path to lattice and k-space info
savepath = f'{path}/SAVE' # ns.db1 database is needed

... and we specify the paths to the necessary databases. Note that, if we want to perform polarization averaging, we have to recompute the excitonic dipoles in python as seen here.

What about <code>bseBARpath</code>? This variable points to the directory where the databases for the optical (zero-momentum) excitons <math>\alpha</math> (which may be computed with <code>Lkind='Lbar'</code>) is located, which can be different from the directory with the full indirect exciton dispersion <math>\beta</math> (usually computed with <code>Lkind='Lfull'</code>, however <code>Lkind='Ltilde'</code> can also be used if one is interested in "irreducible" excitons). This makes it possible to compute the coupling between different exciton kinds.

bands_range=[6,10] # 2 valence, 2 conduction bands
phonons_range=[0,12] # All phonons
nexc_out = 12 # 12 excitonic states at each momentum (Lout)
nexc_in = 12 # 12 excitonic states at Q=0 (Lin)
T_ph = 10 # Lattice temperature
T_exc = 10 # Effective excitonic temperature

emin=4.4 # Energy range and plot details (in eV)
emax=4.7
estep=0.0002
broad = 0.005 # Broadening parameter for peak width (in eV)

Next, we specify the parameters for the calculation. We include valence bands from 7 to 10 (6 to 9 in python index) and the contribution of all 12 phonon modes. We consider 12 excitonic states for the coupling matrix elements.

Here <code>T_ph</code> and <code>T_exc</code> are the lattice and excitonic temperatures, respectively.

# We calculate and load all the inputs:
# * Exciton-phonon matrix elements
# * Excitonic dipole matrix elements
# * Exciton energies
# * Phonon energies
# * We specify bse_path2=bseBARpath meaning we use Lbar calculation for Q=0 excitons
input_data = exc_ph_get_inputs(savepath,elphpath,bsepath,\
bse_path2=bseBARpath,dipoles_path=dipolespath,\
nexc_in=12,nexc_out=12,\
bands_range=[6,10],phonons_range=[0,12])

ph_energies, exc_energies, exc_energies_in, G, exc_dipoles = input_data

The function <code>exc_ph_get_inputs</code> gives us all the input data for luminescence in the correct format: phonon energies, exciton energies, optical exciton energies (if needed), exciton-phonon matrix elements (calculated on the fly and printed to file), unprojected exciton dipoles (optional, calculated on the fly and printed to file).
The exc-ph matrix element calculation is just a wrapper of the same tools that we have tested in the above section on MoS2.

Finally, there is the calculation of the luminescence spectrum via the function <code>exc_ph_luminescence</code>:
# We calculate the luminescence spectrum including the input data from before
w,PL = exc_ph_luminescence(T_ph,ph_energies,exc_energies,exc_dipoles,G,\
exc_energies_in=exc_energies_in,exc_temp=T_exc,\
nexc_out=nexc_out,nexc_in=nexc_in,emin=emin,emax=emax,\
estep=estep,broad=broad)

If you want to print the luminescence data for later plotting, you can also add the following line to the script:
# Save to file
data = np.column_stack((w, PL))
np.savetxt("hBN_luminescence_12x12x1.dat", data, fmt="%.8f")

It is now time to run the whole script:

python luminescence.py

NB: By using Yambopy for Step 6, we have limited the use of the Yambo code to just step 4 in the entire workflow. This option is more flexible, as it allows for a greater degree of control by the user. On the other hand the Yambo postprocessing route features a Yambo-style input that doesn't require python knowledge and the calculation is currently faster in fortran. However, the luminescence expression computed in Yambo is a slightly different than this one: it is more approximated in the description of the satellite oscillator strengths, but it explicitly includes the renormalization of the direct exciton peak. You can check the differences [[Exciton-phonon coupling and luminescence - Yambo postprocessing|here]].

=== Results ===
We can plot the results of the step 6 calculation. If we do it in the same script we can add something like this:
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_xlim(emin,emax)
ax.set_ylim(0,np.max(PL)*1.1)
ax.get_yaxis().set_visible(False)

ax.plot(w, PL, '-',c='red', label="AA' hBN luminescence")

plt.legend()
plt.savefig('hBN_luminescence.png')
plt.show()

[[File:HBN luminescence satellites.png|500px|center]]

Here, the signal corresponds to a finite-momentum exciton that recombines with the help of several different phonon modes, both optical and acoustic. Each phonon mode whose coupling with the exciton is allowed can generate a peak, and the energy shifts of these peaks with respect to the initial exciton energy correspond to the phonon energies. This result is underconverged, but the main features are all there. In the plot, we show a more converged example using a 12x12x4 grid (all the other parameters being equal). These plots can be compared with Fig. 4(a) of reference <ref name="zanfrognini2023" />.

==Exciton lifetimes==

Using the exciton-phonon coupling it is possible also to calculate exciton life-time due to the scattering with phonon. For this kind of calculation input file can be generated with the command ''yambo -excph l'':

excph # [R] Exction-phonon
ExcGkkp # [R][EXCPH] Exciton-Phonon Matrix Elelements
ExcPhLifeT # [R][EXCPH] Exciton-Phonon Life-Times
BoseTemp= 0.000000 eV # Bosonic Temperature
% ELPhExcStates
1 | 4 | # [EXCPH] Incoming (external) exciton states
%
% ELPhExcSum
1 | 8 | # [EXCPH] Outgoing (virtual) exciton states
%
LoutPath= "none" # [EXCPH] Path of the outgoing L
FANdEtresh= 0.100000E-5 eV # [ELPH] Energy treshold for Fan denominator
% ElPhModes
1 | 12 | # [ELPH] Phonon modes included
%
EXCTemp= 0.000000 eV # [EXCPH] Excitonic Temperature (for luminescence spectra)
AlphaQ= 0.000000 # [EXCPH] Excitonic band structure 2D distortion
#ExcPhOffDiago # [EXCPH] Exciton-Phonon off-diagonal self-energy
#NoMatrxEl # [EXCPH] Set all exciton-phonon matrix elements to one (for testing purpose)
LDamping= 0.1 meV # [EXCPH] Damping of exc-ph self-energy

notice that in this case LDamping is the broadening of the exciton-phonon self energy and it should be a very small value of the order of the phonon life-times around the meV.<ref>[https://arxiv.org/abs/1411.5628 First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong an-harmonicity: The example of palladium hydrides], L. Paulatto, et al. Phys. Rev. B '''91''', 054304 (2015) </ref>. Running this input you will get:

<-> Exciton lifetimes:
<-> Exciton [1] Lifetime: 0.697314 meV
<-> Exciton [2] Lifetime: 0.697314 meV
<-> Exciton [3] Lifetime: 20.22155 meV
<-> Exciton [4] Lifetime: 32.95292 meV

the life-time of the first four excitons. Note that this result is not converged respect to the Q-points. In order to speedup convergence that double-grid<ref name='pierre' /> can be used also in the life-time calculations, see tutorial on electron-phonon coupling for tricks on the convergence of the double-grid.
-->

== References ==

<references>
<ref name="toyozawa2003" >Toyozawa, Yutaka, and Chris Oxlade, ''Optical processes in solids'', [https://m.booksee.org/book/1121964?force_lang=en Cambridge University Press, (2003)]. </ref>
<ref name='lechifflart2023'>P. Lechifflart, F. Paleari, D. Sangalli, C. Attaccalite, ''First-principles study of luminescence in hexagonal boron nitride single layer: Exciton-phonon coupling and the role of substrate'',
[https://doi.org/10.1103/PhysRevMaterials.7.024006 Phys. Rev. M, '''7''' (2), 024006 (2023)]; [https://arxiv.org/abs/2212.10407 arXiv2212.1047]</ref>
<ref name='cannuccia2019'>E. Cannuccia, B. Monserrat and C. Attaccalite, ''Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride'', [https://doi.org/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]; [https://arxiv.org/abs/1807.11797 arXiv1807.11797]</ref>
<ref name='paleari2019'>F. Paleari et al., ''Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride'', [https://doi.org/10.1103/PhysRevLett.122.187401 Phys. Rev. Lett. '''122''', 187401 (2019)]; [https://arxiv.org/abs/1810.08976 arXiv1810.089776] </ref>
<ref name='chen2020'>''Exciton-Phonon Interaction and Relaxation Times from First Principles'',
Hsiao-Yi Chen, Davide Sangalli, and Marco Bernardi, [https://doi.org/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]; [https://arxiv.org/abs/2002.08913 arXiv 2002.08913 (2020)]</ref>
<ref name="lechifflart2023_PhD">P. Lechifflart, ''Exciton-phonon coupling and phonon-assisted luminescence in hexagonal Boron Nitride nanostructures'', [https://hal.science/tel-04266805v1 PhD Thesis, University of Marseille (2023)]; [https://www.yambo-code.eu/wiki/images/5/54/These_final.pdf From the yambo website]</ref>
<ref name='paleari2019_PhD'>F. Paleari, ''First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride'', [https://wwwen.uni.lu/research/fstm/dphyms/people/fulvio_paleari PhD thesis, University of Luxembourg (2019)]</ref>
<ref name='zanfrognini2023'>M. Zanfrognini et al., ''Distinguishing different stackings in layered materials via luminescence spectroscopy'', [https://doi.org/10.1103/PhysRevLett.131.206902 Phys. Rev. Lett. '''131''', 206902 (2023)]; [https://arxiv.org/abs/2305.17554 arXiv 2305.17554] </ref>
<ref name='marini2024'>G. Marini, M. Calandra, P. Cudazzo, ''Optical absorption and photoluminescence of single layer boron nitride from a first principles cumulant approach'', [https://doi.org/10.1021/acs.nanolett.4c00669 Nano Lett., '''24''', 20, 6017 (2024)]; [https://arxiv.org/abs/2402.03826 arXiv 2402.03826 (2024)]</ref>
<ref name='antonius2017'>G. Antonius, S. G. Louie, ''Theory of exciton-phonon coupling'', [https://doi.org/10.1103/PhysRevB.105.085111 Phys. Rev. B, '''105''', 085111 (2022)]; [https://arxiv.org/abs/1705.04245 arXiv1705.04245 (2017)]</ref>
<ref name='paleari2022'> F. Paleari, and A. Marini, ''Exciton-phonon interaction calls for a revision of the “exciton” concept'', [https://doi.org/10.1103/PhysRevB.106.125403 Phys. Rev. B, '''106''', 125403 (2022)]; [https://arxiv.org/abs/2205.02783 arXiv 2205.02783]</ref>
<ref name='cudazzo2020'> P. Cudazzo, ''First-principles description of the exciton-phonon interaction: A cumulant approach'', [https://doi.org/10.1103/PhysRevB.102.045136 Phys. Rev. B, '''102''', 045136 (2020)]; [https://orbilu.uni.lu/bitstream/10993/44769/1/main.pdf Open access pdf from Luxembourg University]</ref>
<ref name='chan2023'>Y-h Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, ''Exciton Lifetime and Optical Line Width Profile via Exciton–Phonon Interactions: Theory and First-Principles Calculations for Monolayer MoS2'', [https://doi.org/10.1021/acs.nanolett.3c00732 Nano Lett., '''23''', 9 (2023)]; [https://arxiv.org/abs/2212.08451 arXiv 2212.08451 (2023)]</ref>
<ref name='murali2025'>M. Nalabothula, S. Reichardt, L. Wirtz, ''Origin of Interlayer Exciton–Phonon Coupling in 2D Heterostructures'', [https://doi.org/10.1021/acs.nanolett.5c00355 Nano Lett., '''25''', 15 (2025)], [https://arxiv.org/abs/2407.16111 arXiv 2407.16111 (2025)]</ref>
</references>

Exciton-phonon coupling and luminescence

2026-04-30T07:57:07Z

Wikiadmin:

[[File:Tdgw-phonon-usc-01-1024x829.jpg|thumb|right|400px|(c) Zhenglu Li, University of Southern California]]

In this advanced tutorial, we will calculate exciton-phonon interactions from first principles by interfacing DFPT (for phonon calculations) and BSE (for exciton calculations).

The DFTP calculations are run with Quantum ESPRESSO, while the many-body GW-BSE calculations are run with Yambo. Finally, the exciton-phonon interaction will be obtained by combining and postprocessing the databases computed in the two previous runs. The great advantage of this workflow is that the calculations can be run in the irreducible Brillouin zones both for the electronic momenta ('''k''') and the transfer momenta ('''Q''', '''q''') of excitons and phonons, thus speeding up considerably the jobs while reducing the IO and memory load.

We will first compute the exciton-phonon coupling matrix elements: these are the building blocks needed to construct experimental observables such as phonon-assisted optical spectra (such as luminescence), Raman spectra and exciton lifetimes. We will do this in the case of monolayer MoS2, a 2D system with large spin-orbit interaction.

As an example of application, we will consider the case of phonon-assisted luminescence. We will do this in the case of bulk hBN, a layered indirect insulator with strong electron-phonon coupling.

'''Note:''' this tutorial will be updated when new exc-ph tools become available in Yambopy (including full-python postprocessing, Raman spectra, interpolated lifetimes, etc).

== Requirements ==

This is an advanced topic: we assume that you already know something about the theory<ref name="toyozawa2003" /><ref name="antonius2017" /><ref name="cudazzo2020" /><ref name="paleari2019_PhD" /><ref name="paleari2022" /><ref name="lechifflart2023_PhD" /> and applications<ref name="paleari2019" /><ref name="cannuccia2019" /><ref name="chen2020" /><ref name="zanfrognini2023" /><ref name="lechifflart2023" /><ref name="chan2023" /><ref name="marini2024" /><ref name="murali2025" /> of exciton-phonon physics.

Also, we assume that you already know how to run both a basic '''Yambo''' GW-BSE calculation and a DFPT phonon calculation with '''Quantum ESPRESSO'''.

Besides the QE executables <code>pw.x</code> and <code>ph.x</code>, we also use the yambo phonon-specific executable <code>yambo_ph</code> and the python utility '''Yambopy'''. The auxiliary code '''LetzElPhC''' (executable <code>lelphc</code>) will be used to obtain the electron-phonon matrix elements by reading the same electronic wavefunctions used by Yambo (and stored in the <code>SAVE</code> directory), while also making full use of crystal symmetries. [https://github.com/yambo-code/LetzElPhC LetzElPhC] will be run by Yambopy, but it must nonetheless be installed. Finally, the exciton-phonon properties can be computed either using <code>yambo_ph</code> or using Yambopy itself. Both cases will be covered in this tutorial.

[[File:Workflow scheme.png|800px|center]]

== Step 0: Pseudopotentials, equilibrium structure and convergence ==

In a real calculation, it is important to ensure that both the pseudopotential and the lattice parameters that we are using are compatible and perform well for the electronic excited states and for the lattice vibrations simultaneously. Furthermore, you have to make sure that the wave function cutoff <code>ecutwfc</code> is converged with respect to the DFPT step and not just to the DFT one. This is in addition to the other customary convergence tests for DFT, DFPT, GW and BSE calculations.

This is often the most time-demanding step when starting on a new system.

For the sake of this tutorial, we assume that we have already done all these tests and we are starting the final workflow to get the exciton-phonon properties.

== Step 1: scf calculation ==

First of all, we run a standard scf calculation with <code>pw.x</code> for Yambo. We stick with non-symmorphic symmetries. At the end, we will have the QE <code>save</code> directory.

This is the input <code>mos2.scf</code>:

&control
wf_collect = .true.,
calculation = "scf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Here we are using full relativistic pseudopotentials from the SG-15 database.

We can run it on our machine (for example using 4 MPI tasks) as:

mpirun -np 4 pw.x -inp mos2.scf > scf.out

== Step 2: nscf calculation for Yambo ==

Copy the QE <code>save</code> directory from the scf calculation and run the nscf calculation for any number of empty states, with the correct <code>k</code>-grid we want to use in Yambo. Here we are using a badly underconverged grid of 6x6x1.

This reciprocal-space grid will also match the momentum transfer <code>q</code> grid on which excitons and phonons will be defined!

The electronic wavefunctions computed at this step and stored in the new nscf <code>save</code> directory will be used both by Yambo and by the electron-phonon code: this is important because using different sets of wavefunctions would lead to a phase mismatch issue in the exciton-phonon matrix elements.

The nscf input <code>mos2.nscf</code> is
&control
wf_collect = .true.,
calculation = "nscf",
verbosity = 'high',
pseudo_dir = '$PSEUDO_DIR',
prefix = "mos2",
outdir = '.',
/&end
&system
ecutwfc = 100.0,
occupations = 'fixed',
ibrav = 4,
celldm(1) = 5.9000811881,
celldm(3) = 6.7795677253,
nat = 3,
ntyp = 2,
lspinorb = .true.
noncolin = .true.
nbnd = 250
assume_isolated = '2D'
force_symmorphic = .true.
/&end
&electrons
electron_maxstep = 200,
mixing_beta = 0.7,
conv_thr = 1.d-08,
/&end
ATOMIC_SPECIES
Mo 95.940 Mo_ONCV_PBE_FR-1.0.upf
S 32.065 S_ONCV_PBE_FR-1.1.upf
ATOMIC_POSITIONS { crystal }
Mo 0.333333333 0.666666667 0.000000000
S 0.666666667 0.333333333 0.073413577
S 0.666666667 0.333333333 -0.073413577
K_POINTS { automatic }
6 6 1 0 0 0

Again, we run the calculation

mpirun -np 4 pw.x -inp mos2.nscf > nscf.out

== Step 3: dvscf phonon calculation ==

Now we run the phonon calculation.

Copy the <code>save</code> directory from the '''scf''' calculation and run <code>ph.x</code> for a dvscf calculation with a standard <code>q</code>-grid matching the <code>k</code>-grid we wanna use in Yambo.

At the end, we will have the <code>_ph0</code> directory containing the variation of the self-consistent potential, <math>\Delta V_{SCF}(q)</math>, and the <code>*.dyn</code> files with the phonon energies and eigenvectors.

NB: one could further refine the phonon energies by enforcing the acoustic sum rule, including non-analytic long-range contributions, interpolating to finer grids... all of this can be done within Quantum ESPRESSO and will not be covered in this version of the tutorial.

The input is <code>mos2.dvscf</code>:

mos2_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='mos2',
fildvscf = 'mos2-dvscf',
electron_phonon = 'dvscf',
fildyn='mos2.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=1
/

And now we run as

nohup mpirun -np 8 ph.x -inp mos2.dvscf > dvscf.out &

This time we use nohup and more processes because this calculation may take some time. It is a good idea to set <code>recover=.true.</code> as in a real calculation you will easily breach walltime, and in this way you can safely restart.

== Step 4: create Yambo <code>SAVE</code> directory ==

This is just the standard Yambo initialization: run

p2y

and then

yambo

in the '''nscf''' <code>save</code> folder and then move the newly generated <code>SAVE</code> directory to a convenient place.

== Step 5: run a BSE calculation ==

Now we switch from QE to Yambo. Here, we forgo the GW step for simplicity (we can use a scissor operator to open the band gap).

This calculation has a couple of differences with respect to a standard BSE calculation for optical absorption. We can look at the input file <code>bse.in</code>:

# Runlevels
optics # [R OPT] Optics
rim_cut # [R RIM CUT] Coulomb potential
bss # [R BSS] Bethe Salpeter Equation solver
em1s # [R Xs] Static Inverse Dielectric Matrix
bse # [R BSE] Bethe Salpeter Equation.
bsk # [R BSK] Bethe Salpeter Equation kernel
# RIM and cutoff settings
RandQpts=1000000 # [RIM] Number of random q-points in the BZ
RandGvec= 100 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
# Static screening
X_and_IO_CPU= "1 1 4 2 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
Chimod= "hartree" # [X] IP/Hartree/ALDA/LRC/BSfxc
% BndsRnXs
1 | 200 | # [Xs] Polarization function bands
%
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
NGsBlkXs= 8000 mRy # [Xs] Response block size
# BSE
BS_CPU= "4.1.2" # [PARALLEL] CPUs for each role
BS_ROLEs= "k.eh.t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_diago=4 # [PARALLEL] CPUs for matrix diagonalization
BSEmod= "causal" # [BSE] resonant/causal/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 40000 mRy # [BSK] Exchange components
ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 8000 mRy # [BSK] Screened interaction block size
Lkind="full" #[BSE,X] bar(default)/full/tilde
% KfnQP_E
1.000000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BEnRange
0.00000 | 4.00000 | eV # [BSS] Energy range
%
% BDmRange
0.05000 | 0.05000 | eV # [BSS] Damping range
%
BEnSteps= 2000 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
% BSEQptR
 1 | 7 | # [BSK] Transferred momenta range
%
% BSEBands
25 | 28 | # [BSK] Bands range
%
WRbsWF # [BSS] Write to disk excitonic the FWs

This file was generated using the command: <code> yambo -X s -o b -k sex -y d -r</code> 

First of all, we compute the excitons for all the momenta in the irreducible Brillouin zone for our discrete grid via the <code>BSEQptR</code> variable. This will be a '''finite-momentum''' BSE calculation, analogous to the phonon one.

Second, we change the variable <code>Lkind</code> from <code>bar</code> to <code>full</code>. In Yambo, <code>Lkind="bar"</code>, which is the default for optical absorption, means that we are computing the excitonic response function without the long-range component of the exchange interaction. This cannot be used when computing the exciton momentum dependence, where the long-range exchange interaction can play a role, therefore we have to include it with <code>Lkind="full"</code>. This allows for the calculation of the excitonic longitudinal-transverse splitting (in 3D systems) as well.

We can now run the code:

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

At the end of the calculation, we have obtained the <code>ndb.BS_diago_Q*</code> databases inside the directory <code>bse_Lfull</code>. They contain information on the exciton energies and wavefunctions at each momentum. Do not forget to check the report and logs of your calculation in the same directory to make sure that the code is doing what you want.

== Step 6: obtain the electron-phonon matrix elements ==

We have finished the heavy simulations. Now it's time for the postprocessing. The first order of business is the reconstruction of the electron-phonon coupling matrix elements from the dvscf results and the electronic wavefunctions.

In order to do this, we will run the <code>lelphc</code> executable of the '''LetzElPhC''' code. We will run via command line using yambopy, although it will be instructive to have look at the <code>lelphc</code> input files later.

We run in the same directory where the Yambo <code>SAVE</code> is (remember than you can also virtually move it with a symbolic link).

Type:

yambopy l2y

to see the help for the calculation. Have a look at the various flags and their description.
For example, if we want to do a serial run of LetzElPhC for bands from <math>n_i</math> to <math>n_f</math>, we should type:

yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here <math>n_i</math> and <math>n_f</math> are integers representing the initial and final band indices. '''Please note''': in this case, the band index starts from 1 and the the interval of bands read by the code is <math>[n_i,n_f]</math> including the extrema of the interval.

These should coincide with those used for the Bethe-Salpeter kernel, i.e. those specified in the <code>BSEBands</code> variable of the BSE input file (this is not strictly necessary, but certainly efficient since these calculations use a lot of disk space).

For our system, we want to do a parallel calculation with 4 qpools and 2 kpools. In addition, we want to explicitly specify the path of the <code>lelphc</code> executable and avoid automatically deleting the LetzElPhC data. So we type:

yambopy l2y -ph path/of/dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe

At the end, check your directory: you should find the <code>lelphc.in</code> input file that was run, let's inspect it:

# LetzElPhC input for yambo generated by yambopy
nqpool = 2
nkpool = 4
start_bnd = 25
end_bnd = 28
save_dir = ./SAVE
kernel = dfpt
ph_save_dir = dvscf/ph_save
convention = yambo

Notice the variable <code>convention=yambo</code>: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band <math>n</math> and momentum <math>k-q</math> to band <math>m</math> and momentum <math>k</math>. In the "forward" momentum transfer convention (the more standard one), the transitions go from <math>nk</math> to <math>mk+q</math>. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as <math>\langle mk|dV|nk-q\rangle</math>. This will have consequences also in the formulation of the ''exciton''-phonon coupling matrix element.

You will also find a <code>ndb.elph</code> database that contains the computed electron-phonon matrix elements, in the output format of <code>LetzElPhC</code>.

In addition, if you check the <code>SAVE</code> folder:

ls SAVE/ndb.elph_gkkp*

you will see that yambopy has created the Yambo-compatible electron-phonon databases.

If you want to run LetzElPhC directly, without using yambopy, you can refer to its [https://github.com/yambo-code/LetzElPhC/blob/main/docs/tex_doc/main.pdf|user guide].
Keep in mind that in this case, in order to convert the database to the <code>ndb.elph_gkkp*</code> databases of Yambo, you will then need a couple of lines of python using the Yambopy class <code>ConvertElectronPhononDB</code> in <code>yambopy/letzelph_interface/lelph2y.py</code>.

If you don't need the <code>ndb.elph_gkkp*</code> databases -- e.g., because you are doing the exc-ph postprocessing with yambopy -- you can run <code>yambopy l2y ...</code> with the <code>--no_gkkp</code> option.

If you don't need the <code>ndb.elph</code> databases -- e.g., because you are using <code>yambo_ph</code> for the exc-ph part -- you can run <code>yambopy l2y ...</code> with the <code>--no_lelphc_dbs</code> option.

== Step 7: Obtain the exciton-phonon coupling ==

Now, we can finally access our basic building block for exciton-phonon physics. This could be done entirely in python (using '''Yambopy'''), or by running '''Yambo'''.

* For the '''Yambo postprocessing''' case (running yambo inputs with the <code>yambo_ph</code> executable), '''follow the alternative route to Steps 7-8 [[Exciton-phonon coupling and luminescence - Yambo postprocessing|at this link]]'''.
* For the '''Yambopy postprocessing''' case (using flexible python scripting), '''keep following Steps 7-8 on this page'''.

Our objective is obtaining the following quantity:

<math>\mathcal{G}^\mu_{\alpha\lambda}(0,q)=\sum_{vv^\prime c k} A^{\alpha, *}_{cv^\prime} (k, q) g_{vv^\prime}^\mu (k,q) A^{\lambda}_{cv}(k,0) - \sum_{cc^\prime vk} A^{\alpha, *}_{c^\prime v} (k+q, q) g_{c^\prime c}^\mu (k+q,q) A^{\lambda}_{cv}(k, 0)</math>

Here, <math>A^{\alpha}_{cv}(k,q)</math> are the exciton coefficients extracted from the eigenvector of the two-particles Hamiltonian during the BSE calculation in step 5, while <math>g_{nm}^\mu (k,q)</math> are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton <math>\lambda</math> undergoes phonon-mediated scattering to state <math>\alpha</math> via phonon mode <math>\mu</math>. The scattering can happen for the hole (valence, <math>v</math>) or for the electron (conduction, <math>c</math>).

NB:

(1) This is written in the "backward" momentum transfer convention used by Yambo. The momentum dependence is different in the "forward" transfer convention.

(2) For simplicity, this is written for zero initial exciton momentum. This means that one of the two states involved in the phonon-mediated scattering process will be in the optical limit (and possibly an optically generated exciton), while the other state can have any momentum: this momentum will be the same as the phonon one. This matrix element can be used to describe phonon-assisted absorption and emission spectra.

In order to calculate this quantity using python, we need the <code>ndb.elph</code> databases natively generated by the <code>LetzElPhC</code> code, as well as the BSE databases <code>ndb.BS_diago_Q*</code> containing the information on exciton wavefunctions and energies.

Next, we write a python user script importing the yambopy exciton-phonon tools. You can find a version of this script in <code>yambopy/tutorials/exciton-phonon/calculate_excph.py</code>.

import numpy as np
from yambopy import YamboLatticeDB,YamboWFDB,LetzElphElectronPhononDB #Specific Yambopy dbs loading classes
from yambopy.exciton_phonon.excph_matrix_elements import exciton_phonon_matelem #Specific exc-ph functions

path = '1L_MoS2'
bands_range=[24,28] # 2 valence bands, 2 conduction bands
nexc = 12 # number of excitonic states

bsepath = f'{path}/bse-allq_full' # Path to BSE calculation (Lin=Lout)
savepath = f'{path}/SAVE' # Yambo SAVE
ndb_elph = f'{path}/ndb.elph' # LetzElPhC electron-phonon database (any convention)

# Read lattice
lattice = YamboLatticeDB.from_db_file(filename=f'{savepath}/ns.db1')
# Read electron-phonon
elph = LetzElphElectronPhononDB(ndb_elph,read_all=False)
# Read wave functions
wfcs = YamboWFDB(filename='ns.wf',save=savepath,latdb=lattice,bands_range=bands_range)
# Calculate exciton-phonon matrix elements
exph = exciton_phonon_matelem(lattice,elph,wfcs,BSE_dir=bsepath,nexc_in=nexc,nexc_out=nexc,dmat_mode='save',exph_file='MoS2_Ex-ph.npy')

In this script, we can select the number exciton states <math>\lambda</math> with <code>nexc</code>, the single-particle bands range with <code>bands_range</code>. '''Please note''': in this case, the band index starts from 0 (python indexing) and the the interval of bands read by the code is <math>[n_i,n_f)</math>, therefore the last value (28 in the example) is excluded. That is, in the example we are selecting the 25th, 26th, 27th and 28th bands. Those bands have to be present in both the electron-phonon and BSE calculations.

Here we calculate the couplings of the first twelve states at each finite-<math>q</math> point including <math>q=0</math>. We also include all the nine phonon modes of monolayer MoS2. We also need access to the Yambo <code>SAVE</code> directory in order to read the electronic wavefunctions: they are used to compute the electronic representation matrices (<code>dmat</code> in the code). The argument <code>dmat_mode</code> can be set to <code>'load'</code> for subsequent calculations.

When we are satisfied with the input, we run the code:

python calculate_excph.py

If you check the output, you should find the <code>MoS2_Ex-ph.npy</code> binary file in the directory where you ran.

=== Analysis of the couplings ===

It is a good idea to have a look at what we computed up to now in order to make sure nothing has gone wrong.

It is not easy to know what to expect (apart from symmetry and gauge compliance of the matrix elements), but one can work out the exciton-phonon selection rules in advance, check that the magnitude is reasonable, etc.

It is also not easy to meaningfully plot this quantity. We have to make sure that we are not breaking degenerate states, otherwise the plots will not be invariant.

First of all, we have to know our system: in monolayer MoS2, the first four excitons are all doubly degenerate. The first exciton responsible for a bright peak in the absorption spectrum (the '''A''' peak), is the second state, corresponding to state indices <code>(3,4)</code> in fortran indexing or <code>(2,3)</code> in python indexing.

All these information can be obtained by analyzing the BSE results (this stuff is explained in the BSE tutorials) and by knowledge of the system or class of systems from the literature.

Thus, a good quantity to plot may be the norm of the matrix elements, summed over the degenerate subspace of exciton A, for a certain number of scattered final states mediated by certain phonon modes:

<math>F_A(q)= \sqrt{ \sum_{\alpha \in A,\lambda,\mu} |\mathcal{G}_{\alpha\lambda}^\mu (0,q)|^2 }</math>

In order to do this, we create a python script <code>analyse_excph.py</code> in which we first load the excph dabatases.
You can find a version of this script in the yambopy directory, in <code>tutorials/exciton-phonon</code>.
First, we select the exciton and phonon states to be included in <code>F_A</code>, together with the path of databases and plot details:

# Exciton "in" states
exc_in = [2,3] # First bright peak (A: 2,3 -- B: 6,7)
exc_out = [0,1,2,3] # first 4 states (dispersion of dark triplet state and A)
ph_in = 'all'
# Paths of databases
ns_db1 =f'{path}/SAVE/ns.db1'
ns_ypy = 'MoS2_Ex-ph.npy'

...

Then, we load the data:

# Read lattice and k-space info
ylat = YamboLatticeDB.from_db_file(filename=ns_db1)
print(ylat)

# Load exc-ph database
X_py = np.load(ns_ypy)
G_squared = np.abs(X_py)**2.

The quantity <math>F_A(q)</math> is obtained from a dedicated function as:

if exc_in == 'all': exc_in = range(G_squared.shape[2])
if exc_out == 'all': exc_out = range(G_squared.shape[3])
if ph_in == 'all': ph_in = range(G_squared.shape[1])

G_squared = G_squared[:, ph_in, :, :].sum(axis=(1))
G_squared = G_squared[:, exc_in, :].sum(axis=(1))
G_squared = G_squared[:, exc_out].sum(axis=(1))

F_q = np.sqrt( G_squared )*ha2ev # Switch from Ha to eV

And finally, we have to make a plotting function. For this tutorial we will use a custom scatterplot employing some of the plotting tools provided by yambopy (but you can do whatever you want).

plot_2D_excph(qgrid,G2_to_plot,rlat=ylat.rlat,plt_cbar=True,\
marker='H',s=700,cmap='magma')
...

You can get more experience on using Yambopy for these kinds of visualization by following the [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy Yambopy tutorials]. In fact, remember that these scripts and all the other Yambopy tutorial scripts are just suggestions, not source code written in stone: if you know <code>numpy</code> and <code>matplotlib</code> you can do your own analysis and your own plots, you just need to import the required Yambopy modules to load the data.

In our case, the resulting plot is the following.

[[File:1L MoS2 MoS2 Ex-ph.png|500px|center]]

This can be checked against Fig. 2(d) of reference <ref name="chan2023" />, although you have to keep in mind that our results are badly undersampled in terms of the reciprocal-space grid, as can be easily seen, and the quantity plotted is not exactly the same. However, the main features are already there since they are dictated mostly by crystal symmetries.

Now that we have the exciton-phonon matrix elements, we can use them to build several kinds of observables. Below, we give an example related to phonon-assisted luminescence, but we may update this tutorial in the future to include more cases.

== Step 8: Compute phonon-assisted luminescence ==

[[File:Luminescence scheme.png|250px|right]]

We want to compute the experimental optical signature due to the phonon-assisted recombination of an exciton (as sketched in the figure).

The signal from the phonon replicas can be modeled as a second-order scattering process involving one phonon and one photon:

<math>I^{Sat}(\omega)=\frac{1}{ N_q} \frac{1}{3} \sum_{\epsilon s\beta \mu q} \frac{1}{E_{\beta q}-s\Omega_{\mu q}}\left|\sum_\alpha\frac{ D^{\epsilon}_\alpha \mathcal{G}_{\beta \alpha}^{\mu,*}(q)}{E_\alpha -E_{\beta q} +s\Omega_{\mu q}+\mathrm{i}\eta}\right|^2 \frac{N^{exc}_{\beta q}(T_{exc})[\frac{1+s}{2}+n_{\mu q}(T)]}{\omega -[E_{\beta q}-s\Omega_{\mu q}]+\mathrm{i}\eta}</math>

In this equation, the oscillator strength of the peak is given by the exciton-phonon coupling matrix elements <math>\mathcal{G}</math> multiplied by the exciton dipoles <math>D</math> (they are called "residuals" in Yambo). Here <math>E_\lambda</math> and <math>E_{\alpha q}</math> are the energies of the optical and finite-momentum excitons, respectively, while <math>\Omega_{\mu q}</math> are the phonon energies.

Here, <math>n_{\mu q}(T)</math> is the temperature-dependent phonon Bose-Einstein occupation function. As it can be seen, <math>s=1</math> corresponds to processes of phonon ''emission'' (<math>\propto n(T)+1</math>), while <math>s=-1</math> corresponds to processes of phonon ''absorption'' (<math>\propto n(T)</math>). Therefore, <math>I^{Sat}_{PL}(\omega;T)</math> describes ''light'' emission by recombining excitons mediated by either ''phonon'' absorption or emission.

The quantity <math> N_{\alpha q}(T_{exc})</math> is the exciton occupation function. Luminescence is technically an out-of-equilibrium process, but we can assume that for very low density of excitations and in steady-state conditions, the exciton population can be approximately described by an equilibrium distribution evaluated at an effective temperature. Here, we use the Boltzmann distribution. Experimentally, <math>T_{exc}</math> tends to coincide with the lattice temperature <math>T</math> more or less above 100 K, while at very low temperature (< 10 K), <math>T_{exc}</math> may vary between 10-50 K. It goes without saying that this needs to carefully be checked in your realistic calculations.

Finally, the spectrum is averaged over the polarization directions of the emitted photons (<math>\epsilon=x,y,z</math> representing the respective dipole components).

=== Running the jobs ===

In order to study luminescence in a paradigmatic system, we switch to bulk hexagonal boron nitride and we repeat the workflow. As you can easily see, one can think about automatizing the execution of all these calculations via scripting or more advanced tools. However, in the case of very large simulations (memory-limited or disk-space limited) or for systems whose electronic and lattice properties are fragile with respect to tiny calculation details, one must be very careful and run many basic tests.

Fortunately, we are running a fast underconverged example. We use LDA pseudopotentials from the pseudo-dojo library and the following are the calculations steps.

1. Input <code>hbn.scf</code>:

&control
calculation='scf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR',
outdir = './tmp'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.scf > scf.out

2. Input <code>hbn.nscf</code>:

&control
calculation='nscf',
prefix='hBN',
restart_mode='from_scratch'
pseudo_dir = '$PSEUDO_DIR'
outdir = './'
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 4,
celldm(1) = 4.703675849
celldm(3) = 2.603711434
nat= 4,
ntyp= 2,
force_symmorphic=.true.
ecutwfc = 100,
nbnd = 120
/
&electrons
diago_david_ndim = 2
diago_full_acc=.true.
diago_thr_init=5.0e-6
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
B 10.81100 B_LDA_dojo.UPF
N 14.00674 N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N 0.6666666670 0.3333333330 -0.250000000000
B 0.3333333330 0.6666666670 -0.250000000000
B 0.6666666670 0.3333333330 0.25000000000
N 0.3333333330 0.6666666670 0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.nscf > nscf.out

3. Input <code>hbn.dvscf</code>:

hbn_dvscf
&inputph
tr2_ph=1.0d-12,
verbosity='high'
prefix='hBN',
fildvscf = 'hBN-dvscf',
electron_phonon = 'dvscf',
fildyn='hBN.dyn',
epsil=.false.,
ldisp=.true.,
recover=.true.,
nq1=6,
nq2=6,
nq3=2
/

nohup mpirun -np 8 pw.x -inp hbn.dvscf > dvscf.out &

4. Input <code>bse.in</code> (we include 2 valence and 2 conduction bands):

optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
em1s
DIP_CPU= "1 8 1" # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_and_IO_CPU= "1 1 1 8 1" # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
BS_CPU= "8 1 1" # [PARALLEL] CPUs for each role
BS_ROLEs= "k eh t" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
BS_nCPU_LinAlg_DIAGO=-1 # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
% QpntsRXs
1 | 14 | # [Xs] Transferred momenta
%
% BndsRnXs
1 | 120 | # [Xs] Polarization function bands
%
NGsBlkXs= 10 Ry # [Xs] Response block size
% LongDrXs
1.000000 | 1.000000 | 1.000000 | # [Xs] [cc] Electric Field
%
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Lkind= "Lfull" # [BSE] Lbar (default) / full
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
% DipBands
1 | 120 | # [DIP] Bands range for dipoles
%
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
DipComputed= "R V P" # [DIP] [default R P V; extra P2 Spin Orb]
BSENGexx= 30000 Ry # [BSK] Exchange components
#ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 9000 Ry # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% KfnQP_E
1.25997 | 1.08816 | 1.12683 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BSEQptR
 1 | 14 | # [BSK] Transferred momenta range
%
% BSEBands
 7 | 10 | # [BSK] Bands range
%
% BEnRange
0.50000 | 8.00000 | eV # [BSS] Energy range
%
% BDmRange
0.050000 | 0.050000 | eV # [BSS] Damping range
%
BEnSteps= 1000 # [BSS] Energy steps
% BLongDir
1.000000 | 1.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
WRbsWF # [BSS] Write to disk excitonic the WFs

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

Importantly, since we want to describe the phonon-assisted recombination process of an *optical* exciton (i.e., emitting a transverse photon), this time we also run an additional calculation at `Q=0` omitting the nonanalytic long-range Coulomb exchange. Make a second input <code>bse_Lbar.in</code> with the following changes:

Lkind= "Lbar" # [BSE] Lbar (default) / full
% BSEQptR
 1 | 1 | # [BSK] Transferred momenta range
%

4b. So now we make a second BSE run in a different directory specified by <code>-J</code>. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases <code>ndb.em1s*</code> that we do not want to recompute from scratch.

mpirun -np 8 yambo -F bse_Lbar.in -J bse_Lbar,bse_Lfull -C bse_Lbar

5. Now we run <code>lelphc</code> (with or without yambopy: in the latter case remember the option <code>-D</code>) to get the el-ph matrix elements, particularly the <code>ndb.elph</code> databases:

yambopy l2y -ph path/of/dvscf/hbn.dvscf -b 7 10 -par 4 2 -D

=== Luminescence calculation ===

6. And finally we calculate exciton-phonon matrix elements and the luminescence spectrum in one go using a python script importing the Yambopy exciton-phonon tools: in order to do this, we need to take a look at all the necessary input variables for the formula written above.

We can start from the script <code>luminescence.py</code> available in <code>tutorials/exciton-phonon</code>.

import numpy as np
import matplotlib.pyplot as plt
from yambopy.exciton_phonon.excph_luminescence import exc_ph_luminescence
from yambopy.exciton_phonon.excph_input_data import exc_ph_get_inputs

We import the necessary tools...

path = '3D_hBN'
# Path to BSE calculation (Lout--> response is Lfull)
bsepath = f'{path}/bse_Lfull' # ndb.BS_diago_Q* databases are needed
# Path to BSE calculation for optically active exciton (Lin --> response is Lbar)
bseBARpath = f'{path}/bse_Lbar' # ndb.BS_diago_Q1 database is needed
# Path to electron-phonon calculation
elphpath = path # ndb.elph is needed
# Path to unprojected dipoles matrix elements (optional)
dipolespath = bsepath # ndb.dipoles is needed (optional)
# Path to lattice and k-space info
savepath = f'{path}/SAVE' # ns.db1 database is needed

... and we specify the paths to the necessary databases. Note that, if we want to perform polarization averaging, we have to recompute the excitonic dipoles in python as seen here.

What about <code>bseBARpath</code>? This variable points to the directory where the databases for the optical (zero-momentum) excitons <math>\alpha</math> (which may be computed with <code>Lkind='Lbar'</code>) is located, which can be different from the directory with the full indirect exciton dispersion <math>\beta</math> (usually computed with <code>Lkind='Lfull'</code>, however <code>Lkind='Ltilde'</code> can also be used if one is interested in "irreducible" excitons). This makes it possible to compute the coupling between different exciton kinds.

bands_range=[6,10] # 2 valence, 2 conduction bands
phonons_range=[0,12] # All phonons
nexc_out = 12 # 12 excitonic states at each momentum (Lout)
nexc_in = 12 # 12 excitonic states at Q=0 (Lin)
T_ph = 10 # Lattice temperature
T_exc = 10 # Effective excitonic temperature

emin=4.4 # Energy range and plot details (in eV)
emax=4.7
estep=0.0002
broad = 0.005 # Broadening parameter for peak width (in eV)

Next, we specify the parameters for the calculation. We include valence bands from 7 to 10 (6 to 9 in python index) and the contribution of all 12 phonon modes. We consider 12 excitonic states for the coupling matrix elements.

Here <code>T_ph</code> and <code>T_exc</code> are the lattice and excitonic temperatures, respectively.

# We calculate and load all the inputs:
# * Exciton-phonon matrix elements
# * Excitonic dipole matrix elements
# * Exciton energies
# * Phonon energies
# * We specify bse_path2=bseBARpath meaning we use Lbar calculation for Q=0 excitons
input_data = exc_ph_get_inputs(savepath,elphpath,bsepath,\
bse_path2=bseBARpath,dipoles_path=dipolespath,\
nexc_in=12,nexc_out=12,\
bands_range=[6,10],phonons_range=[0,12])

ph_energies, exc_energies, exc_energies_in, G, exc_dipoles = input_data

The function <code>exc_ph_get_inputs</code> gives us all the input data for luminescence in the correct format: phonon energies, exciton energies, optical exciton energies (if needed), exciton-phonon matrix elements (calculated on the fly and printed to file), unprojected exciton dipoles (optional, calculated on the fly and printed to file).
The exc-ph matrix element calculation is just a wrapper of the same tools that we have tested in the above section on MoS2.

Finally, there is the calculation of the luminescence spectrum via the function <code>exc_ph_luminescence</code>:
# We calculate the luminescence spectrum including the input data from before
w,PL = exc_ph_luminescence(T_ph,ph_energies,exc_energies,exc_dipoles,G,\
exc_energies_in=exc_energies_in,exc_temp=T_exc,\
nexc_out=nexc_out,nexc_in=nexc_in,emin=emin,emax=emax,\
estep=estep,broad=broad)

If you want to print the luminescence data for later plotting, you can also add the following line to the script:
# Save to file
data = np.column_stack((w, PL))
np.savetxt("hBN_luminescence_12x12x1.dat", data, fmt="%.8f")

It is now time to run the whole script:

python luminescence.py

NB: By using Yambopy for Step 6, we have limited the use of the Yambo code to just step 4 in the entire workflow. This option is more flexible, as it allows for a greater degree of control by the user. On the other hand the Yambo postprocessing route features a Yambo-style input that doesn't require python knowledge and the calculation is currently faster in fortran. However, the luminescence expression computed in Yambo is a slightly different than this one: it is more approximated in the description of the satellite oscillator strengths, but it explicitly includes the renormalization of the direct exciton peak. You can check the differences [[Exciton-phonon coupling and luminescence - Yambo postprocessing|here]].

=== Results ===
We can plot the results of the step 6 calculation. If we do it in the same script we can add something like this:
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_xlim(emin,emax)
ax.set_ylim(0,np.max(PL)*1.1)
ax.get_yaxis().set_visible(False)

ax.plot(w, PL, '-',c='red', label="AA' hBN luminescence")

plt.legend()
plt.savefig('hBN_luminescence.png')
plt.show()

[[File:HBN luminescence satellites.png|500px|center]]

Here, the signal corresponds to a finite-momentum exciton that recombines with the help of several different phonon modes, both optical and acoustic. Each phonon mode whose coupling with the exciton is allowed can generate a peak, and the energy shifts of these peaks with respect to the initial exciton energy correspond to the phonon energies. This result is underconverged, but the main features are all there. In the plot, we show a more converged example using a 12x12x4 grid (all the other parameters being equal). These plots can be compared with Fig. 4(a) of reference <ref name="zanfrognini2023" />.

== References ==

<references>
<ref name="toyozawa2003" >Toyozawa, Yutaka, and Chris Oxlade, ''Optical processes in solids'', [https://m.booksee.org/book/1121964?force_lang=en Cambridge University Press, (2003)]. </ref>
<ref name='lechifflart2023'>P. Lechifflart, F. Paleari, D. Sangalli, C. Attaccalite, ''First-principles study of luminescence in hexagonal boron nitride single layer: Exciton-phonon coupling and the role of substrate'',
[https://doi.org/10.1103/PhysRevMaterials.7.024006 Phys. Rev. M, '''7''' (2), 024006 (2023)]; [https://arxiv.org/abs/2212.10407 arXiv2212.1047]</ref>
<ref name='cannuccia2019'>E. Cannuccia, B. Monserrat and C. Attaccalite, ''Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride'', [https://doi.org/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]; [https://arxiv.org/abs/1807.11797 arXiv1807.11797]</ref>
<ref name='paleari2019'>F. Paleari et al., ''Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride'', [https://doi.org/10.1103/PhysRevLett.122.187401 Phys. Rev. Lett. '''122''', 187401 (2019)]; [https://arxiv.org/abs/1810.08976 arXiv1810.089776] </ref>
<ref name='chen2020'>''Exciton-Phonon Interaction and Relaxation Times from First Principles'',
Hsiao-Yi Chen, Davide Sangalli, and Marco Bernardi, [https://doi.org/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]; [https://arxiv.org/abs/2002.08913 arXiv 2002.08913 (2020)]</ref>
<ref name="lechifflart2023_PhD">P. Lechifflart, ''Exciton-phonon coupling and phonon-assisted luminescence in hexagonal Boron Nitride nanostructures'', [https://hal.science/tel-04266805v1 PhD Thesis, University of Marseille (2023)]; [https://www.yambo-code.eu/wiki/images/5/54/These_final.pdf From the yambo website]</ref>
<ref name='paleari2019_PhD'>F. Paleari, ''First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride'', [https://wwwen.uni.lu/research/fstm/dphyms/people/fulvio_paleari PhD thesis, University of Luxembourg (2019)]</ref>
<ref name='zanfrognini2023'>M. Zanfrognini et al., ''Distinguishing different stackings in layered materials via luminescence spectroscopy'', [https://doi.org/10.1103/PhysRevLett.131.206902 Phys. Rev. Lett. '''131''', 206902 (2023)]; [https://arxiv.org/abs/2305.17554 arXiv 2305.17554] </ref>
<ref name='marini2024'>G. Marini, M. Calandra, P. Cudazzo, ''Optical absorption and photoluminescence of single layer boron nitride from a first principles cumulant approach'', [https://doi.org/10.1021/acs.nanolett.4c00669 Nano Lett., '''24''', 20, 6017 (2024)]; [https://arxiv.org/abs/2402.03826 arXiv 2402.03826 (2024)]</ref>
<ref name='antonius2017'>G. Antonius, S. G. Louie, ''Theory of exciton-phonon coupling'', [https://doi.org/10.1103/PhysRevB.105.085111 Phys. Rev. B, '''105''', 085111 (2022)]; [https://arxiv.org/abs/1705.04245 arXiv1705.04245 (2017)]</ref>
<ref name='paleari2022'> F. Paleari, and A. Marini, ''Exciton-phonon interaction calls for a revision of the “exciton” concept'', [https://doi.org/10.1103/PhysRevB.106.125403 Phys. Rev. B, '''106''', 125403 (2022)]; [https://arxiv.org/abs/2205.02783 arXiv 2205.02783]</ref>
<ref name='cudazzo2020'> P. Cudazzo, ''First-principles description of the exciton-phonon interaction: A cumulant approach'', [https://doi.org/10.1103/PhysRevB.102.045136 Phys. Rev. B, '''102''', 045136 (2020)]; [https://orbilu.uni.lu/bitstream/10993/44769/1/main.pdf Open access pdf from Luxembourg University]</ref>
<ref name='chan2023'>Y-h Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, ''Exciton Lifetime and Optical Line Width Profile via Exciton–Phonon Interactions: Theory and First-Principles Calculations for Monolayer MoS2'', [https://doi.org/10.1021/acs.nanolett.3c00732 Nano Lett., '''23''', 9 (2023)]; [https://arxiv.org/abs/2212.08451 arXiv 2212.08451 (2023)]</ref>
<ref name='murali2025'>M. Nalabothula, S. Reichardt, L. Wirtz, ''Origin of Interlayer Exciton–Phonon Coupling in 2D Heterostructures'', [https://doi.org/10.1021/acs.nanolett.5c00355 Nano Lett., '''25''', 15 (2025)], [https://arxiv.org/abs/2407.16111 arXiv 2407.16111 (2025)]</ref>
</references>

Shift current

2026-03-17T13:40:10Z

Wikiadmin: /* Analysis of the results */

[[File:Shift current.png|right|200px | Shift_current]]
== Introduction ==
In this tutorial we will show how to calculate Shift Current in bulk materials. 
We suppose you are already familiar with the non-linear response using the Yambo code.
If it is not the case please study the previous tutorials: [[Linear response using Dynamical Berry Phase]] and [[Real time approach to non-linear response (SHG)]]. 
This theoretical implementation of this tutorial is described in the paper <ref name="yuncheng"></ref> , please cite it if you use it. 
The tutorial was created with the help of [https://www.linkedin.com/in/yuncheng-mao-8aba6279/ Yuncheng Mao].

== Setup calculations ==

In this tutorial we will take as example the two dimensional GeS monolayer. 
DFT inputs can be downloaded here: [http://www.attaccalite.com/tutorials_yambo/GeS_DFT.tgz GeS_DFT.tgz]. 

First of all run the setup, then remove symmetries along the '''y''' direction, as explained in the tutorial above.

== Real-time setup and calculations ==

In order to generate input file for shift current you do: yambo_nl -u n -V par -F input.in 

nloptics # [R] Non-linear spectroscopy
NLogCPUs=0 # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced" # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= "8 1" # [PARALLEL] CPUs for each role
NL_ROLEs= "w k" # [PARALLEL] CPUs roles (w,k)
DIP_CPU= "" # [PARALLEL] CPUs for each role
DIP_ROLEs= "" # [PARALLEL] CPUs roles (k,c,v)
OSCLL_CPU= "" # [PARALLEL] CPUs for each role
OSCLL_ROLEs= "" # [PARALLEL] CPUs roles (k,b)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
NL_Threads=0 # [OPENMP/NL] Number of threads for nl-optics
OSCLL_Threads=0 # [OPENMP/X] Number of threads for Oscillators
% NLBands
17 | 24 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX/LSEX/LHF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 8.000000 | eV # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
NLEnSteps= 24 # [NL] Energy steps for the loop on frequencies
% NLrotaxis
0.000000 | 0.000000 | 0.000000 | # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0 # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
Field1_Freq= 0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "SOFTSIN" # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time

Notice that in this input we turned one the evaluation of current EvalCurrent and force the parallelization on the frequencies, <code>NL_CPU= "8 1"</code>, <code> NL_ROLEs= "w k</code>", that is much more efficient than the one on k-points. 
Then you can run simulation by doing: yambo_nl -F input.in

== Analysis of the results ==

In order to analyze the result we use Yambopy: [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy], and the following script for the post-processing:

from yambopy import *
from yambopy.units import fs2aut

NLDB=YamboNLDB()

SIN = Xn_from_sine(NLDB,X_order=4,l_out_current=True)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)
SIN.reconstruct_signal(OUT)

This script will produce different file containing the response function of the current and polarization respect to the total field.
From the second order response of the current at zero frequency, file <code>o.YamboPy-Sigma_probe_order_0</code> we can extract the shift current coefficient along the 'y' direction
that correspond to the 5th column of the file. In the file <code>o.YamboPy-Sigma_probe_order_0:</code>

the first column ti the energy

columns 2,3 are the imaginary and real part along x-direction

columns 4,5 are the imaginary and real part along y-direction

columns 6,7 are the imaginary and real part along y-direction

Notice that the imaginary part is zero within the numerical precision.

 
Hereafter we report the shift current spectrum compared with the results of Ibañez-Azpiroz et al. <ref name="Azpiroz"></ref>.

[[File:Sigma2.png|center| 800px | Shift current]]

The python script to plot the results, the data from the Ibañez-Azpiroz paper can be downloaded here: [https://www.attaccalite.com/tutorials_yambo/GeS_plot.tgz GeS_plot]. 
In this example we used very few k-points, increasing the number of k-points bring the result very close to the one of Ibañez-Azpiroz, see also Supp. Mat. of Ref. <ref name="yuncheng"></ref> . 
Notice that we rescale the shift-currect for an effective thickness of 2.65 Angstrom, you can find all the rescaling factor in the python script.

Correlation effects can be included in the SC response in the same way we do for other non-linear response functions, for more info look at the tutorial:
[https://wiki.yambo-code.eu/wiki/index.php?title=SHG_within_the_TD-aGW_level_(also_called_TD-HSEX,_TD-BSE) SHG within the TD-aGW level (also called TD-HSEX, TD-BSE)]

== References ==
<references>
<ref name="yuncheng">Yuncheng Mao, Myrta Grüning, Ju Zhou, and Claudio Attaccalite, ''Shift current in 2D Janus Transition-Metal Dichalcogenides: the role of excitons
'', [https://arxiv.org/abs/2506.16067 arXiv:2506.16067(2025)]</ref>
<ref name="Azpiroz">Julen Ibañez-Azpiroz, Stepan S. Tsirkin, and Ivo Souza, ''Ab initio calculation of the shift photocurrent by Wannier interpolation'', [https://doi.org/10.1103/PhysRevB.107.205204 Phys. Rev. B '''107''', 205204 (2023)]</ref>
</references>

Sum frequency generation

2026-03-06T15:37:21Z

Wikiadmin: /* Analysis of the single run */

[[File:Sum frequency.jpg|right|500px |Sum frequency generation]]

'''This tutorial is for internal use only, these response functions are not implemented/tested in yambo/yambopy suite.'''

In this tutorial we will show you how to calculate Sum Frequency Generation(SFG) and also Difference Frequency Generation (DFG) in bulk materials. 
We assume you are already familiar with the non-linear response using the Yambo code.
If it is not the case please refer to the previous tutorials: [[Linear response using Dynamical Berry Phase]] and [[Real time approach to non-linear response (SHG)]]
and Ref.<ref name=nl>[https://arxiv.org/abs/1309.4012 Nonlinear optics from an ab initio approach by means of the dynamical Berry phase: Application to second- and third-harmonic generation in semiconductors], C. Attaccalite and M. Grüning, Phys. Rev. B '''88''', 235113(2013)</ref>. 
This tutorial has been made possible thanks to the work of [https://www.uni-giessen.de/de/fbz/fb07/fachgebiete/physik/institute/theorie/agsanna/people/mike-pionteck Mike N. Pionteck] published in Ref. <ref name=sfg>[https://arxiv.org/abs/2503.07095 Sum frequency generation from real-time simulations in two-dimensional crystals], M. N. Pionteck, M. Grüning, S. Sanna, C. Attaccalite</ref>, results are also presented in this Libreoffice presentation [https://www.attaccalite.com/wp-content/uploads/2025/04/sum_frequency.odp sum_frequency.odp]. 

== DFT calculations ==
In this example, we will consider a single layer of hexagonal boron nitride (hBN).
If you haven't already done so, you can download the input files and Yambo databases for this tutorial here: [https://media.yambo-code.eu/educational/tutorials/files/hBN-2D-RT.tar.gz hBN-2D-RT.tar.gz].
and/or follow the instructions to generate the databases here: [[Prerequisites for Real Time propagation with Yambo]]

== Removing symmetries ==
In this tutorial we will calculate the SFG along when both fields are in the 'x' direction,
therefore we remove symmetries not compatible with an external field along this direction, with the command <code> ypp_nl -y</code>:

fixsyms # [R] Remove symmetries not consistent with an
external perturbation
% Efield1
 1.000000 | 0.000000 | 0.000000 | # First external Electric Field
%
% Efield2
0.000000 | 0.000000 | 0.000000 | # Additional external Electric Field
%
BField= 0.000000 T # [MAG] Magnetic field modulus
Bpsi= 0.000000 deg # [MAG] Magnetic field psi angle [degree]
Btheta= 0.000000 deg # [MAG] Magnetic field theta angle [degree]
#RmAllSymm # Remove all symmetries
RmTimeRev # Remove Time Reversal
#RmSpaceInv # Remove Spatial Inversion

== Real-time simulation with two external fields ==
You go in the <code>FixSymm</code> folder and run again the setup. Then you can put the following input file, that has been generated with the command <code>yambo_nl -u n </code>, in the folder with the name <code>yambo.in_sfg</code>:

nloptics # [R] Non-linear spectroscopy
NLogCPUs=0 # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced" # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= "10 1" # [PARALLEL] CPUs for each role
NL_ROLEs= "w k" # [PARALLEL] CPUs roles (w,k)
DIP_CPU= "" # [PARALLEL] CPUs for each role
DIP_ROLEs= "" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
NL_Threads=0 # [OPENMP/NL] Number of threads for nl-optics
% NLBands
3 | 6 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime= 70.00000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
2.000000 | 8.000000 | eV # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
NLEnSteps= 30 # [NL] Energy steps for the loop on frequencies
% NLrotaxis
0.000000 | 0.000000 | 0.000000 | # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0 # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 18475 RL # [HA] Hartree RL components
EXXRLvcs= 18475 RL # [XX] Exchange RL components
Field1_Freq= -0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= " SOFTSIN" # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time
Field2_Freq= 3.000000 eV # [RT Field2] Frequency
Field2_Int= 1000.00 kWLm2 # [RT Field2] Intensity
Field2_Width= 0.000000 fs # [RT Field2] Width
Field2_kind= " SOFTSIN" # [RT Field2] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field2_pol= "linear" # [RT Field2] Pol(linear|circular)
% Field2_Dir
 1.000000 | 0.000000 | 0.000000 | # [RT Field2] Versor
%
Field2_Tstart= 0.010000 fs # [RT Field2] Initial Time

If you run this input file with the command <code>yambo_nl -F yambo.in_sfg</code>, the code will run 30 simulations for a first laser field with frequency between 2.0 and 8.0 while the frequency of the second field is fixed at 3.0 eV. For this reason we provide you a simple python script to change also the frequency of the second field (in blue in the input). 
In the next section we will show how to analyse the non-linear response generated by the presence of two external fields.

== Analysis of the single run ==
We will use new approach to extract analyze time-dependent polarization in presence of two laser fields <ref name=sfg></ref> and used for the tutorials on second/third harmonic generation: [[Real time approach to non-linear response (SHG)]]. These new subroutine are available in the new YamboPy package. Hereafter the python script for the analysis:

from yambopy import *
from yambopy.nl.sum_frequencies import SF_Harmonic_Analysis
from yambopy.units import fs2aut
X_order=4
NLDB=YamboNLDB()
SIN2=Xn_from_freqmix(NLDB,X_order=X_order,T_range=[50.0*fs2aut,-1.0])
print(SIN2)
OUT=SIN2.perform_analysis()
SIN2.output_analysis(OUT)
SIN2.reconstruct_signal(OUT)

Notice that we set the beginning of time range where analysis is performed with the variable <code>T_range=[50.0*fs2aut,-1.0]</code>, you can also omit this variable and the code automatically will set the optimal range for the analysis using the damping given in input in the real-time simulation. 
This script will analyse the real-time response and produce all the different <math>\chi^{(n)}(\omega)</math> in the form:

o.YamboPy-X_probe_order_1_0 ==> <math>\chi^{(1)}(\omega=\omega_1)</math> 
o.YamboPy-X_probe_order_1_1 ==> <math>\chi^{(2)}(\omega=\omega_1+\omega_2)</math> 
o.YamboPy-X_probe_order_2_0 ==> <math>\chi^{(2)}(\omega=2\omega_1)</math> 
o.YamboPy-X_probe_order_0_2 ==> <math>\chi^{(2)}(\omega=2\omega_2)</math> 
o.YamboPy-X_probe_order_1_2 ==> <math>\chi^{(3)}(\omega=\omega_1+2\omega_2)</math> 
o.YamboPy-X_probe_order_2_2 ==> <math>\chi^{(4)}(\omega=2\omega_1+2\omega_2)</math> 
etc....

for example we can plot the 4 and 5 columns of the file <code>o.YamboPy-X_probe_order_1_1</code> that corresponds to the <math>\chi^{(2)}_{xyy} (\omega=\omega_1+\omega_2)</math> for <math>\omega_2 =3.0~eV</math>:

[[File:Xhi2 11.png|center | 700px |xhi11]]

Here the python script used to generate the previous figure: [https://www.attaccalite.com/tutorials_yambo/plot_xhi2_monolayer.py plot_xhi2_monolayer.py]

=== Checking the Fourier analysis===
In the previous script we turned on the flag <code>prn_Peff=True</code> that automatically will produce a series of file to check if the Fourier was successful.
The files <code> o.YamboPy-sampling_F1</code> contain the sampling points used to fix the real-time polarization, you can for example plot <code>o.YamboPy-sampling_F1</code> vs <code>o.polarization_F1</code> together to see which part of <math>P(t)</math> is used to extract the response functions:
[[File:Sampling vs p.png|center| 700px |sampling_p]]

Then you can check if real-time polarization is well reconstructed by the Fourier coefficients calculated by the script by plotting <code> o.YamboPy-pol_reconstructed_F1</code> vs <code>o.polarization_F1</code>

[[File:Reconstructed p vs p.png|center |700px |Reconstructed polarization]]

As you see the polarization is well reconstructed starting from the extracted Fourier coefficients, expect for the beginning, but this is normal because the initial signal is affected by the internal eigen--modes of the systems that later are washed out by the de-phasing. For this reason we use only the last part of the signal for the analysis, for more details see Ref. <ref name=nl></ref>

== Analysis of all runs ==
For this reason we provide you a simple python script to change also the frequency of the second field (in blue in the input).

[https://www.attaccalite.com/tutorials_yambo/run_many_sfg.py run_many_sfg.py]

please modify the python script by set the correct path of your executable, the parallelization and the number of frequency steps you are interested in. In this example we will do '''30''' frequency steps for both the first (in the Yambo input) and the second field (in the python script).

Calculations will take some time, in the above example we parallelized on 10 cores to speed up them. 
When the calculations are terminated you will find many folders called EF1, EF2, ... EF29. These folders are generated from the loop on the second field frequency(Field2_Freq), and each one contains the response real-time polarization at 30 different frequencies of the first external field (Field1_Freq).

Now we will use a script to read all the real-time polarization, Fourier analyze them as it was done in the previous section for a single one, and the write the <math>\chi^{(2)}(\omega=\omega_1 + \omega_2)</math> in a file in the form: <math>\omega_1 , \omega_2, Re \{\chi^{(2)}\} , Im \{\chi^{(2)}\}</math>.

[https://www.attaccalite.com/tutorials_yambo/analysis_sfg.py analysis_sfg.py]

then we can interpolate the results <math>\chi^{(2)}(\omega)</math> and plot results in 3D or 2D: [https://www.attaccalite.com/tutorials_yambo/plot_xhi_2d.py plot_xhi_2d.py], [https://www.attaccalite.com/tutorials_yambo/plot_xhi_3d.py plot_xhi_3d.py]

[[File:SFG 3D.png|center | 700px|SFG 3D]]

[[File:SFG 2D.png|center | 500px|SFG in 2D]]

You can compare these results with those of Ref.<ref name=sfg></ref>. Note that in order to get fully converged results, you need to increase: the number of k-points, the length of the simulations, and the number of frequencies used in the plots. Converged parameters for this system can be found in Ref.<ref name=sfg></ref>, including all details to include excitonic effects. 
'''Important''': In case of degenerate or multiple frequencies <math>\omega_1 = n \omega_2</math> where <math>n</math> is an integer, the algorithm used to extract the SFG/DFG coefficient could not converge. In this case you can discard these points and use the neighbour solution for interpolation, see discussion in reference <ref name=sfg></ref>. Another possibility is to use a least-squares fit algohrithm that starts from the neighbour solution to find the SFG/DFG coefficents, this is implemented in YamboPy in the sum_frequency.py file, but is commented.

== Comparison with SHG ==

Finally you can compare the diagonal of <math>\chi^{(2)}(\omega=\omega_1 +\omega_2)</math>, when <math>\omega_1 =\omega_2</math>, with the standard second harmonic generation <math>\chi^{(2)}(\omega=2 \omega_1)</math> calculated as explained in the tutorial [[Real time approach to non-linear response (SHG)]]:

[[File:Xhi2 vs xhi2.png|center | 700px|xhi2_vs_xhi2]]

Notice that there is factor two between the two response functions, this is the so-called degeneracy factor, for more details see Sec. 1.3 in Ref.<ref>Boyd, Robert W., Alexander L. Gaeta, and Enno Giese. "Nonlinear optics." Springer Handbook of Atomic, Molecular, and Optical Physics. Cham: Springer International Publishing, 1097-1110.(2008)</ref>. The small differences between the two spectra are given the different sampling strategies of the <code>SF_Harmonic_Analysis</code> function respect to the <code>Harmonic_Analysis</code> one.

== Reference ==

Shift current

2026-03-06T14:31:01Z

Wikiadmin: /* Analysis of the results */

[[File:Shift current.png|right|200px | Shift_current]]
== Introduction ==
In this tutorial we will show how to calculate Shift Current in bulk materials. 
We suppose you are already familiar with the non-linear response using the Yambo code.
If it is not the case please study the previous tutorials: [[Linear response using Dynamical Berry Phase]] and [[Real time approach to non-linear response (SHG)]]. 
This theoretical implementation of this tutorial is described in the paper <ref name="yuncheng"></ref> , please cite it if you use it. 
The tutorial was created with the help of [https://www.linkedin.com/in/yuncheng-mao-8aba6279/ Yuncheng Mao].

== Setup calculations ==

In this tutorial we will take as example the two dimensional GeS monolayer. 
DFT inputs can be downloaded here: [http://www.attaccalite.com/tutorials_yambo/GeS_DFT.tgz GeS_DFT.tgz]. 

First of all run the setup, then remove symmetries along the '''y''' direction, as explained in the tutorial above.

== Real-time setup and calculations ==

In order to generate input file for shift current you do: yambo_nl -u n -V par -F input.in 

nloptics # [R] Non-linear spectroscopy
NLogCPUs=0 # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced" # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= "8 1" # [PARALLEL] CPUs for each role
NL_ROLEs= "w k" # [PARALLEL] CPUs roles (w,k)
DIP_CPU= "" # [PARALLEL] CPUs for each role
DIP_ROLEs= "" # [PARALLEL] CPUs roles (k,c,v)
OSCLL_CPU= "" # [PARALLEL] CPUs for each role
OSCLL_ROLEs= "" # [PARALLEL] CPUs roles (k,b)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
NL_Threads=0 # [OPENMP/NL] Number of threads for nl-optics
OSCLL_Threads=0 # [OPENMP/X] Number of threads for Oscillators
% NLBands
17 | 24 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX/LSEX/LHF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 8.000000 | eV # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
NLEnSteps= 24 # [NL] Energy steps for the loop on frequencies
% NLrotaxis
0.000000 | 0.000000 | 0.000000 | # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0 # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
Field1_Freq= 0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "SOFTSIN" # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time

Notice that in this input we turned one the evaluation of current EvalCurrent and force the parallelization on the frequencies, <code>NL_CPU= "8 1"</code>, <code> NL_ROLEs= "w k</code>", that is much more efficient than the one on k-points. 
Then you can run simulation by doing: yambo_nl -F input.in

== Analysis of the results ==

In order to analyze the result we use Yambopy: [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy], and the following script for the post-processing:

from yambopy import *
from yambopy.nl.harmonic_analysis import Harmonic_Analysis
from yambopy.units import fs2aut

NLDB=YamboNLDB()

SIN = Xn_from_sine(NLDB,X_order=4,l_out_current=True)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)
SIN.reconstruct_signal(OUT)

This script will produce different file containing the response function of the current and polarization respect to the total field.
From the second order response of the current at zero frequency, file <code>o.YamboPy-Sigma_probe_order_0</code> we can extract the shift current coefficient along the 'y' direction
that correspond to the 5th column of the file. In the file <code>o.YamboPy-Sigma_probe_order_0:</code>

the first column ti the energy

columns 2,3 are the imaginary and real part along x-direction

columns 4,5 are the imaginary and real part along y-direction

columns 6,7 are the imaginary and real part along y-direction

Notice that the imaginary part is zero within the numerical precision.

 
Hereafter we report the shift current spectrum compared with the results of Ibañez-Azpiroz et al. <ref name="Azpiroz"></ref>.

[[File:Sigma2.png|center| 800px | Shift current]]

The python script to plot the results, the data from the Ibañez-Azpiroz paper can be downloaded here: [https://www.attaccalite.com/tutorials_yambo/GeS_plot.tgz GeS_plot]. 
In this example we used very few k-points, increasing the number of k-points bring the result very close to the one of Ibañez-Azpiroz, see also Supp. Mat. of Ref. <ref name="yuncheng"></ref> . 
Notice that we rescale the shift-currect for an effective thickness of 2.65 Angstrom, you can find all the rescaling factor in the python script.

Correlation effects can be included in the SC response in the same way we do for other non-linear response functions, for more info look at the tutorial:
[https://wiki.yambo-code.eu/wiki/index.php?title=SHG_within_the_TD-aGW_level_(also_called_TD-HSEX,_TD-BSE) SHG within the TD-aGW level (also called TD-HSEX, TD-BSE)]

== References ==
<references>
<ref name="yuncheng">Yuncheng Mao, Myrta Grüning, Ju Zhou, and Claudio Attaccalite, ''Shift current in 2D Janus Transition-Metal Dichalcogenides: the role of excitons
'', [https://arxiv.org/abs/2506.16067 arXiv:2506.16067(2025)]</ref>
<ref name="Azpiroz">Julen Ibañez-Azpiroz, Stepan S. Tsirkin, and Ivo Souza, ''Ab initio calculation of the shift photocurrent by Wannier interpolation'', [https://doi.org/10.1103/PhysRevB.107.205204 Phys. Rev. B '''107''', 205204 (2023)]</ref>
</references>

Shift current

2026-03-06T14:30:24Z

Wikiadmin: /* Analysis of the results */

[[File:Shift current.png|right|200px | Shift_current]]
== Introduction ==
In this tutorial we will show how to calculate Shift Current in bulk materials. 
We suppose you are already familiar with the non-linear response using the Yambo code.
If it is not the case please study the previous tutorials: [[Linear response using Dynamical Berry Phase]] and [[Real time approach to non-linear response (SHG)]]. 
This theoretical implementation of this tutorial is described in the paper <ref name="yuncheng"></ref> , please cite it if you use it. 
The tutorial was created with the help of [https://www.linkedin.com/in/yuncheng-mao-8aba6279/ Yuncheng Mao].

== Setup calculations ==

In this tutorial we will take as example the two dimensional GeS monolayer. 
DFT inputs can be downloaded here: [http://www.attaccalite.com/tutorials_yambo/GeS_DFT.tgz GeS_DFT.tgz]. 

First of all run the setup, then remove symmetries along the '''y''' direction, as explained in the tutorial above.

== Real-time setup and calculations ==

In order to generate input file for shift current you do: yambo_nl -u n -V par -F input.in 

nloptics # [R] Non-linear spectroscopy
NLogCPUs=0 # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced" # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= "8 1" # [PARALLEL] CPUs for each role
NL_ROLEs= "w k" # [PARALLEL] CPUs roles (w,k)
DIP_CPU= "" # [PARALLEL] CPUs for each role
DIP_ROLEs= "" # [PARALLEL] CPUs roles (k,c,v)
OSCLL_CPU= "" # [PARALLEL] CPUs for each role
OSCLL_ROLEs= "" # [PARALLEL] CPUs roles (k,b)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
NL_Threads=0 # [OPENMP/NL] Number of threads for nl-optics
OSCLL_Threads=0 # [OPENMP/X] Number of threads for Oscillators
% NLBands
17 | 24 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX/LSEX/LHF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 8.000000 | eV # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
NLEnSteps= 24 # [NL] Energy steps for the loop on frequencies
% NLrotaxis
0.000000 | 0.000000 | 0.000000 | # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0 # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
Field1_Freq= 0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "SOFTSIN" # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time

Notice that in this input we turned one the evaluation of current EvalCurrent and force the parallelization on the frequencies, <code>NL_CPU= "8 1"</code>, <code> NL_ROLEs= "w k</code>", that is much more efficient than the one on k-points. 
Then you can run simulation by doing: yambo_nl -F input.in

== Analysis of the results ==

In order to analyze the result we use Yambopy: [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy], and the following script for the post-processing:

from yambopy import *
from yambopy.nl.harmonic_analysis import Harmonic_Analysis
from yambopy.units import fs2aut
X_order=4
NLDB=YamboNLDB()
SIN = Xn_from_sine(NLDB,l_out_current=True)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)
SIN.reconstruct_signal(OUT)

This script will produce different file containing the response function of the current and polarization respect to the total field.
From the second order response of the current at zero frequency, file <code>o.YamboPy-Sigma_probe_order_0</code> we can extract the shift current coefficient along the 'y' direction
that correspond to the 5th column of the file. In the file <code>o.YamboPy-Sigma_probe_order_0:</code>

the first column ti the energy

columns 2,3 are the imaginary and real part along x-direction

columns 4,5 are the imaginary and real part along y-direction

columns 6,7 are the imaginary and real part along y-direction

Notice that the imaginary part is zero within the numerical precision.

 
Hereafter we report the shift current spectrum compared with the results of Ibañez-Azpiroz et al. <ref name="Azpiroz"></ref>.

[[File:Sigma2.png|center| 800px | Shift current]]

The python script to plot the results, the data from the Ibañez-Azpiroz paper can be downloaded here: [https://www.attaccalite.com/tutorials_yambo/GeS_plot.tgz GeS_plot]. 
In this example we used very few k-points, increasing the number of k-points bring the result very close to the one of Ibañez-Azpiroz, see also Supp. Mat. of Ref. <ref name="yuncheng"></ref> . 
Notice that we rescale the shift-currect for an effective thickness of 2.65 Angstrom, you can find all the rescaling factor in the python script.

Correlation effects can be included in the SC response in the same way we do for other non-linear response functions, for more info look at the tutorial:
[https://wiki.yambo-code.eu/wiki/index.php?title=SHG_within_the_TD-aGW_level_(also_called_TD-HSEX,_TD-BSE) SHG within the TD-aGW level (also called TD-HSEX, TD-BSE)]

== References ==
<references>
<ref name="yuncheng">Yuncheng Mao, Myrta Grüning, Ju Zhou, and Claudio Attaccalite, ''Shift current in 2D Janus Transition-Metal Dichalcogenides: the role of excitons
'', [https://arxiv.org/abs/2506.16067 arXiv:2506.16067(2025)]</ref>
<ref name="Azpiroz">Julen Ibañez-Azpiroz, Stepan S. Tsirkin, and Ivo Souza, ''Ab initio calculation of the shift photocurrent by Wannier interpolation'', [https://doi.org/10.1103/PhysRevB.107.205204 Phys. Rev. B '''107''', 205204 (2023)]</ref>
</references>

Third Harmonic Generation (THG)

2026-03-06T14:28:42Z

Wikiadmin: /* Analysis of the results with YamboPy */

==Third Harmonic Generation in bulk Silicon==

In this tutorial we will calculate third harmonic generation (THG) in bulk Silicon.
[[File:Thg silicon.jpg|center| 150px]]

 Calculation of third harmonic generation proceeds in similar way to the SHG calculations. DFT input can be found here ([http://www.attaccalite.com/lumen/tutorials/silicon_dft.tgz abinit] or [http://www.attaccalite.com/lumen/tutorials/silicon_qe.tgz QuantumEspresso]). We import the wave-functions as explained in the [[Real time approach to non-linear response]], then perform the setup using only 1000 plane waves. 

==Removing symmetries==
We are interested in the <math> \chi^{(3)}_{1111} </math> therefore we consider an external field in direction [1, 0, 0] (equivalent to [x,0,0]). 
We remove all the symmetries not compatible with the external field plus the time-reversal symmetry by doing <code>ypp_nl (or ypp) -y</code>: 
 
fixsyms # [R] Reduce Symmetries
% Efield1
1.00 | 0.00 | 0.00 | # First external Electric Field
%
% Efield2
0.00 | 0.00 | 0.00 | # Additional external Electric Field
%
#RmAllSymm # Remove all symmetries
RmTimeRev # Remove Time Reversal

==Real-time simulations==
Then you go in the FixSymm folder, run again the setup (as explained in the [[Linear response using Dynamical Berry Phase]]) and then generate the input file with the command <code>yambo_nl -u n -V qp -F input.in</code>: 

nloptics # [R NL] Non-linear optics
% NLBands
1 | 7 | # [NL] Bands
%
NLstep= 0.0100 fs # [NL] Real Time step length
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "CRANKNIC" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/HF/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
0.200000 | 2.000000 | eV # [NL] Energy range
%
NLEnSteps= 6 # [NL] Energy steps
NLDamping= 0.200000 eV # [NL] Damping
%Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [RT Field1] Versor
Field1_kind= "SOFTSIN" # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_Int= 0.1000E+8 kWLm2 # [NL ExtF] Intensity
% GfnQP_E
0.600000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v) eV|adim|adim
%
We calculate the <math>\chi^{(3)}_{1111}</math> only on 6 frequencies, but calculations can be extended to larger frequency range. We introduce a scissor operator of 0.6 eV in such a way to reproduce the gap of bulk silicon. Notice that we increased the field intensity in such a way to improve the ratio between non-linear response and numerical noise. 
Now that you have created the file <code>input.in</code> you can run the non-linear optics calculation with the command: 
 

<pre>yambo_nl -F input.in</pre>
 
when the simulation will end, it will have produced the file <code>SAVE/ndb.Nonlinear</code> with all information on the simulation and the files <code>SAVE/ndb.Nonlinear_fragment_x</code> with the polarization and current for each laser frequency. 

==Analysis of the results==
Results can be analyzed in the same way of the SHG, using the command <code>ypp_nl -u</code>: 
 

nonlinear # [R] NonLinear Optics Post-Processing
Xorder= 5 # Max order of the response functions
% TimeRange
40.000000 | -1.00000 | fs # Time-window where processing is done
%
ETStpsRt= 200 # Total Energy steps
% EnRngeRt
0.00000 | 10.00000 | eV # Energy range
%
DampMode= "NONE" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.10000 eV # Damping parameter

This time we increase the de-phasing time to 40 fs in such a way to have a clean third-harmonic response. 
Hereafter we report the result with the 8x8x8 k-points sampling and the converged result with the 24x24x24 k-points sampling:
[[File:Si THG results.png|center| 600px]]
Plot data can be downloaded [http://www.attaccalite.com/lumen/tutorials/Si_thg_results.tgz here].

==Analysis of the results with YamboPy==
'''This part works only with Yambo 5.3 or the last version available on Github:''' [https://github.com/yambo-code/yambo https://github.com/yambo-code/yambo]
 
 

The previous analysis can be performed using <code>YamboPy</code>, see these other tutorials for more details ([[Real time approach to non-linear response#Analysis of the results using YamboPy]])
From the folder of the calculations in python you use the script:

from yambopy import *
from yambopy.plot import *
from yambopy.units import fs2aut

NLDB=YamboNLDB()

pol =NLDB.Polarization[0]
time=NLDB.IO_TIME_points

SIN = Xn_from_sine(NLDB,X_order=5,T_range=[40.0*fs2aut,-1.0])
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)

it will produce files with the non-linear response ''o.YamboPY-X_probe_order_x'' equivalent to the ones produced by <code>ypp_nl</code>.

==How to choose the external field direction for the different X3==
In the X3abcd response you have four field directions. The three components ''b'',''c'' and ''d'' are determined by the external field and the four ''a''
is the direction where you measure the response. 

The external field in Yambo is set using the variable
%Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [NL ExtF] Versor .
%

For example if you set in Yambo field direction as:

1.0 | 0.0 | 0.0 | ---> you calculate all components X3*xxx
0.0 | 1.0 | 0.0 | ---> you calculate all components X2*yyy
0.0 | 0.0 | 1.0 | ---> you calculate all components X2*zzz

the * correspond to the different direction in the response file <code>o.YPP-X_probe_order_3</code>.
In the <code>o.YPP-X_probe_order_3</code> you have seven columns:

The first column is the energy [[omega]],
the 2nd and 3rd columns are the imaginary and real part of X3x***
the 4th and 5th columns are the imaginary and real part of X3y***
the 6th and 7th columns are the imaginary and real part of X3z***

where the *** are the direction of the incoming field, see above.

Instead for an external field in Yambo in the direction
%Field1_Dir
1.000000 | 1.000000 | 1.000000 | # [NL ExtF] Versor .
%

you can get in output in the <code>o.YPP-X_probe_order_3</code> file the components:

the 2nd and 3rd columns are the imaginary and real part of X3xxyz
the 4th and 5th columns are the imaginary and real part of X3yxyz
the 6th and 7th columns are the imaginary and real part of X3zxyz

Other terms of the X3 tensor can be obtained combining the different columns of the <code>o.YPP-X_probe_order_3</code> file and using external fields along linear combination of the Cartesian axes.

'''Nota bene:''' if you change the direction of the external field you have to remove the corresponding symmetries as explained in the tutorial
on [[Prerequisites for Real Time propagation with Yambo#Reduce symmetries real-time linear response tutorial]].

Real time approach to non-linear response (SHG)

2026-03-06T14:25:18Z

Wikiadmin: /* Analysis of the results using YamboPy */

== Second Harmonic Generation in AlAs ==
In this tutorial, we will calculate the second harmonic generation of bulk AlAs, the Lumen databases can be downloaded here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_DBs.tar.gz AlAs_DBs.tar.gz] (10 MB). The first steps of this tutorials are the same as the one on [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo] with the only difference that we will work on AlAs and we will consider an external field in the direction <code> 1.000000 | 1.000000 | 0.000000 |</code>. The DFT input are available here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_abinit.tar.gz ABINIT] or [https://media.yambo-code.eu/educational/tutorials/files/AlAs_pwscf.tar.gz QuantumEspresso].
You can run the DFT calculation with ABINIT with the command:

abinit < AlAs.in > output_AlAs

or using QuantumEspresso:

pw.x -inp AlAs.scf.in > output_scf
pw.x -inp AlAs.nscf.in > output_scf

and then import the ABINIT wave-function with the command:

a2y -F AlAso_DS2_WFK.nc

and the QuantumEspresso one, in the folder AlAs.save with the command:

p2y
Now we consider an external field in the <code>[1,1,0]</code> direction and remove symmetries not compatible with this field, as explained in the tutorial [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo].

==Real-time simulation for the SHG==

You can generate the input file with the command <code>yambo_nl -u n</code>:

nloptics # [R NL] Non-linear optics
% NLBands
 3 | 6 | # [NL] Bands
%
NLstep= 0.0100 fs # [NL] Real Time step length
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/SEX/HF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 5.000000 | eV # [NL] Energy range
%
NLEnSteps= 5 # [NL] Energy steps
NLDamping= 0.150000 eV # [NL] Damping
% Field1_Dir
 1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
Field1_kind= "SOFTSIN" # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

The line <code>|1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor</code> referees to the direction of the external field (x,y,0). 
The default parameters of Lumen are already tuned for second-harmonic generation, so the only thing you have to change is the band range, between 3 and 6 and the energy range between 1.0-5.0 eV and the number of energy steps in this interval that we set to 10. Notice that you cannot set to zero the lowest value of the energy range because this will requires a simulation that lasts infinite time, see below. Finally, consider that Yambo performs a separate calculation for each frequency, so if you set many energy steps the computational time grows linearly with this number. 

Notice that we set <code>NLverbosity= "high"</code> in this way the code will produce a file for each laser frequency containing the time dependent polarization. 

Run <code>yambo_nl</code>. Calculations will take about fifteen minutes on a single processor PC, you have time to study the next section that explains how non-linear response is extracted from the real-time simulations.
In order to speed up calculations, you can run them in parallel.

==Non-linear response with Yambo==

In order to calculate the non-linear response, the system is excited with different laser fields with a sinusoidal shape at frequencies <math>\omega_1, \omega_2, .... , \omega_n </math> determined by the parameters [[Variables#NLEnRangev|NLEnRange]] and [[Variables#NLEnSteps|NLEnSteps]]. A dephasing term is added to the Hamiltonian <math>\gamma = </math>[[Variables#NLDamping|NLDamping]] to simulate a finite broadening and to remove the eigenmodes that are excited by the sudden turning on of the external field. After the dephasing time the outgoing signal is analyzed to extract the non-linear coefficients as shown in the figure below:

[[File:Pt analysis.png|600px|center|Non-linear response analysis]]

In the report file <code>r_nlinear</code> you will find the length of the dephasing part and of the sampling one:

Dephasing Time [fs]: 52.65695
Sampling Time [fs]: 4.18566
Total simulation time [fs]: 56.84262
The length of the dephasing interval is inversely proportional to the damping term <math>T_{depth} \simeq 1/\gamma</math> while the length of the sampling is dictated by the smallest frequency we are interested in: <math>T_{samp} \simeq 1/\omega_1 </math>. For this reason if <math> \omega_1=0 </math> or <math> \gamma = 0 </math> simulation time goes to infinity. The response at zero frequency can be calculated as limit of small frequency perturbation.
You can have a look to the file <code>src/nloptics/NL_initialize.F</code> to see how simulation lengths are defined.
Calculations can take some time, run them in parallel, the best number of processors = number of frequencies or reduce the number of frequencies step in the SHG.

==Analysis of the results==

In the sampling region we suppose that the polarization can be written as: <math>\bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{p}_n e^{-i\omega_n t} </math> where the coefficient <math>\bf{p}_1,...,\bf{p}_n </math> are related to <math>\chi^{(1)},...,\chi^{(n)} </math>. We sample the polarization signal at different times and invert the previous equation by truncating the sum at a finite order <ref>[https://arxiv.org/abs/1309.4012 Nonlinear optics from an ab initio approach by means of the dynamical Berry phase: Application to second- and third-harmonic generation in semiconductors], C. Attaccalite and M. Grüning, Phys. Rev. B 88, 235113(2013)</ref>.
This is done with the command <code>ypp_nl -u</code> that automatically produce an input with the correct values:

nonlinear # [R] NonLinear Optics Post-Processing
Xorder= 4 # Max order of the response functions
% TimeRange
52.65695 | -1.00000 | fs # Time-window where processing is done
%
ETStpsRt= 200 # Total Energy steps
% EnRngeRt
0.00000 | 10.00000 | eV # Energy range
%
DampMode= "NONE" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.10000 eV # Damping parameter

where <code>Xorder</code> is the order where the previous sum is truncated, and the TimeRange specifies the sampling region. Notice that differently from the first tutorial, in this case, we do not need Damping in <code>ypp_nl</code> because we already included it in the real-time dynamics. Run <code>ypp_nl</code> and it will produce a file called <code>o.YPP-X_probe_order_2</code>. This file contains the different components of <math>\chi^2_{xy*} </math> and in particular the columns 6 and 7 correspond to the imaginary and real part of <math>\chi^2_{xyz} </math>. You can plot the result with gnuplot and the command:

p 'o.YPP-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

the result should look like

[[File:AlAs results.png|700px|center|SHG in AlAs]]

and compare it with the <math> \chi^2_{zxy} </math> calculated with more frequencies and with the converged result on a <math> 18x18x18</math> k-point grid and bands between 2 and 10. In the figure you can find also the comparison with the results of Ref. <ref>Luppi et al. [https://arxiv.org/abs/1006.2649 PRB B, '''82''', 235201(2010)]</ref>.
Notice that the QuantumEspresso results are slightly different from the Abinit ones. This is due to the different pseudo-potential employed in the calculations if pseudopotentials were the same calculation would have been identical. You can download the script to generate this plot and the converged results [http://www.attaccalite.com/lumen/tutorials/AlAs/AlAs_results.tgz here].

'''Notice that''' : the following parameters are not used in the non-linear response analysis: <code> EnRngeRt</code>, <code>ETStpsRt</code>, <code>DampMode</code>, <code>DampFactor</code>. The damping factor is set directly in the real-time simulation.

'''System of units in Non-linear Optics''': notice that in Yambo the non-linear response functions are in the gaussian system of units, for example the <math> \chi^2 (\omega) </math> is in cm/statvolt. In order to convert the non-linear coefficients between different system of units, you can have a look at the [https://www.attaccalite.com/tutorials_yambo/Appendix_C_Boyd.pdf Appendix C] of the book [https://www.elsevier.com/books/nonlinear-optics/boyd/978-0-12-369470-6 "Non-linear optics" by Robert W. Boyd].

More tutorials on non-linear response (Third harmonic generation, parallelization, spin-orbit etc..) with Yambo can be found here: [http://www.attaccalite.com/lumen/tutorials.html More Tutorials]. 

==Analysis of the results using YamboPy==
The analysis of the result performed in the previous section using <code>ypp_nl</code> can be performed using python script [https://github.com/yambo-code/yambopy YamboPy]. 
In YamboPy we implemented all necessary functions to reads the Yambo databases and post-process the results. Here we show how to get the SHG signal from the previous
simulations. We suppose you correctly installed <code>YamboPy</code> on your PC. Go in the folder where you ran non-linear calculation and type in the python:
 
from yambopy import *
NLDB=YamboNLDB()
pol =NLDB.Polarization[0]
time=NLDB.IO_TIME_points

in this way you read all the Non-linear databases. If your runs are in a different folder then the 'SAVE' one you can specify it using the command:
NLDB=YamboNLDB(calc='MYJOB')

now in the array <code>pol</code> you have all the polarization for all laser frequencies in the three Cartesian directions, while the variable <code>time</code> contains all the time series of your simulation.
Now you can get the non-linear response with the command:

SIN = Xn_from_sine(NLDB,X_order=5)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)

this command will perform a Fourier analysis of the results in the same way of <code>ypp_nl</code> and generate new files with all the requested harmonics ''o.YamboPy-X_probe_order_1'', ''o.YamboPy-X_probe_order_2'' etc...
equivalent to the one generated by <code>ypp_nl</code>. Results can be plot with gnuplot using the command:

p 'o.YamboPy-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

The script that performs non-linear analysis can be found in <code>yambopy/nl/harmonic_analysis.py</code>. We strongly advice you to have a look to this script and modify it according
to your needs and in order to extract other non-linear response functions.

==How to choose the external field direction for the different X2==

In the X2abc response you have three field directions. The two components ''b'' and ''c'' are determined by the external field and the third ''a''
is the direction where you measure the response. 

The external field in Yambo is set using the variable
% Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [NL ExtF] Versor .
%

For example if you set in Yambo field direction as:

1.0 | 0.0 | 0.0 | ---> you calculate all components X2*xx
0.0 | 1.0 | 0.0 | ---> you calculate all components X2*yy
0.0 | 0.0 | 1.0 | ---> you calculate all components X2*zz

the * correspond to the different direction in the response file <code>o.YPP-X_probe_order_2</code>.
In the <code>o.YPP-X_probe_order_2</code> you have seven columns:

The first column is the energy [[omega]],
the 2nd and 3rd columns are the imaginary and real part of X2x**
the 4th and 5th columns are the imaginary and real part of X2y**
the 6th and 7th columns are the imaginary and real part of X2z**

where the ** are the direction of the incoming field, see above.

For example if you want to calculate the <math>\chi^{(2)}_{zxx}</math> you put the external field in the direction <code> 1.0 | 0.0 | 0.0 |</code>
and then plot the column 6 and 7 of the <code>o.YPP-X_probe_order_2</code> file. 

In order to get the other components of the X2 the calculation are more involved. For example if you put field in the

1.0 | 1.0 | 0.0 |

this corresponds to <math>E(t)=\frac{1}{\sqrt{2}} [ E_x(t) + E_y(t) ]</math> and so at the second order you will have <math>E^2(t)=\frac{1}{2} \left [ E^2_x(t) + E^2_y(t) + E_x(t) E_y(t) + E_y(t) E_x(t) \right]</math> 
this means that the X2 in the <code>o.YPP-X_probe_order_2</code> will correspond to a linear combination of the following different response functions. 
For example in the the 2nd and 3rd columns are the imaginary and real part of: 

<math>\chi^{(2)}_{x**}=\frac{1}{2} \left [\chi^2_{xxx}+\chi^{(2)}_{xyy} + \chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right ]</math>

the same for the other column replacing the first <code>x</code> for <code>y,z</code>.
If you know that some components are zero this allow you to extract the corresponding off-diagonal part of <math>\chi^{(2)}_{abc}</math> otherwise you have to perform different calculations to get the component you are interested in. 
Now we present '''a simple example of monolayer hBN (m-hBN)'''. In m-hBN there is only a non-zero component of the <math>\chi^{(2)}</math>and all the others are related by symmetry: 
<math>\chi^{(2)}_{yyy}=-\chi^{(2)}_{yxx}=-\chi^{(2)}_{xxy}=-\chi^{(2)}_{xyx}</math> 
so if you put the field in the <code>|1|1|0|</code> direction you will get: 
<math>\chi^{(2)}_{x**} = \frac{1}{2} \left[\chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right] = - \chi^{(2)}_{yyy} </math> 
<math>\chi^{(2)}_{y**} =\frac{1}{2} \left[ \chi^{(2)}_{yyy} + \chi^{(2)}_{yxx} \right] = 0.0 </math>
 
In the figure below you can see the different <math>\chi^{(2)}</math> of monolayer hBN.

[[File:Xhi2 xy.png|center| 800px|xhi2_xy]]

as you can see in the second panel, when the field is in the <code>xy</code> you get a linear combination of <math>\chi^{(2)}</math> corresponding to the equation above, in this case the <code>y</code> is zero while the <code>x</code> component is equal to <math>-\chi^{(2)}_{yyy}</math>.

'''Nota bene:''' if you change the direction of the external field you must to remove the corresponding symmetries as explained in the tutorial [[Prerequisites for Real Time_propagation with Yambo#Reduce symmetries real-time linear response tutorial]]. If you have to calculate <math>\chi^{(2)}</math> along many directions it can be convenient to remove all symmetries at once.

==Non-linear response of low-dimensional structures: 2D or 1D ==

The Yambo code is a 3D periodic code, so when you want to study a low-dimensional system 2D or 1D, you have to use a supercell approach.
For example to simulate the non-linear response in a 2D-crystal one, can put the crystal in the xy plane and choose a larger distance in the z-direction
in order to reduce the interaction between the periodic replica. Clearly in this case k-point sampling will be only in the x,y directions.
For calculations on this type of systems, these two warnings must be taken into account:

* if you calculate GW corrections or include electron-hole interaction in the linear/non-linear response it is a good idea to use a '''Coulomb cutoff''', similar to the 2D BSE case [http://www.yambo-code.org/wiki/index.php?title=How_to_treat_low_dimensional_systems How to treat low dimensional systems]. You can add the cutoff just adding the "-r" in the input generation.

* the calculation of SHG and THG are always performed respect to the supercell, therefore you have to '''re-scale''' the result to the effective thickness of the layer, usually the inter-layer distance of the corresponding bulk material for 2D systems: <math>\chi_{rescaled}(\omega) = L_z/d_{eff} \cdot \chi(\omega) </math> where <math>L_z</math> is the z-dimension of the supercell, and <math>d_{eff}</math> is the effective thickness of the 2D system.

'''Nota bene''': Non-linear response along non-periodic directions is not implemented in Yambo, the code will print zero in these directions for all response functions. We plan to implement it in a future version of the code.

==References==
<references />

==Links for schools==
* Back to [[Rome 2023#Tutorials]]
* Back to [[ICTP2020]]

 
{| style="width:100%" border="1"
|style="width:25%; text-align:left"|Prev: [[Linear response from real time simulations]]
|style="width:45%; text-align:center"|Now: [[ICTP2020|ICTP Tutorials]] --> [[Real time approach to non-linear response|Non Linear Response]]
|style="width:50%; text-align:left"|Next: [[Correlation effects in the non-linear response]]
|-
|}

Main Page

2026-03-04T17:06:23Z

Wikiadmin:

Lumen wiki

* [[Installation]]
* [[Tutorials]]<nowiki/>
* [https://enccs.github.io/max-coe-workshop/day5-tutorials/ Yambo tutorials on Markdown]
* [[Lumen release notes|Lumen Releases]]
* [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy]
* [[Configure library externally]]
* [https://www.lumen-code.org/ford/index.html Ford documentation]
* [[Internal wiki]]
* [[Test-suite]]

Magneto-electric reponse

2026-02-16T09:44:25Z

Wikiadmin: Created page with "Sum frequency generation In this tutorial we show step by step how to calculate magneto-electric coefficents using the Lumen code. '''This tutorial has been supported by EU funding via the TIMES MASC-DN''' == References =="

[[File:Quantizedmag.jpg|right|300px |Sum frequency generation]]

In this tutorial we show step by step how to calculate magneto-electric coefficents using the Lumen code.

'''This tutorial has been supported by EU funding via the TIMES MASC-DN'''
== References ==

2026-01-29T14:16:11Z

Wikiadmin:

2026-01-29T14:10:55Z

Wikiadmin: