Self-trapped polarons & Variational Polaron Equations¶

This tutorial provides a guide on using the Variational Polaron Equations framework to compute self-trapped polarons in semiconductors and insulators.

We start with a brief overview of the self-trapped polaron formation problem, followed by an introducion to the theoretical formalism of the Variational Polaron Equations. Then we demonstrate the practical aspects of computation of various polaronc properties using small hole and large electron polarons in LiF as examples. We also showcase the capabilities of AbiPy to post-process the results and perform typical convergence studies.

It is assumed the user has already completed the two tutorials RF1 and RF2, and that they are familiar with the calculation of ground state (GS) and response properties in particular phonons, Born effective charges and the high-frequency dielectric tensor. The user should have read the introduction tutorial for the EPH code before running these examples.

This lesson should take about 2 hours.

Formalism¶

Self-trapped polarons¶

In many condensed matter systems, interaction with phonons leads to the formation of polarons, when a particle couples to the lattice vibrations of a system. Instead of a free particle, it becomes a quasiparitcle with altered physical properties – a polaron. In general, this notion applies to any particle irrespective of its type and statistics: electrons, holes, excitions, magnons and so on. In the present formalism, we treat only electron and hole polarons, which formation is mediated by the electron-phonon coupling.

Depending on the strenght of the electron-phonon coupling, one can distinguish different regimes. In the weak-coupling limit, a charge carrier is coherent with phonon oscillations, dragging the associated phonon cloud. This primarily leads to the band renormalization, and the physics can be captured by the ZPR calculations. In the limit of strong-coupling, the picture is different: a carrier is trapped in the potential well, created by lattice deformation – a phenomena called self-trapping or autolocalization. This leads to the occurence of a localized polaronic state inside the bandgap, and may switch the transport regime from a band-like to hopping transport. The intermediate regime between the strong-coupling and the weak-coupling limit is the most challenging to address and requires advanced techniques. In the present formalism, we treat only the strong-coupling regime of polaron formation, i.e. self-trapped polarons.

The probabilty of a self-trapped electron or hole polaron formation can be described by its binding energy:

\[\begin{equation} \label{eq:epol} E_\mathrm{pol} = E\left( N \pm 1, \Delta \mathbf{\tau} \right) - E\left( N \pm 1, \Delta \mathbf{\tau} \equiv 0 \right), \end{equation}\]

which is the difference of total energy of a charged system with $N \pm 1$ electrons in the distorted polaronic geometry and the pristine one. Here, the lattice distortion is denoted by $\Delta \tau$ – collective atomic displacements from their equilibrium positions. A stable polaronic configuration corresponds to $E_\mathrm{pol} < 0$.

Schematically, it can be illustrated by the following configurational diagram for total energy of a system with $N\pm1$ electrons as a function of lattice distortion.

This figure also contains some foreshadowing as it uses LiF to illustrate both configurations of an extra electron in the system: delocalized and polaronic one. By the end of the tutorial, the user will learn how to generate such figures themselves.

From the configurational diagram, it is seen that $E_\mathrm{pol}$ can be decomposed into two contributions: the phonon energy $E_\mathrm{ph}$ arising purely from the lattice distortion, and the charge localization energy $\varepsilon_\mathrm{loc}$, which accounts for both electronic contribution and electron-phonon coupling. Physically, the latter term corresponds to the energy of the localized polaronic state inside the bandgap, computed from the VBM/CBM in case of hole/electron polaron.

During the tutorial, the aformentioned quantities and other properties, characterizing the polaron formation, are computed within a Variational Polaron Equations, implemented in ABINIT.

Variational Polaron Equations¶

In DFT, the common approach to the polaronic problem is a direct geometry relaxation of a charged system with a compensating jellium background. This approach is powerful and straightforward, but requires supercell calculations, making it prohibitive for polarons with large spatial extension. Also, standard DFT favors delocalization due the self-interaction, and self-interaction corrections and special approaches have to be used to localize polarons [Falletta2022].

However, it possible to avoid the need of charged supercell relaxation by reformulating the problem of polaron formation in reciprocal space. This formalism was presented by Sio and coworkers in 2019 [Sio2019], and required the solution of an effective polaronic Hamiltionian, parametrized from first-principles calculations. For the underlying theory behind this formalism and derivations, the user may be interested to consult the original paper.

We take on the approach presented by Sio and coworkers, albeit we reformulate it into a variational form. We refer to this formalism as Variational Polaron Equations or, simply, VarPEq. The derivations of our formalism and relevant discussions can be found in [Vasilchenko2022]. In this approach, we utilize a variational expression for the polaron binding energy, enabling the application of efficient gradient-based optimization techniques. Moreover, as we demonstrate in the tutorial, explicit treatment of the gradient allows us to probe systems for multiple degenerate polaronic states.

Important

In what follows, we assume uniform, equal and $\Gamma$-centred $\kk$- and $\qq$-meshes for electron and phonon subspace. However, in principle, it is enough for these meshes to be commensurate.

The polaron formation is modeled in reciprocal space, where the BvK supercell is represented by the correspnding $\kk$-sampling of the Brillouin Zone, i.e. a $2\times2\times2$ $\kk$-mesh corresponds to a $2\times2\times2$ supercell. The polaron binding energy in a supercell of size $N_p$ is given by the variational expression

\[\begin{equation} \label{eq:varpeq} E_\mathrm{pol} \left( \boldsymbol{A}, \boldsymbol{B} \right) = E_\mathrm{el} \left( \boldsymbol{A} \right) + E_\mathrm{ph} \left( \boldsymbol{B} \right) + E_\mathrm{el\text{-}ph} \left( \boldsymbol{A}, \boldsymbol{B} \right), \end{equation}\]

where

\[\begin{equation} \label{eq:varpeq_el} E_\mathrm{el} \left( \boldsymbol{A} \right) = \frac{1}{N_p} \sum_{n\mathbf{k}} | A_{n\mathbf{k}} |^2 \left( \varepsilon_{n\mathbf{k}} - \varepsilon_\mathrm{loc} \right) + \varepsilon_\mathrm{loc}, \end{equation}\]

\[\begin{equation} \label{eq:varpeq_ph} E_\mathrm{ph} \left( \boldsymbol{B} \right) = \frac{1}{N_p} \sum_{\mathbf{q}\nu} | B_{\mathbf{q}\nu} |^2 \omega_{\mathbf{q}\nu}, \end{equation}\]

\[\begin{equation} \label{eq:varpeq_elph} E_\mathrm{el-ph} \left( \boldsymbol{A}, \boldsymbol{B} \right) = -\frac{1}{N_p^2} \sum_{\substack{mn\nu \\ \mathbf{kq}}} A^\ast_{m\mathbf{k+q}} B^\ast_{\mathbf{q}\nu} g_{mn\nu}(\mathbf{k, q}) A_{n\mathbf{k}} + \mathrm{(c.c)} \end{equation}\]

are the electronic, phonon and electron-phonon contribution to the binding energy, respectively.

In Eqs. (\ref{eq:varpeq}-\ref{eq:varpeq_elph}), $\boldsymbol{A} \equiv \left\{ A_{n\mathbf{k}} \right\}$ and $\boldsymbol{B} \equiv \left\{ B_{\mathbf{q}\nu} \right\}$ are the variational coefficients, defining the polaron wavefunction $\phi(\mathbf{r})$ (charge localization) and lattice distortion $\Delta \tau_{\kappa \alpha p}$ as

\[\begin{equation} \label{eq:polaron_phi} \phi (\mathbf{r}) = \frac{1}{\sqrt{N_p}} \sum_{n\mathbf{k}} A_{n\mathbf{k}} \psi_{n\mathbf{k}} (\mathbf{r}), \end{equation}\]

and

\[\begin{equation} \label{eq:polaron_tau} \Delta \tau_{\kappa \alpha p} = - \frac{2}{N_p} \sum_{\mathbf{q}\nu} B^\ast_{\mathbf{q}\nu} \sqrt{ \dfrac{1}{2M_{\kappa}\omega_{\mathbf{q}\nu}} } e_{\kappa \alpha, \nu}(\mathbf{q}) e^{i\mathbf{q}\cdot\mathbf{R}_p}. \end{equation}\]

Charge localization $\varepsilon_\mathrm{loc}$ is represented by the Lagrange multiplier that takes into account the normalization of the polaron wavefunction

\[\begin{equation} \sum_{n\mathbf{k}} |A_{n\mathbf{k}}|^2 = N_p \end{equation}\]

Eqs. (\ref{eq:varpeq}-\ref{eq:polaron_tau}) depend on the input paramters, that can be computed within a unit cell of a system:

electronic bands $\varepsilon_{n\mathbf{k}}$ and Kohn-Sham wavefunctions $\psi_{n\mathbf{k}}$;
phonon dispersion $\omega_{\mathbf{q}\nu}$ and eigenmodes $e_{\kappa \alpha, \nu}(\mathbf{q})$;
electron-phonon matrix elements $g_{mn\nu}(\mathbf{k, q})$.

The variational expression can be optimized self-consistently. At any charge distribution $\boldsymbol{A}$, one can calculate the optimal deformation potential as

\[\begin{equation} \label{eq:phgrad} B_{\mathbf{q}\nu} \left( \boldsymbol{A} \right) = \frac{1}{N_p}\sum_{mn\mathbf{k}} A^\ast_{m\mathbf{k+q}} \frac{g_{mn\nu}(\mathbf{k, q})}{\omega_{\mathbf{q}\nu}} A_{n\mathbf{k}}. \end{equation}\]

Then, the optimal charge localization is given by the Lagrange multiplier

\[\begin{equation} \varepsilon_{\rm loc} \left( \boldsymbol{A}, \boldsymbol{B} \right) = \frac{1}{N_p} \sum_{n\mathbf{k}} | A_{n\mathbf{k}} |^2 \varepsilon_{n\mathbf{k}} - \frac{1}{N_p^2} \sum_{\substack{mn\nu \\ \mathbf{kq}}} \left( A^*_{m\mathbf{k+q}} B^*_{\mathbf{q}\nu} g_{mn\nu}(\mathbf{k}, \mathbf{q}) A_{n\mathbf{k}} + \text{c.c.} \right). \end{equation}\]

The gradient of $E_\mathrm{pol}$ with respect to the electronic degrees of freedom is then obtained as

\[\begin{equation} \begin{split} D_{n\mathbf{k}} \bigl( \boldsymbol{A}, \boldsymbol{B}, \varepsilon_{\rm loc} \bigr) = & \frac{2}{N_p} A_{n\mathbf{k}} \left( \varepsilon_{n\mathbf{k}} - \varepsilon_{\rm loc} \right) \\ & - \frac{2}{N_p^2} \sum_{m\nu\mathbf{q}} A_{m\mathbf{k-q}} B^\ast_{\mathbf{q}\nu} g_{nm\nu}(\mathbf{k-q}, \mathbf{q}) + A_{m\mathbf{k+q}} B_{\mathbf{q}\nu} g^\ast_{mn\nu}(\mathbf{k}, \mathbf{q}). \end{split} \label{eq:elgrad} \end{equation}\]

From this gradient, a new charge distribution $\boldsymbol{A}$ can be determined. Eqs. (\ref{eq:phgrad}-\ref{eq:elgrad}) constitute the building blocks of any gradient-based algorithm for the binding energy optimization. In ABINIT, we employ the preconditioned conjugate gradient algorithm [Payne1992]. After the variational expression is minimized, the optimal $\boldsymbol{A}$ and $\boldsymbol{B}$ coefficients are determined and the polaron in real space can be reconstructed via Eqs. (\ref{eq:polaron_phi}),(\ref{eq:polaron_tau}). Schematically, the optimization process is described in the following chart:

Typical workflow for self-trapped polaron calculations¶

As the description of the VarPEq formalism suggests, prior to the actual polaron calculations, one needs to obtain the required input parameters, describing the electronic and phonon subspaces, and the electron-phonon coupling.

A typical workflow for self-trapped polaron calculations involves the following steps (see the introductory e-ph tutorial):

GS calculation to obtain the WFK and the DEN file. The $\kk$-mesh should be dense enough to converge both electronic and vibrational properties.
DFPT calculations for all the IBZ $\qq$-points corresponding to the ab-initio ddb_ngqpt mesh that will be used to perform the Fourier interpolation of the dynamical matrix and of the DFPT potentials. In the simplest case, one uses a $\qq$-mesh that is equal to the GS $\kk$-mesh (sub-meshes are also fine) and the DFPT calculations can directly start from the WFK produced in step #1. Remember to compute $\bm{\epsilon}^{\infty}$, $\bm{Z}^*$ (polar materials) and the dynamical quadrupoles $\bm{Q}^*$ as these quantities are needed for an accurate interpolation of phonon frequencies and DFPT potentials.
NSCF computation of a WFK file on a much denser $\kk$-mesh representing a BvK supercell hosting a polaron. The NSCF run uses the DEN file produced in step #1. Compute at least one unnocupied electronic state, so later on the code can check if a system has a gap.
Merge the partial DDB and POT files with mrgddb and mrgdv, respectively.
Compute a GSTORE file with eph_task 11. This step performs the interpolation of phonon frequencies and DFPT potentials and computes the electron-phonon matrix elements $g_{mn\nu}(\mathbf{k,q})$ on an arbitrary eph_ngqpt_fine $\mathbf{q}$-mesh for a set of bands gstore_brange or within an energy window gstore_erange near the band edge. The resulting GSTORE.nc is mandatory for VarPEq calculations.

Warning

GSTORE computation requires ABINIT to be configured with support for parallel I/O for HDF5 and NetCDF libraries. See, for example, the corresponding section in the installation tutorial

To verify, whether the installed libraries and ABINIT support the parallel I/O, try issuing

    h5pcc -showconfig | grep "Parallel HDF5"
    nc-config --has-parallel
    abinit -b | grep "Parallel I/O"

If any of these commands gives you a negative prompt, your installation may not work properly.

Start from the WFK, DDB and GSTORE.nc files produced in steps #3-5 to perform VarPEq calculations with eph_task 13.
Visualize polaron wavefuction and displacements using the VPQ.nc file produced in step #6 and eph_task -13.

Getting started¶

Note

Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.

In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:

export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH      # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/             # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Pspdir/  # Do not change this line: path to pseudos dir

Examples in this tutorial use these shell variables: copy and paste the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively, source the set_abienv.sh script located in the ~abinit directory:

source ~abinit/set_abienv.sh

The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.

To execute the tutorials, create a working directory (Work*) and copy there the input files of the lesson.

Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.

Before beginning, you might consider to work in a different subdirectory as for the other tutorials. Why not create Work_eph4vpq in $ABI_TESTS/tutorespfn/Input?

cd $ABI_TESTS/tutorespfn/Input
mkdir Work_eph4vpq
cd Work_eph4vpq

In this tutorial, we prefer to focus on the use of the EPH code hence we will be using precomputed DDB and POT files to bypass the DFPT part. We also provide a DEN.nc file to initialize the NSCF calculations and a GSR.nc file with electronic bands on a $\kk$-path for AbiPy post-processing.

If git is installed on your machine, one can easily fetch the entire repository (29 MB) with:

git clone https://github.com/abinit/LiF_eph_varpeq.git

Alternatively, use wget:

wget https://github.com/abinit/LiF_eph_varpeq/archive/master.zip

or curl:

curl -L https://github.com/abinit/LiF_eph_varpeq/archive/master.zip -o master.zip

or simply copy the tarball by clicking the “download button” available in the github web page, unzip the file and rename the directory with:

unzip master.zip
mv LiF_eph_varpeq-master LiF_eph_varpeq

Warning

The directory with the precomputed files must be located in the same working directory in which you will be executing the tutorial and must be named LiF_eph_varpeq.

The script used to precompute the files is available here, with the necessary *abi files provided in the same repository.

Note that several parameters have been tuned to reach a reasonable compromise between accuracy and computational cost so do not expect the results obtained at the end of the lesson to be fully converged. More specifically, we use norm-conserving pseudopotentials with a cutoff energy ecut of 30 Ha (too low, it should be ~45 Ha). The DFPT computations is done for the set of irreducible $\qq$-points corresponding to a $\Gamma$-centered 3x3x3 $\qq$ mesh (again, too coarse, it should be 6x6x6 for ab initio $\mathbf{k/q}$-meshes). $\bm{Z}^*$ and $\bm{\ee}^\infty$ are also computed with the same underconverged settings.

Important

Note that there is no computation of quadrupoles $\boldsymbol{Q}$, as we use LiF as a prototypical insulator for polaron calculations. In systems with Fm$\overline{3}$m space group, $\boldsymbol{Q} \equiv 0$ by symmetry.

To check if a system has non-zero quadrupole, one can use the following command:

abistruct.py has_quad FILE

where FILE is any file providing a structure in ABINIT format (e.g. netcdf output files or .abi files)

How to extract useful info from the output files¶

For post-processing, this tutorial will rely heavily on the AbiPy toolkit. The installation notes and some examples can be found in the AbiPy repository.

As an example, to plot the band structure using the GS eigenvalues stored in the GSR.nc file of the tutorial’s repo, use the abiopen.py script provided by AbiPy with the -e option:

abiopen.py LiF_eph_varpeq/ebandso_GSR.nc -e

Note the VBM is triply degenerate if SOC is not used. This suggests that we will have to use last three occupied bands for hole polaron calculations. As the CBM has no degeneracy, only the first unoccupied band will be considered in electron polaron calculations.

How to merge partial DDB files with mrgddb¶

Normally, after DFPT calculations, we need to merge partial DDB files. This has already been done in the repository at the end of the run script via the

mrgddb out_DDB dfpt*DDB.nc

This command takes advantage of globbing and merges all the partial DDB files into a single out_DDB file, containing information on all the calculated perturbations.

Since we are dealing with a polar material, it is worth checking whether our final DDB contains Born effective charges and the electronic dielectric tensor. Instead of running anaddb or abinit and then checking the output file, we can simply use abiopen.py with the -p option:

abiopen.py LiF_eph_varpeq/out_DDB -p

================================== DDB Info ==================================

Number of q-points in DDB: 4
guessed_ngqpt: [3 3 3] (guess for the q-mesh divisions made by AbiPy)
ecut = 30.000000, ecutsm = 0.000000, nkpt = 14, nsym = 48, usepaw = 0
nsppol 1, nspinor 1, nspden 1, ixc = 11, occopt = 1, tsmear = 0.010000

Has total energy: False
Has forces: False
Has stress tensor: False

Has (at least one) atomic pertubation: True
Has (at least one diagonal) electric-field perturbation: True
Has (at least one) Born effective charge: True
Has (all) strain terms: False
Has (all) internal strain terms: False
Has (all) piezoelectric terms: False
Has (all) dynamical quadrupole terms: False

We can also invoke anaddb directly from python to have a quick look at the phonon dispersion:

abiview.py ddb LiF_eph_varpeq/out_DDB

Computing phonon bands and DOS from DDB file with:
nqsmall: 10, ndivsm: 20,
asr: 2, chneut: 1, dipdip: 1, lo_to_splitting: automatic, dos_method: tetra

that produces the following figures:

The results seem reasonable: the acoustic modes go to zero linearly for $\qq \rightarrow 0$ as we are dealing with a 3D system, no instability is present and the phonon dispersion shows the LO-TO splitting typical of polar materials.

Note, however, that the acoustic sum-rule is automatically enforced by the code so it is always a good idea to compare the results with/without asr as this is an indirect indicator of the convergence/reliability of our calculations. We can automate the process with the ddb_asr command of abiview.py :

abiview.py ddb_asr LiF_eph_varpeq/out_DDB

that produces the following figure:

Important

In this case, the small breaking of the acoustic sum rule is mainly due the too low cutoff energy employed in our calculations. In production, one should perform convergence tests wrt ecut and $\qq$-mesh but since the goal of this lesson is to teach you how to run polaron calculations, we ignore this problem and continue with the other examples.

PS: If you want to compute the phonons bands with/without dipdip, use:

abiview.py ddb_dipdip LiF_eph_varpeq/out_DDB

How to merge partial POT files with mrgdv¶

Now we can merge the DFPT potential with the mrgdv tool using the command.

mrgdv < teph4vpq_1.abi > teph4vpq_1.stdout

and the following input file:

This can also be achieved by calling mrgdv with the merge argument, output filename and wildcard for the POT files

mrgdv merge teph4vpq_1_DVDB LiF_eph_varpeq/dfpt*POT*nc

Tip

The number at the end of the POT file corresponds to the (idir, ipert) perturbation for that particular $\qq$-point. The pertcase index is computed as:

pertcase = idir + 3 * (ipert-1)

where idir gives the direction ([1, 2, 3]) and ipert specifies the perturbation type:

ipert in [1, …, natom] corresponds to atomic perturbations (reduced directions)
ipert = natom + 1 corresponds d/dk (reduced directions)
ipert = natom + 2 corresponds the electric field
ipert = natom + 3 corresponds the uniaxial stress (Cartesian directions)
ipert = natom + 4 corresponds the shear stress. (Cartesian directions)

All DFPT POT files with 1 <= pertcase <= 3 x natom therefore correspond to atomic pertubations for a given $\qq$-point.

The value of pertcase and qpt are reported in the ABINIT header. To print the header to terminal, use abitk with the hdr_print command

 abitk hdr_print LiF_eph_varpeq/dfpto_DS2_POT4.nc

 ===============================================================================
 ECHO of part of the ABINIT file header

 First record :
.codvsn,headform,fform = 10.3.4.6   80  111

 Second record :
 bantot,intxc,ixc,natom  =    40     0    11     2
 ngfft(1:3),nkpt         =    30    30    30     8
 nspden,nspinor          =     1     1
 nsppol,nsym,npsp,ntypat =     1    48     2     2
 occopt,pertcase,usepaw  =     1     4     0
 ecut,ecutdg,ecutsm      =  3.0000000000E+01  3.0000000000E+01  0.0000000000E+00
 ecut_eff                =  3.0000000000E+01
 qptn(1:3)               =  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
 rprimd(1:3,1)           =  0.0000000000E+00  3.8375000000E+00  3.8375000000E+00
 rprimd(1:3,2)           =  3.8375000000E+00  0.0000000000E+00  3.8375000000E+00
 rprimd(1:3,3)           =  3.8375000000E+00  3.8375000000E+00  0.0000000000E+00
 stmbias,tphysel,tsmear  =  0.0000000000E+00  0.0000000000E+00  1.0000000000E-02

 The header contain   4 additional records.

Use --prtvol 1 to output more records.

Now we discuss in more detail the output file produced by mrgdv

For each $\qq$-point found in the partial POT files, the code prints a lists with the atomic perturbations that have been merged in the database.

 qpoint: [ 0.0000E+00,  0.0000E+00,  0.0000E+00] is present in the DVDB file
 The list of irreducible perturbations for this q vector is:
    1)  idir= 1, ipert=   1, type=independent, found=Yes
    2)  idir= 2, ipert=   1, type=symmetric, found=No
    3)  idir= 3, ipert=   1, type=symmetric, found=No
    4)  idir= 1, ipert=   2, type=independent, found=Yes
    5)  idir= 2, ipert=   2, type=symmetric, found=No
    6)  idir= 3, ipert=   2, type=symmetric, found=No

The term symmetric means that this particular (idir, ipert) perturbation can be reconstructed by symmetry from the other independent entries with the same $\qq$-point. If all the independent entries are available, the code prints the following message at the end of the output file:

 All the independent perturbations are available
 Done

Warning

If you don’t get this message, the DVDB cannot be used by the EPH code. In this case, check carefully your DFPT input files and the list of POT files that have been merged.

Also, note that it is not possible to change the value of nsym at the level of the EPH calculation as symmetries are automatically inherited from the previous GS/DFPT calculations. To perform calculations without symmetries, one should set nsym to 1 (no spatial symmetries) and use kptopt 4 to disable time-reversal symmetry in all the input files.

Computing the WFK files¶

At this point we have all the ingredients (DDB and DVDB) required to compute/interpolate the e-ph scattering potentials, and we can finally start to generate the WFK files.

For our first NSCF calculation, we compute wavefunctions on $\Gamma$-centered $\kk$-meshes (3x3x3, 4x4x4, 5x5x5, and 13x13x13) with 6 bands (5 occupied + 1 unnocupied). The first three wavefunctions will be used in calculations of small hole polaron in LiF. The latter will be used for the large electron polaron in LiF. Note the use of getden_filepath to read the precomputed DEN.nc file instead of getden or irdden.

You may now run the NSCF calculation by issuing:

abinit teph4vpq_2.abi > teph4vpq_2.log 2> err &

with the input file given by:

Small hole polaron: computing the GSTORE files¶

For our first GSTORE calculation, we compute the electron-phonon matrix elements required for the small hole polaron optimization. With this example, we discuss the most important gstore_ input variables and the concept of the GSTORE file. First of all, you may want to start immediately the computation by issuing:

abinit teph4vpq_3.abi > teph4vpq_3.log 2> err &

with the following input file:

Tip

To run the examples in parallel with e.g 2 MPI processes use:

mpirun -n 2 abinit teph4vpq_3.abi > teph4vpq_3.log 2> err &

The EPH code will automatically distribute the workload using a predefined distribution scheme (not necessarily the most efficient in terms of memory and wall-time). The input variable eph_np_pqbks can be used to specify manually the distribution scheme.

Warning

GSTORE calculations imply various level of the MPI parallelism and parallel IO. At the time of writing this tutorial, the bottleneck for number of MPI processes is the number of $\mathbf{k}/\mathbf{q}$-points in the IBZ. Specifing more may lead to idle processes that are handled differently by different compileres.

As we use ndtset mode to compute three GSTORE files in a single run, the first dataset (3x3x3 $\kk/\qq$-mesh) will be the limiting factor, and here the maximum number of allowed MPI processes is 4. In production, it is highly recommended to compute seperate GSTORE files in individual runs.

Let’s now discuss the meaning of the different variables in more detail. We use optdriver 7 to enter the EPH code while eph_task 11 activates the computation of the electron-phonon matrix elements via interpolation of the DFPT scattering potentials.

The paths to the external files (DDB, WFK, DVDB) are specified with the three variables:

This is an excerpt of the input file:

# Electron-phonon matrix elements parameters
 gstore_cplex  2     # Complex e-ph matrix elements   (mandatory for VarPEq)
 gstore_kzone  "bz"  # Use full BZ for k and q meshes (mandatory for VarPEq)
 gstore_qzone  "bz"
 gstore_brange 3 5   # Bands participating in the hole polaron formation:
                     # from 3 to 5 (triply degenerate valence band maximum)
# Wavefunction
# Read the WFK file from the preparatory run
 getwfk_filepath1 "teph4vpq_2o_DS1_WFK"  # 3x3x3 k-mesh
 getwfk_filepath2 "teph4vpq_2o_DS2_WFK"  # 4x4x4 k-mesh
 getwfk_filepath3 "teph4vpq_2o_DS3_WFK"  # 5x5x5 k-mesh

 ngkpt1   3 3 3  # These parameters have to be consistent with the WFK file
 ngkpt2   4 4 4
 ngkpt3   5 5 5
 ...

# Phonons
 getddb_filepath
 "LiF_eph_varpeq/out_DDB"  # Read the precomputed DDB stored in git submodule
                           # NOTE: for the input parameters used to produce
                           # the DDB file, consult the LiF_eph_varpeq/dfpt.abi

# Scattering potentials
 getdvdb_filepath
 "teph4vpq_1_DVDB"  # Read the DVDB file from the preparatory run

# Phonons and scattering potentials interpolation
 ddb_ngqpt  3 3 3       # This has to be consistent with the DDB/DVDB files
 ...

 eph_ngqpt_fine1  3 3 3  # Interpolate to 3x3x3 q-grid
 eph_ngqpt_fine2  4 4 4  #                4x4x4 q-grid
 eph_ngqpt_fine3  5 5 5  #                5x5x5 q-grid

The mesh for electrons (ngkpt, nshiftk and shiftk) must correspond to the one used for the input WFK file. ddb_ngqpt is set to 3x3x3 as this is the $\qq$-mesh used in the DFPT part to generate the DDB and DVDB files, but the integration in $\qq$-space is performed with the eph_ngqpt_fine mesh. As eph_ngqpt_fine differs from ddb_ngqpt, the code will automatically activate the interpolation of the DFPT potentials as discussed in the introduction to the EPH code.

Polaron calculations require complex electron-phonon matrix elements, hence gstore_cplex 2 is used. As we will deal with small hole polaron, and the VBM in LiF is triply-degenerate, we select the last three occupied bands via gstore_brange variable.

Important

At present, at the level of polaron computations, no symmetry is assumed for the electron-phonon matrix elements. Hence these terms have to be computed in the full BZ for $\kk$ and $\qq$ vectors. This behaviour is specified by the gstore_kzone and gstore_qzone variables. As a result, for some systems GSTORE files may require a relatively large amount of phyical memory.

In principle, it is possible to utilize symmetries for matrix elements, albeit this requires the correct treatment of the scattering between degenerate states via degenerate phonons [Li2024]. This feature, however, is not yet implemented.

We can now have a look at the main output file:

For each DATASET we find a section that summarizes the gstore parameters:

 === Gstore parameters ===
 kzone: bz
 kfilter: none
 nkibz: 4
 nkbz: 27
 glob_nk_spin: [27]
 qzone: bz
 nqibz: 4
 nqbz: 27
 glob_nq_spin: [27]
 kptopt: 1
 qptopt: 1
 with_vk: 1

 gqk_cplex: 2
 gqk_bstart: 3
 gqk_bstop: 5
 gqk_nb: 3
 gqk_my_npert: 6
 gqk_my_nk: -1
 gqk_my_nq: -1

Then we find another section related to the MPI parallelism. In this case we are running in sequential but the output will change if we run in parallel (see also eph_np_pqbks). The final message informs the user that the EPH code will either read the qpts from file (if the DVDB contains all of them, in case eph_ngqpt_fine is not defined in the input) or interpolate the scattering potentials from ddb_ngqpt to eph_ngqpt_fine.

  === MPI distribution ===
P Number of CPUs for parallelism over perturbations: 1
P Number of perturbations treated by this CPU: 6
P Number of CPUs for parallelism over q-points: 1
P Number of CPUs for parallelism over k-points: 1

 DVDB file contains all q-points in the IBZ --> Reading DFPT potentials from file.

or

 Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.

Finally, we have the section with the squared absolute values of specific matrix elements are printed. This is mostly done for regression tests of ABINIT. Note, however, the gstore_completed: 1 tag. This indicates that all the matrix elments have been computed for the corresponding GSTORE file.

 gstore_completed: 1
 gstore_done_qbz_spin: 125
 gqk%nb: 3
 gqk%glob_nq: 125
 gqk%glob_nk: 125
 E-PH matrix elements:
    iq    ik  mode im_kq  in_k    |g|^2 in Ha^2
 |g(k,q)|^2 in Ha^2 for iq: 1 ik: 1 mode: 4
    iq    ik  mode im_kq  in_k    |g|^2 in Ha^2
     1     1     4     1     1     3.315387E-18
     1     1     4     1     2     5.564496E-19
     1     1     4     1     3     6.911046E-19
     1     1     4     2     1     5.564496E-19
     1     1     4     2     2     3.620046E-18
     1     1     4     2     3     1.066932E-19
     1     1     4     3     1     6.911046E-19
     1     1     4     3     2     1.066932E-19
     1     1     4     3     3     3.147723E-18
...

The calculation has produced the following output files (e.g. for the first DATASET):

$ ls teph4vpq_3o_DS1*

teph4vpq_3o_DS1_EBANDS.agr teph4vpq_3o_DS1_GSTORE.nc  teph4vpq_3o_DS1_PHDOS.nc

where:

EBANDS.agr → Electron bands in xmgrace format. See also prtebands
PHDOS.nc → Phonon DOS in netcdf format (see prtphdos is given by ph_ngqpt).
GSTORE.nc → Netcdf file with electron-phonon matrix elements $g_{mn\nu}(\mathbf{k,q})$ and other related datastructures.

All the matrix elements are stored in the GSTORE.nc NetCDF file for all $\kk/\qq$-points, electronic states and phonon modes. As usual, one can use abiopen.py with the -p option (--print) to print a summary to terminal:

$ abiopen.py teph4vpq_3o_DS1_GSTORE.nc -p
...
nsppol: 1
gstore_completed: True
gstore_cplex: 2
gstore_kzone: bz
gstore_kfilter: none
gstore_gmode: phonon
gstore_qzone: bz
gstore_with_vk: 1
gstore_kptopt: 1
gstore_qptopt: 1
gstore_brange_spin[0]: [2 5]
gstore_erange_spin[0]: [0. 0.]
gstore_glob_spin_nq[0]: 27

Important

Pay attention to the gstore_completed tag. If set to False, some matrix elements may not have been computed. This typically occurs if the job was interrupted, for example due to a timeout.

To restart a gstore computation, simply instruct the code to read the incomplete GSTORE.nc file using getgstore_filepath, and rerun the same input file with optdriver = 11. If needed, you can disable this automatic restart capability by setting eph_restart to 0.

Small hole polaron: variational optimization of binding energy¶

At this point we have all the ingredients to start our first VarPEq calcualation. This will be done for three $\kk/\qq$-meshes of incrementally increasing sampling: 3x3x3, 4x4x4, and 5x5x5, which represent a supercell hosting the polaron. For each mesh, we have already computed the mandatory inputs: WFK and GSTORE files. We will also need the same DDB file used in the GSTORE calculations.

With this example, we can discuss the most important vpq_ input variables and the concept of VarPEq calculations. First of all, you may want to start the computation immediately by issuing:

abinit teph4vpq_4.abi > teph4vpq_4.log 2> err &

with the following input file:

Tip

In this example, the GSTORE files have already been computed and will be reloaded from disk. The EPH code will automatically distribute the workload and memory using a predifined distribution scheme.

The manual parallelism is possible over $\kk$-, $\qq$-points, pertubations and spin indices. See the related eph_np_pqbks input variable.

Important

We use ndtset mode to perform three polaron calculations in a single run. Hence, the first dataset (3x3x3 $\kk/\qq$-mesh) will be limiting factor for parallelism: max 6 procs over pertubations (2*natom), and max 9 procs for $\kk/\qq$-vectors each. In production, it is highly recommended to compute seperate GSTORE files in individual runs.

Let’s now discuss the meaning of the different variables in more detail. As usual, we use optdriver 7 to enter the EPH code while eph_task 13 activates the solver for the variational polaron equations.

The paths to the mandatory external files (DDB, WFK, GSTORE) are specified by the three variables:

This is an excerpt of the input file:

# Variational Polaron Equations parameters
 vpq_pkind  "hole"  # Hole polaron
 vpq_tolgrs  1e-5   # Tolerance on the gradient norm

 # these parameters are default but given here explicitly for clarity
 vpq_aseed      "gau_energy"  # Initial charge localization is Guassian-like
 vpq_gpr_energy  0. 1.        # Gaussian parameters for initial charge
 vpq_nstates     1            # Number of polaronic states
 vpq_nstep       50           # Maximum number of optimization steps

# Frohlich correction
 eph_frohl_ntheta  32   # Number of divisions for the spherical integration of
                        # the long-range part of electron-phonon energy.
                        # Non-zero value activates the computation of Frohlich
                        # corrections.
                        # ntheta=32 gives the optimal number of divisions for the
                        # spherical integration.

 vpq_avg_g  1  # Apply Frohlich correction to e-ph coupling


# Electron-phonon matrix elements
 getgstore_filepath1 "teph4vpq_3o_DS1_GSTORE.nc" # 3x3x3 k/q-mesh
 getgstore_filepath2 "teph4vpq_3o_DS2_GSTORE.nc" # 4x4x4 k/q-mesh
 getgstore_filepath3 "teph4vpq_3o_DS3_GSTORE.nc" # 5x5x5 k/q-mesh

# Wavefunction
# Read the WFK file from the preparatory run
 getwfk_filepath1 "teph4vpq_2o_DS1_WFK"  # 3x3x3 k-mesh
 getwfk_filepath2 "teph4vpq_2o_DS2_WFK"  # 4x4x4 k-mesh
 getwfk_filepath3 "teph4vpq_2o_DS3_WFK"  # 5x5x5 k-mesh

 ngkpt1   3 3 3  # These parameters have to consistent with the WFK file
 ngkpt2   4 4 4
 ngkpt3   5 5 5
 ...

# Phonons
 ddb_ngqpt  3 3 3          # This has to be consistent with the DDB file
 getddb_filepath
 "LiF_eph_varpeq/out_DDB"  # Read the precomputed DDB stored in git submodule
                           # NOTE: for the input parameters used to produce
                           # the DDB file, consult the LiF_eph_varpeq/dfpt.abi

The mesh for electrons (ngkpt, nshiftk and shiftk) must correspond to the one used for the input WFK file. ddb_ngqpt is set to 3x3x3 as this is the $\qq$-mesh used in the DFPT part to generate the DDB and DVDB files used previously in GSTORE.

Polaron calculations require an explicit vpq_pkind variable specification which tells the code which polaron we are dealing with (here, "hole"). In principle, this is the only mandatory vpq_ variable, and other inputs take default values if not specified.

In our example, the initial guess for the electronic coefficients $\boldsymbol{A}$ is given by the combination of vpq_aseed and vpq_gpr_energy variables. The vpq_nstates variable gives the number of polaronic states. For each state, the solver will iterate until it converges on the gradient norm specified by the vpq_tolgrs variable OR reaches the maximum number of iterations vpq_nstep.

Here, we also use a Fr"ohlich correction for the electron-phonon matrix divergence at $\Gamma$. For a fixed $\kk/\qq$-mesh, this trick brings the solution closer to the infinite-mesh limit, and also facilitates the convergence of our algorithm. One activates this correction via vpq_avg_g tag and the specification of the eph_frohl_ntheta variable.

We can now have a look at the main output file:

Since first we load the GSTORE files, we find a section that summarizes the gstore parametes and MPI parallelism:

 Computing variational polaron equations from pre-existent GSTORE file: teph4vpq\_3o\_DS1_GSTORE.nc
 Initializing gstore object from: teph4vpq_3o_DS1_GSTORE.nc

 === Gstore parameters ===
  ...

  === MPI distribution ===
  ...

The following section promts the parameters of the VarPEq optimization and other related quantities:

  === Variational Polaron Equations ===

 Polaron:
   * Polaron kind: hole
   * BvK supercell: 3x3x3
   * Number of independent spin polarizations: 1
   * Number of polaronic states: 1
   * Filtering of electronic states: no

 Long-range corrections:
   * Frohlich correction: yes
   * Frohlich correction value: 9.265067E-01 eV
   * Frohlich correction included in matrix elements: yes

 Optimization parameters:
   * Initial seed: charge localization A_nk
   * Initial seed type: Gaussian, based on the electroic energies
   * Tolerance on the gradient norm: 1.00E-05
   * Maximum number of iterations per state: 50

And, finally, the code enters the solver and prints the optimization log at the end of computation.

 Solving the variational polaron equations for each state...
   * spin 1/1, pstate 1/1...
   Done

 Printing the optimization logs
   --------------------------------------------------------------------------------------
   * spin 1/1, pstate 1/1
   * values in the optimization log are in (a.u.)
   --------------------------------------------------------------------------------------
   Step        E_pol         E_el         E_ph       E_elph      epsilon     ||grad||
      1  -3.0118E-03   4.1093E-02   4.4105E-02  -8.8210E-02   4.7117E-02   1.0900E-02
      2  -2.4528E-02   1.9417E-02   4.3945E-02  -8.7891E-02   6.8473E-02   2.2198E-03
      3  -3.9895E-02   3.0152E-03   4.2910E-02  -8.5819E-02   8.2804E-02   4.9899E-03
      4  -4.1219E-02   1.7941E-03   4.3013E-02  -8.6026E-02   8.4232E-02   3.2342E-03
      5  -4.3429E-02   5.7133E-04   4.4000E-02  -8.8000E-02   8.7429E-02   1.8175E-03
      6  -4.7595E-02   2.6635E-03   5.0259E-02  -1.0052E-01   9.7854E-02   4.6156E-03
      7  -5.2013E-02   5.3655E-03   5.7378E-02  -1.1476E-01   1.0939E-01   6.3385E-03
      8  -5.6356E-02   8.0023E-03   6.4359E-02  -1.2872E-01   1.2072E-01   6.1036E-03
      9  -6.1355E-02   1.0870E-02   7.2225E-02  -1.4445E-01   1.3358E-01   2.7021E-03
     10  -6.3184E-02   1.1570E-02   7.4754E-02  -1.4951E-01   1.3794E-01   4.1324E-03
     11  -6.9063E-02   1.5656E-02   8.4719E-02  -1.6944E-01   1.5378E-01   6.5822E-03
     12  -7.5399E-02   2.2493E-02   9.7892E-02  -1.9578E-01   1.7329E-01   3.7146E-03
     13  -7.6048E-02   2.4094E-02   1.0014E-01  -2.0028E-01   1.7619E-01   2.3119E-03
     14  -7.6484E-02   2.6017E-02   1.0250E-01  -2.0500E-01   1.7899E-01   3.8701E-04
     15  -7.6489E-02   2.6203E-02   1.0269E-01  -2.0538E-01   1.7918E-01   2.9740E-04
     16  -7.6495E-02   2.6384E-02   1.0288E-01  -2.0576E-01   1.7937E-01   5.8608E-05
     17  -7.6495E-02   2.6413E-02   1.0291E-01  -2.0582E-01   1.7940E-01   4.1064E-05
     18  -7.6496E-02   2.6449E-02   1.0294E-01  -2.0589E-01   1.7944E-01   4.5160E-05
     19  -7.6499E-02   2.6426E-02   1.0293E-01  -2.0585E-01   1.7942E-01   6.9919E-05
     20  -7.6504E-02   2.6408E-02   1.0291E-01  -2.0582E-01   1.7942E-01   1.6567E-04
     21  -7.6524E-02   2.6393E-02   1.0292E-01  -2.0583E-01   1.7944E-01   2.8985E-04
     22  -7.6587E-02   2.6434E-02   1.0302E-01  -2.0604E-01   1.7961E-01   3.2632E-04
     23  -7.6771E-02   2.6522E-02   1.0329E-01  -2.0658E-01   1.8006E-01   5.2643E-04
     24  -7.7070E-02   2.6540E-02   1.0361E-01  -2.0722E-01   1.8068E-01   6.7279E-04
     25  -7.7399E-02   2.6568E-02   1.0397E-01  -2.0793E-01   1.8137E-01   5.2144E-04
     26  -7.7402E-02   2.6586E-02   1.0399E-01  -2.0798E-01   1.8139E-01   3.8771E-04
     27  -7.7485E-02   2.6715E-02   1.0420E-01  -2.0840E-01   1.8168E-01   7.3990E-04
     28  -7.7588E-02   2.6682E-02   1.0427E-01  -2.0854E-01   1.8186E-01   2.9923E-04
     29  -7.7616E-02   2.6588E-02   1.0420E-01  -2.0841E-01   1.8182E-01   2.7317E-04
     30  -7.7633E-02   2.6620E-02   1.0425E-01  -2.0850E-01   1.8189E-01   1.2044E-04
     31  -7.7633E-02   2.6648E-02   1.0428E-01  -2.0856E-01   1.8191E-01   5.7835E-05
     32  -7.7636E-02   2.6711E-02   1.0435E-01  -2.0869E-01   1.8198E-01   9.0984E-05
     33  -7.7637E-02   2.6646E-02   1.0428E-01  -2.0857E-01   1.8192E-01   5.3464E-05
     34  -7.7637E-02   2.6617E-02   1.0425E-01  -2.0851E-01   1.8189E-01   4.7310E-05
     35  -7.7638E-02   2.6658E-02   1.0430E-01  -2.0859E-01   1.8193E-01   1.9031E-05
     36  -7.7638E-02   2.6661E-02   1.0430E-01  -2.0860E-01   1.8194E-01   1.8077E-05
     37  -7.7638E-02   2.6652E-02   1.0429E-01  -2.0858E-01   1.8193E-01   1.9269E-05
     38  -7.7638E-02   2.6654E-02   1.0429E-01  -2.0858E-01   1.8193E-01   7.3431E-06
   --------------------------------------------------------------------------------------
   Converged: ||grad||=7.3431E-06 < vpq_tolgrs=1.0000E-05
   E_pol (eV):  -2.11263475E+00
     eps (eV):   4.95056023E+00
   --------------------------------------------------------------------------------------

 Saving results to: teph4vpq_4o_DS1_VPQ.nc

In the log above, one can see the convergence of the polaron binding energy E_pol, its individual electronic, phonon and electron-phonon terms E_el, E_ph, and E_elph, as well as charge localization energy epsilon. The ||grad|| column shows the evolution of the gradient norm over the iterations.

Tip

Note, however, that this log is given to get a first glimpse on the behaviour of the solver. Complete results, including optimized $\boldsymbol{A}$ and $\boldsymbol{B}$ coefficients are saved in the VPQ NetCDF at the end of computation. We will heavily rely on these netcdf files for post-processing.

The calculation has produced the following output files (e.g. for the first DATASET):

$ ls teph4vpq_4o_DS1*

teph4vpq_4o_DS1_EBANDS.agr teph4vpq_4o_DS1_PHDOS.nc   teph4vpq_4o_DS1_VPQ.nc

where the most important is VPQ.nc – a NetCDF file with the optimization log, optimized $\boldsymbol{A}$ and $\boldsymbol{B}$ coefficients and other relevant data. We will heavily rely on these files and the AbiPy toolkit for post-processing our results. As usual, one can use abiopen.py with the -e option (--expose) to print a summary to terminal and visualize the optimization log:

The subplot on the left shows the evolution of the gradient norm. Two subplots on the right show the convergence of energy terms and their log-scale difference wrt final values. As the binding energy $E_\mathrm{pol}$ is the only variational quantity, on these two subplots it is shown by solid lines.

Small hole polaron: extrapolation to the infinite-size supercell¶

As seen from our calculations, the small hole polaron is already being localized at 3x3x3 $\kk/\qq$-mesh (3x3x3 supercell). However, to evalute its binding energy in the infinite-size limit $N_p \to \infty$, we need to perform the extrapolation wrt inverse supercell size (that’s why we have optimized polaron for multiple meshes!):

\[\begin{equation} \label{eq:makovpayne} E_\mathrm{pol}(N_p) = E^\infty_\mathrm{pol} + a N_p^{-1/3} + \mathcal{O}(N_p^{-1}), \end{equation}\]

To fit the parameters of Eq. (\ref{eq:makovpayne}) to the optimized data, one can use any suitable software or library of choice. In this tutorial, however – just as in other examples – we will rely on the capabilities of AbiPy as a Python module.

AbiPy provides a Robot class that can be used to analzye multiple NetCDF files produced in ABINIT calculations. In particular, we will be interested in the VpqRobot. The following Python script gives an example on how to use this class with the relevant methods:

from abipy.eph.vpq import VpqRobot

vpq_files = [
    "teph4vpq_4o_DS1_VPQ.nc",
    "teph4vpq_4o_DS2_VPQ.nc",
    "teph4vpq_4o_DS3_VPQ.nc",
]
robot = VpqRobot.from_files(vpq_files)

df = robot.get_final_results_df()
print(df)

robot.plot_kconv(nfit=3)

It imports the VpqRobot class from the module and initializes it with the list of VPQ files. Then, a pandas dataframe with optimized data for each VPQ file is returned via the get_final_results_df method. One can use this dataframe to perform further data analysis and create custom publication-ready figures. We then use a plot_kconv method that peforms the fitting of Eq. (\ref{eq:makovpayne}) parameters to the data in the robot, and plots the following figure:

Tip

For formal arguments of the methods and their meaning, user can invoke the Python help() funtion. E.g. calling help(VpqRobot.plot_kconv) gives the method’s docstring:

...
Plot the convergence of the results with respect to k-point sampling.

Args:
    nfit (int, optional): Number of points used in linear extrapolation. Defaults to 3.
...

The extrapolated values $E_\mathrm{pol}$ = -2.03 eV and $\varepsilon_\mathrm{loc}$ = 4.99 eV differ from a fully converged values of $E_\mathrm{pol}$ = -1.94 eV and $\varepsilon_\mathrm{loc}$ = 4.68 eV. This is mainly due to the low ecut and sampling of $\kk/\qq$-meshes used to produce the input parameters for the solver. Also, one has to impose lower vpq_tolgrs convergence threshold and use denser meshes for extrapolation.

Now, let’s try to perform the same calculations as in teph4vpq_4.abi but without the long-range corrections for the matrix elements at $\Gamma$. The input is the same, but with the vpq_avg_g tag set to 0:

As usual, run this example with:

abinit teph4vpq_5.abi > teph4vpq_5.log 2> err &

Upon completion, we can modify our Python script to include the new VPQ files in our robot.

from abipy.eph.vpq import VpqRobot

vpq_files = [
    # calculations with LR correction
    "teph4vpq_4o_DS1_VPQ.nc",
    "teph4vpq_4o_DS2_VPQ.nc",
    "teph4vpq_4o_DS3_VPQ.nc",
    # calculations without LR correction
    "teph4vpq_5o_DS1_VPQ.nc",
    "teph4vpq_5o_DS2_VPQ.nc",
    "teph4vpq_5o_DS3_VPQ.nc",
]
robot = VpqRobot.from_files(vpq_files)

robot.plot_kconv(nfit=3)

The robot is smart enough (but stupid nonetheless, so trust your own judgemenet) to detect 2 different kind of calculations and peform separate extrapolation for each case. It produces the following figure:

In principle, extrapolated values for $E_\mathrm{pol}$ and $\varepsilon_\mathrm{loc}$ have to agree whether or not the LR correction is employed. The discrepancy here is mainly due to the relatively high value of vpq_tolgrs and the coarse $\kk/\qq$-meshes used for extrapolation.

Exercise

Try changing the nfit argument in plot_kconv to 2 (i.e. throw away the 3x3x3 point) from the extrapolation data and see if the extrapolated values agree better.

Note, however, how adding the LR correction brings the values of $E_\mathrm{pol}$ at fixed $\kk/\qq$-meshes closer to the infinite-size value.

Warning

For delocalized polarons, i.e. there is no binding, adding LR correction will anyway give the negative value of $E_\mathrm{pol}$, solely due to the long-range correction that vanishes in the limit of $N_p \to \infty$. Without the LR correction, $E_\mathrm{pol}$ is zero for a delocalized particle. As a result, one alway needs to check the extrapolation to the infinite-size limit and reach $\kk/\qq$-meshes (supercells) large enough to host a polaron.

Small hole polaron: multiple states and visualiziation¶

For now, we have probed our system for a single polaronic state. However, our method supports the search for multiple polaronic states: once solutions with charge localization $\boldsymbol{A}^i, i < n$ are found, we find the next one $\boldsymbol{A}^{n}$ by imposing an orthogonality constraint wrt $\boldsymbol{A}^i, i < n$.

This is similar to the approach used to compute multiple electornic bands [Payne1992]. In polaron problem, however, this procedure does not give an optimized solution for $\boldsymbol{A}^n$, as polaronic states are not necessarily orthogonal to each other. Rather, it brings the next solution far away from all the previous ones. After this is achieved, i.e. the solver is stuck at a constrained-induced stationary point, the constraint can be removed, and a new solution is found.

Here we will utilize this approach to demonstrate that there exist multiple variational solutions for hole polaron in LiF.

Our input file is

As usual, run this example with:

abinit teph4vpq_6.abi > teph4vpq_6.log 2> err &

Note that vpq_nstates variable is now 3. Also, two new vpq_ variables are introduced, vpq_nstep_ort and vpq_translate. The former defines how many orthogonalization steps are to be performed. After vpq_nstep_ort steps, the orthogonality constraint is lifted, and the algorithm relaxes to a solution from a starting point for $\boldsymbol{A}$ that is far away from the previous ones. The vpq_translate tag turns on the translation of polaronic solutions, i.e. we peform $n$-th orthogonalization wrt computed $\boldsymbol{A}^i, i < n$ and their translated images within the supercell. This is rather important point, as polarons can localize at any equivalent atomic site within a supercell, and we have to eliminate these trivial solutions.

Important

At large $\kk/\qq$ meshes the translation of polaronic solutions can be computationally prohibitive due to the large supercell size.

Now let’s examine the main output file

Note the difference in the initial prompt

...
 Solving the variational polaron equations for each state...
   * spin 1/1, pstate 1/3...
   Done
   * spin 1/1, pstate 2/3...
   Done
   * spin 1/1, pstate 3/3...

Also, there are now three optimization logs for each state, which are a bit different when the orthogonalization is peformed. For example, for state 2:

   --------------------------------------------------------------------------------------
   * spin 1/1, pstate 2/3
   * values in the optimization log are in (a.u.)
   (o) - orthogonal all pstates < 2
   --------------------------------------------------------------------------------------
   Step        E_pol         E_el         E_ph       E_elph      epsilon     ||grad||
      1  -9.2204E-04   4.4614E-02   4.5536E-02  -9.1071E-02   4.6458E-02   1.0902E-02  (o)
      2  -2.8019E-02   1.8675E-02   4.6694E-02  -9.3388E-02   7.4713E-02   5.0189E-03  (o)
      3  -4.0570E-02   7.1368E-03   4.7707E-02  -9.5414E-02   8.8277E-02   3.2275E-03  (o)
      4  -4.4896E-02   4.8512E-03   4.9747E-02  -9.9494E-02   9.4642E-02   4.2843E-03  (o)
      5  -5.1065E-02   5.2935E-03   5.6358E-02  -1.1272E-01   1.0742E-01   6.1297E-03  (o)
      6  -5.7578E-02   7.8267E-03   6.5405E-02  -1.3081E-01   1.2298E-01   5.4966E-03  (o)
      7  -6.1360E-02   1.0590E-02   7.1950E-02  -1.4390E-01   1.3331E-01   2.7006E-03  (o)
      8  -6.1924E-02   1.1454E-02   7.3377E-02  -1.4675E-01   1.3530E-01   3.4064E-03  (o)
      9  -6.3122E-02   1.3396E-02   7.6519E-02  -1.5304E-01   1.3964E-01   5.3188E-03  (o)
     10  -6.5647E-02   1.7329E-02   8.2976E-02  -1.6595E-01   1.4862E-01   7.5766E-03  (o)
     11  -6.8053E-02   2.2522E-02   9.0575E-02  -1.8115E-01   1.5863E-01   9.0101E-03  (o)
     12  -6.8691E-02   2.5699E-02   9.4390E-02  -1.8878E-01   1.6308E-01   9.6510E-03  (o)
     13  -6.8721E-02   2.6332E-02   9.5054E-02  -1.9011E-01   1.6378E-01   9.8228E-03  (o)
     14  -6.8721E-02   2.6327E-02   9.5049E-02  -1.9010E-01   1.6377E-01   9.8214E-03  (o)
     15  -6.8721E-02   2.6359E-02   9.5080E-02  -1.9016E-01   1.6380E-01   9.8300E-03  (o)
     16  -6.8721E-02   2.6349E-02   9.5070E-02  -1.9014E-01   1.6379E-01   9.8272E-03  (o)
     17  -6.8721E-02   2.6352E-02   9.5074E-02  -1.9015E-01   1.6380E-01   9.8282E-03  (o)
     18  -6.8721E-02   2.6351E-02   9.5072E-02  -1.9014E-01   1.6379E-01   9.8279E-03  (o)
     19  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     20  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     21  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     22  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     23  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     24  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     25  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     26  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     27  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     28  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     29  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     30  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03  (o)
     31  -6.8721E-02   2.6351E-02   9.5073E-02  -1.9015E-01   1.6379E-01   9.8280E-03
     32  -6.8842E-02   2.7693E-02   9.6535E-02  -1.9307E-01   1.6538E-01   1.0118E-02
     33  -7.6401E-02   2.7006E-02   1.0341E-01  -2.0681E-01   1.7981E-01   1.1018E-03
     34  -7.6461E-02   2.6598E-02   1.0306E-01  -2.0612E-01   1.7952E-01   2.8937E-04
     35  -7.6466E-02   2.6433E-02   1.0290E-01  -2.0580E-01   1.7937E-01   4.2615E-05
     36  -7.6466E-02   2.6420E-02   1.0289E-01  -2.0577E-01   1.7935E-01   1.1954E-05
     37  -7.6466E-02   2.6411E-02   1.0288E-01  -2.0575E-01   1.7934E-01   6.4220E-06
   --------------------------------------------------------------------------------------
   Converged: ||grad||=6.4220E-06 < vpq_tolgrs=1.0000E-05
   E_pol (eV):  -2.08075740E+00
     eps (eV):   4.88019624E+00

Note the (o) tag for optimization step with the orthogonality constraint. Also, check the gradient norm, which is indeed stuck at a constrained-induced stationary point before the constraint is lifted.

All results are gathered in a single VPQ NetCDF file. As before, we can visualize the convergence of our method using abiopen.py with the -e option (--expose):

The $E_\mathrm{pol}$ values for each state differ slighly:

$ grep "E_pol (eV)" teph4vpq_6.abo

   E_pol (eV):  -2.11263475E+00
   E_pol (eV):  -2.08075740E+00
   E_pol (eV):  -2.11071583E+00`

To visualize these states, one can compute the real-space charge distribution and displacements via Eqs. (\ref{eq:polaron_phi}),(\ref{eq:polaron_tau}).

In this tutorial, this is done with the following input

As usual, run this example with:

abinit teph4vpq_7.abi > teph4vpq_7.log 2> err &

This calculations requires the VPQ NetCDF file produced in the previous calculations, as well as the corresponding WFK, DDB, DVDB and GSTORE files. Note the new getvpq_filepath variable.

The code responsible for producing .xsf files with visualiziation data is activated via the eph_task -13 variable. The variable vpq_mesh_fact specifies the resolution downsampling of the charge distribution isosurface. Increasing its value will reduce memory requirements by losing some fine feautres of the isosurface. For small polarons vpq_mesh_fact 1 (no downsampling) is perfectly fine, but for the sake of tutorial we will use the downsampling anyaway.

This calculation has produced the following .xsf files:

$ ls *xsf

teph4vpq_7o_pstate_1_POLARON.xsf               teph4vpq_7o_pstate_2_POLARON.xsf               teph4vpq_7o_pstate_3_POLARON.xsf
teph4vpq_7o_pstate_1_POLARON_DISPL.xsf         teph4vpq_7o_pstate_2_POLARON_DISPL.xsf         teph4vpq_7o_pstate_3_POLARON_DISPL.xsf
teph4vpq_7o_pstate_1_POLARON_DISPL_VECTORS.xsf teph4vpq_7o_pstate_2_POLARON_DISPL_VECTORS.xsf teph4vpq_7o_pstate_3_POLARON_DISPL_VECTORS.xsf

For each state, we have

POLARON.xsf → polaron charge distribution in pristine geometry;
POLARON_DISPL.xsf → polaron charge distribution with distorted geometry;
POLARON_DISPL_VECTORS.xsf → pristine geometry with displacement vectors shown as forces.

One can use software that can read .xsf file format (e.g. VESTA) to visualize these outputs. The following figure shows the charge distribution of each polaronic state:

Note that state 1 and 3 are instances of a triply-degenerate hole polaron, and oriented along two different directions of the <100> cubic family. That’s why they are close in binding energy (~ 0.1 % differnece), and the difference may decrease if we ask for stricter convergence with vpq_tolgrs. These polarons have $D_{4h}$ point group symmetry.

However, state 2 possesses different symmetry ($D_{2h}$ point group) and is oriented along one of 6 direction of the <110> cubic family.

Having multiple and degenerate polaronic solutions for a single system is not surprising. This happens due to the degeneracy of electronic states and can be seen already at the level of generalized Frohlich model [Vasilchenko2024]. In ab initio calculations, this behaviour may be even more complex.

Also, note that these polarons localize on fluorine atoms and can occupy any fluorine site – that’s why we had to additionaly treat the translational invariance.

Tip

If you want to manually translate the polarons at visualization step, use the vpq_trvec variable.

Tip

Obtained solutions can be used as a starting point for calculations with other $\kk/\qq$-meshes. For this, one has to provide the VPQ file via getvpq_filepath, select a pstate via vpq_select and activate the interpolation using the vpq_interp tag.

This way, one can try to propagate a solution through a series of calculations on progressively increasing meshes and achieve consistent localization of polaron. However, as the potential energy surfaces of polarons possess multiple minima, it is still recommended to check the consistency of the solutions between meshes.

Finally, let’s look at the AbiPy capabilities to analyze individual contributions of the electronic bands and phonon modes to the polaron formation. This can be done with the VpqFile class which loads the VPQ NetCDF files and provides post-process tools. With the following Python scripts

from abipy.eph.vpq import VpqFile

ebands_file = "LiF_eph_varpeq/ebandso_GSR.nc"

vpq = VpqFile("teph4vpq_6o_VPQ.nc")

for polaron in vpq.polaron_spin:
    polaron.plot_ank_with_ebands(ebands_file, normalize=True)

from abipy.eph.vpq import VpqFile

ddb_file = "LiF_eph_varpeq/out_DDB"

vpq = VpqFile("teph4vpq_6o_VPQ.nc")

for polaron in vpq.polaron_spin:
    polaron.plot_bqnu_with_ddb(ddb_file)

we can generate the following figures (it may take some time):

For each state, left and right figure show the contribution of the electronic bands and phonon modes to the polaron formation, respectively. $\boldsymbol{A}$ and $\boldsymbol{B}$ coefficients are interpolated on $\kk/\qq$ paths – that’s why we see such dense sampling although our coefficients have been obtained on 3x3x3 $\kk/\qq$ meshes. These results are unconverged mainly due to the coarse mesh of the calculation, but still provide qualitative insights.

Note the difference between the $\boldsymbol{A}$ contribution for each state. Even for the degenerate states 1 and 3 it is different, which is mainly due to the fact that we used high-symmetry paths of the pristine crystal to plot them. With polarons, we have a symmetry-breaking effect, and one have to use special high-symmetry paths for each polaron. For states 1 and 3, for example, one of the three $\Gamma-X$ direction becomes inequivalent with the two others. For state 2, the same will hold for the set of $\Gamma-K$ directions.

This inequivalency can also be seen in the A-DOS, which have been computed with (red) and without symmetries (golden). If $\boldsymbol{A}$ coefficients have possesed the same symmetry as the crystal, these DOS would have been identical. However, the detailed symmetry analysis of polarons is beyond the scope of this tutorial.

Large electron polaron¶

Conceptually, calculations of large polarons follow the same flow as for small polarons. The main difference is in the size of $\kk/\qq$-meshes we have to employ as large polarons require large hosting supercells. Here, we will obtain large electron polaron in LiF on 13x13x13 $\kk/\qq$-mesh.

Again, we will need DDB, WFK and GSTORE files to launch the VarPEq optimization. The first two files are already available: DDB has been precomputed, and we also have obtained WFK on 13x13x13 $\kk$-mesh at the NSCF step.

First, let’s launch the gstore calculation with the input file

As usual, run this example with:

abinit teph4vpq_8.abi > teph4vpq_8.log 2> err &

While the calculation is running, let’s examine the input file. The main difference with our previous gstore runs, is the gstore_erange variable. We use it to select the energy window near the CBM, so matrix elements $g_{mn\nu}(\kk,\qq)$ will be computed ONLY for wavevectors $\kk$ and scattering $\kk+\qq$ inside this window. Energy filtering may significantly reduce memory requirements and computational time.

Of course, we could have used gstore_brange to select the entire conduction band. However, as large polarons are localized in reciprocal space, states contributing to the polaron formation, are confined within an energy window near the band edge.

Important

Of course, in production calculations, one has to perform convergence studies wrt gstore_erange. This can be done at a fixed $\kk/\qq$-mesh by increasing the energy window and monitoring the convergence of the polaron energy terms.

In this example, 1 eV window is a bit too small for convergence.

Tip

One can also use energy filtering at the step of computing the WFK files to achieve extra-dense meshes. See the corresponding section of the mobility tutorial.

Once the GSTORE file has been computed, we can launch the variational optimization (note that the vpq_pkind is now “electron”, mandatory):

abinit teph4vpq_9.abi > teph4vpq_9.log 2> err &

Exercise

Try using AbiPy to analyze to results stored in the teph4vpq_9o_VPQ.nc file the same we have demonstrated for small polarons. Do these results match your expectations? Keep in mind that our calculations are underconverged.

Of course, skip the extrapolation part, as we have not computed large polarons on progressively increasing meshes. Or do these calculations yourself if you are intrested!

Finally, we can visualize the large polaron via the following input

abinit teph4vpq_10.abi > teph4vpq_10.log 2> err &

This produces usual *xsf files that can be used for viauslization. If we open teph4vpq_10o_pstate_1_POLARON.xsf in VESTA, the following charge distribution can bee seen: