Calculation of polarons
Author: J. Lafuente-Bartolome' (v1, XX/01/2024) Revision: F. Giustino (v1.1, 07/23/2024)
In this notebook, we perform calculations of polarons. The theory and computational method are described in Phys. Rev. Lett. 122, 246403 (2019) and Phys. Rev. B 99, 235139 (2019).
Electrons and phonons are computed with Quantum ESPRESSO (QE), maximally-localized Wannier functions are computed with Wannier90, and polarons are computed with EPW.
Theory
[FG: We could remove this theory part to be consistent with the other notebooks]
Here, we briefly present the main concepts and equations that are solved by EPW to obtain polaron formation energy, wavefunction, and atomic displacements.
The ground state wave function \(\psi(\mathbf{r})\) and atomic displacements \(\Delta \tau_{\kappa\alpha p}\) of the polarons are found by minimizing the total energy functional of an excess electron added to a crystal. This minimization problem translates into the solution of the following coupled system of equations:
\[\hat{H}_{\mathrm{KS}}^{0} \, \psi(\mathbf{r})
+
\sum_{\kappa\alpha p} \frac{\partial V_{\mathrm{KS}}^{0}}{\partial \tau_{\kappa \alpha p}}
\Delta \tau_{\kappa \alpha p} \, \psi(\mathbf{r})
=
\varepsilon \, \psi(\mathbf{r}) ~,\]
\[\Delta \tau_{\kappa \alpha p}
=
-\sum_{\kappa'\alpha' p'} (C^{0})^{-1}_{\kappa\alpha p,\kappa'\alpha'p'}
\int d\mathbf{r} \, \frac{\partial V_{\mathrm{KS}}^{0}}{\partial \tau_{\kappa'\alpha'p'}}
\, |\psi(\mathbf{r})|^2 ~.\]
In these expressions, \(\tau_{\kappa\alpha p}\) represents the Cartesian coordinate of atom \(\kappa\) in unit cell \(p\) along the direction \(\alpha\), \(C^{0}_{\kappa\alpha p,\kappa'\alpha'p'}\) is the matrix of interatomic force constants, and \(\hat{H}_{\mathrm{KS}}^{0}\) and \(V_{\mathrm{KS}}^{0}\) represent the Kohn-Sham Hamiltonian and the self-consistent potential, respectively. The superscript \(^{0}\) indicates that these quantities are evaluated in the ground state without extra electron. The integral is performed over a Born-Von Karman supercell of the crystal containing \(N_p\) unit cells. We will refer to \(\varepsilon\) as the polaron eigenvalue.
By expanding the polaron wave function in terms of the single-particle Kohn-Sham states \(\psi_{n\mathbf{k}}\) with eigenvalues \(\varepsilon_{n\mathbf{k}}\):
\[\psi(\mathbf{r}) = \frac{1}{\sqrt{N_p}}
\sum_{n\mathbf{k}} A_{n\mathbf{k}} \psi_{n\mathbf{k}} ~,\]
and the atomic displacements in terms of the phonon eigenmodes \(e_{\kappa\alpha,\nu}(\mathbf{q})\) with frequencies \(\omega_{\mathbf{q}\nu}\):
\[\Delta\tau_{\kappa\alpha p } = -\frac{2}{N_p}
\sum_{\mathbf{q}\nu} B^{*}_{\mathbf{q}\nu}
\left( \frac{\hbar}{2M_\kappa \omega_{\mathbf{q}\nu}} \right)^{1/2}
e_{\kappa\alpha,\nu}(\mathbf{q}) \, e^{i\mathbf{q}\cdot\mathbf{R}_p},\]
where \(M_\kappa\) is the mass of the atom \(\kappa\) and \(\mathbf{R}_p\) is the lattice vector of the unit cell \(p\), we can transform the first two equations into a coupled set of equations for the expansion coefficients in reciprocal space:
\[\frac{2}{N_p} \sum_{\mathbf{q}m\nu} B_{\mathbf{q}\nu}
\, g_{mn\nu}^{*}(\mathbf{k},\mathbf{q}) \, A_{m\mathbf{k+q}}
=
(\varepsilon_{n\mathbf{k}}-\varepsilon) A_{n\mathbf{k}} ~,\]
\[B_{\mathbf{q}\nu} = \frac{1}{N_p}
\sum_{mn\mathbf{k}} A^{*}_{m\mathbf{k+q}}
\frac{g_{mn\nu}(\mathbf{k},\mathbf{q})}{\hbar\omega_{\mathbf{q}\nu}} A_{n\mathbf{k}} ~.\]
The polaron formation energy \(\Delta E_{f}\), defined as the energy required to trap a conduction band state with eigenvalue \(\varepsilon_{\mathrm{CBM}}\) into a localized polaron, can be obtained from the expansion coefficients that solve the coupled set of equations by:
\[\label{eq:eq3}
\Delta E_{f}
=
\frac{1}{N_p} \sum_{n\mathbf{k}} |A_{n\mathbf{k}}|^2
(\varepsilon_{n\mathbf{k}}-\varepsilon_{\mathrm{CBM}})
-
- \frac{1}{N_p} \sum_{\mathbf{q}\nu} |B_{\mathbf{q}\nu}|^2 \hbar\omega_{\mathbf{q}\nu} ~.\]
We will refer to the first and second terms on the right hand side as the electron and phonon parts of the formation energy, respectively.
From these expressions, we observe that the necessary ingredients to solve the polaron equations are the Kohn-Sham eigenvalues, phonon frequencies, and electron-phonon matrix elements on the \(\mathbf{k}\)- and \(\mathbf{q}\)-grids corresponding to the equivalent Born-Von Karman supercell in which the first two equations are defined. In order to obtain the formation energy of an isolated polaron, we will need to solve the coupled set of equations in increasingly denser grids and extrapolate the results to the infinite supercell limit.
Set up working environment
Load required modules:
[7]:
import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import MultipleLocator,FormatStrFormatter
import os
import pandas as pd
#
import sys
print(str(os.getcwd()))
sys.path.insert(0,str(os.getcwd())+'/../')
from EPWpy import EPWpy
from EPWpy.plotting import plot_polaron
from EPWpy.utilities import EPW_util
/workspace/Sabya/codes/EPW_py/TACC/EPWpy_release/epwpy/notebooks_basic
Set paths to relevant directories:
[8]:
# Set QE binaries directory
QE = '/workspace/Sabya/codes/EPW_develop/q-e/bin'
prefix='lif'
Set the number of cores to be used in the calculations:
[9]:
# Maximum number of cores to be used
cores = 16
print('Maximum number of cores to be used:', cores)
Maximum number of cores to be used: 16
Define general calculation parameters to be used throughout the workflow:
[4]:
lif=EPWpy.EPWpy({'prefix':prefix,
'restart_mode':'\'from_scratch\'',
'calculation':'\'scf\'',
'ibrav':2,
'nat':2,
'ntyp':2,
'atomic_species':['Li', 'F'],
'atomic_pos':np.array([[0.0, 0.0, 0.0], [-0.5, 0.5, 0.5]]), # in alat
'mass':[6.941, 18.9984],
'atoms':['Li', 'F'],
'ecutwfc':'80',
'celldm(1)':'7.67034',
'verbosity':'high',
'pseudo_auto':True},
code=QE,
env='ibrun',
system='lif')
lif.run_serial = True
lif.env = 'mpirun'
#'pseudo_dir':'\''+str(pseudo)+'\''},
# Print relevant info
os.system('module list')
-*#*- ...............-
.+*= .+%*-=%%: .=#*- -===============-:.
:*%=*%%- *%* #%* :+%+-%%+ .:. -=. :==-.
-%S -%%*: :#%. -%%-. -##: #%* -=. :==-
.. .%S: +%%%%*. :*%%%#= %%= -=. :==-
:=#%%*- .#S- .. .%%= :*#*: -=. :==-
-%S:.=#%%*==%# *%%=::=##-+%%. . .=-. :==- .=
:%%- .-+++: -+##*=. =%S :-::==: .==- --.
*%# #%+ -=--:. .----:.
:%%- -%S.
.-=*####SS%#########: -###########*+: +#######= =%SS####+::.
=+#**%SSSSSSSSSSSSSSSS= =SSSSSSSSSSSSSSS= #SSSSSSS*. +SSSSSSSS%%%%-
*%% =SSS.. .SSS= SSS=. .:+SSS+ -SSS:. .-::. .SSS+. #%*
#%#. =SSS. *S# *S%: SSS= %SS%. .SSS= #SSS%. :SSS: =%#.
*%%: =SSS#*#SSS- SSS= .+SSS+. %SS*.=SSSSS+ +SS%.:+%+
+%%:=SSSSSSSSS- SSSSSSSSSSSSS=. +SSS:SSS%SSS:%SS*-#%=
....#S==SSS:..SSS: =+- SSS%######+=. -SSS%SS%.#SS%SSS==%%=.
.:+##%%#*- =SSS. ::. :SSS. SSS=. SSSSSS:..SSSSSS: :*%%%*=-
#%+. -+#SSS*+++++++*SSS: .=+SSS#++++: %SSSS+. =SSSS%. :-+#%%+
#%* .SSSSSSSSSSSSSSSSSS: *SSSSSSSSSSS. +SSS%. %SSS*. . .%%=
=%S :::::::::::::::::. .:::::::::. :::. :::. =+=-.#%+
-%S: =+===#S:
==*------------------------------=========+++++++++++++++++++++++========++-+##
=+++++++++++*******++++++++++++++++========------------------==========+++++++-
-- -- -- -- -- -- -- -- -- -- -- -- Structure Info -- -- -- -- -- -- -- -- -- -- -- -- --
2
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Li/Li.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Li/Li_r.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Li/Li-sp_r.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Li/Li-d_r.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Li/Li-s_r.upf
pseudo found at pseudodojo : ONCVPSP-PBE-FR-PDv0.4/Li-s_r.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/F/F.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/F/F_r.upf
pseudo found at pseudodojo : ONCVPSP-PBE-FR-PDv0.4/F_r.upf
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
sh: 1: module: not found
[4]:
32512
Workflow
We will consider a hole polaron in LiF as an example.
Preliminary scf, ph, and nscf calculations
[7]:
# Prepare scf input file
lif.scf(electrons={'conv_thr':'1.0d-12'}, kpoints={'kpoints':[[6, 6, 6]]}, name='scf')
lif.prepare(4,type_run='scf')
# Run scf calculation
lif.run(4)
-- -- -- -- -- -- -- -- -- -- -- Calculation: scf -- -- -- -- -- -- -- -- -- -- --
Running scf |████████████████████████████████████████| in 1.6s (1.91/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[8]:
# Prepare ph input file
#lif.ph(phonons={'fildyn':'\'lif.dyn\'','nq1':6,'nq2':6,'nq3':6,'fildvscf':'\'dvscf\''})
lif.ph(phonons={'nq1':6,
'nq2':6,
'nq3':6})
# Run ph calculation
lif.prepare(20,type_run='ph')
lif.run(4,type_run='ph')
-- -- -- -- -- -- -- -- -- -- -- Calculation: ph -- -- -- -- -- -- -- -- -- -- --
Running ph |████████████████████████████████████████| in 15:27.1 (0.00/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[9]:
# Prepare nscf input file
lif.nscf(system={'nbnd':15},kpoints={'grid':[6, 6, 6],'kpoints_type': 'crystal'})
lif.prepare(20,type_run='nscf')
# Run nscf calculation
lif.run(4,type_run='nscf')
-- -- -- -- -- -- -- -- -- -- -- Calculation: nscf -- -- -- -- -- -- -- -- -- -- --
Running nscf |████████████████████████████████████████| in 14.1s (0.10/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Coarse mesh epw calculations
[6]:
# Prepare epw input for coarse mesh calculation
lif.reset()
lif.epw(epwin={'proj':['\'F : p\''],
'elph':'.true.',
'epwwrite':'.true.',
'epwread':'.false.',
'wannier_plot':'.true.',
'band_plot':'.true.',
'filkf':'\'./path.kpt\'',
'filqf':'\'./path.kpt\'',
'num_iter':500,
'epbwrite':'.false.',
'nbndsub':'3',
'bands_skipped':'\'exclude_bands = 1:2, 6:15\'',
'wannier_plot':'.true.',
'wannier_plot_supercell':'6, 6, 6',
'lpolar':'.true.',
'eig_read':'.false.',
'wannierize':'.true.'},
name='epw1')
# Generate filkf if needed
W=[0.5,0.75,0.25]
L=[0.0,0.5,0.0]
G=[0.0,0.0,0.0]
X=[0.5,0.5,0.0]
K=[0.375,0.75,0.375]
sympoints=[W,L,G,X,W,K]
#
lif.filkf(path=sympoints, length=[41, 41, 41, 41, 41], name='path.kpt')
#
lif.prepare(1,type_run='epw1')
obtaining nscf and ph attributes
nscf ph
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: .true. (related to epsil in ph)
[17]:
print(lif.code)
#lif.env='mpirun -np 4'
lif.run(4,type_run='epw1')
/workspace/Sabya/codes/EPW_develop/q-e/bin
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw1 -- -- -- -- -- -- -- -- -- -- --
Running epw1 |████████████████████████████████████████| in 5:40.8 (0.00/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Interpolation to fine mesh and solution of polaron equations
[7]:
# Prepare epw input for fine mesh interpolation and polaron calculation
lif.epw(epwin={'bands_skipped':'\'exclude_bands = 1:2, 6:15\'',
'nbndsub':3,
'plrn':'.true.',
'restart_plrn':'.false.',
'type_plrn':1,
'init_plrn':1,
'init_sigma_plrn':1.0,
'niter_plrn':500,
'conv_thr_plrn':1.0E-4,
'ethrdg_plrn':1.0E-6,
'adapt_ethrdg_plrn':'.true.',
'init_ethrdg_plrn':1.0E-4,
'nethrdg_plrn':20,
'io_lvl_plrn':0,
'nkf1':6,
'nkf2':6,
'nkf3':6,
'nqf1':6,
'nqf2':6,
'nqf3':6,
'lpolar':'.true.',
'eig_read':'.false.'},
name='epw2')
#
lif.prepare(1,type_run='epw2')
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: .true. (related to epsil in ph)
[8]:
# Run
lif.run(4,type_run='epw2')
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw2 -- -- -- -- -- -- -- -- -- -- --
Running epw2 |████████████████████████████████████████| in 10.2s (0.14/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Post-processing: Polaron wave function and single-particle expansion coefficients
[9]:
# Prepare epw input for fine mesh interpolation and polaron calculation
lif.epw(epwin={'bands_skipped':'\'exclude_bands = 1:2, 6:15\'',
'plrn':'.true.',
'restart_plrn':'.true.',
'type_plrn':1,
'cal_psir_plrn':'.true.',
'step_wf_grid_plrn':4,
'interp_Ank_plrn':'.true.',
'interp_Bqu_plrn':'.true.',
'filkf':'\'../path.kpt\'',
'filqf':'\'../path.kpt\'',
'nkf1':6,
'nkf2':6,
'nkf3':6,
'nqf1':6,
'nqf2':6,
'nqf3':6,
'eig_read':'.false.',
'lpolar':'.true.'},
name='epw3')
lif.prepare(1,type_run='epw3')
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: .true. (related to epsil in ph)
cp: cannot stat '../path.kpt': No such file or directory
[10]:
# Run post-processing
lif.run(4,type_run='epw3')
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw3 -- -- -- -- -- -- -- -- -- -- --
Running epw3 |████████████████████████████████████████| in 5.3s (0.29/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[11]:
# Plot Ank coefficients
print(os.getcwd())
prefix=lif.prefix
print(prefix+'/path.kpt')
plot_polaron.plot_Ank_Bqv('Ank', prefix+'/epw/path.kpt', sympoints, prefix+'/epw/Ank.band.plrn')
/workspace/Sabya/codes/EPW_py/TACC/EPWpy_release/epwpy/notebooks_basic
lif/path.kpt
[12]:
# Plot Bqv coefficients
plot_polaron.plot_Ank_Bqv('Bqv', prefix+'/path.kpt', sympoints, prefix+'/epw/Bmat.band.plrn')
OPTIONAL: Plot wavefunction; needs installation of VESTA (TODO improve this, plot with python?)
[6]:
Data=EPW_util.read_psir_plrn('./lif/epw/psir_plrn.xsf')
plot_polaron.plot_psir_plrn(Data)
['-12.1769074', '0.0000000', '12.1769074']
['0.0000000', '12.1769074', '12.1769074']
['-12.1769074', '12.1769074', '0.0000000']
432
Perform polaron calculations on increasingly dense meshes, and extrapolate to \(N_p\rightarrow\infty\) limit
[21]:
# List of k-point meshes to consider
klist = [8, 10, 12]
for k in klist:
lif.epw(epwin={'bands_skipped':'\'exclude_bands = 1:2, 6:15\'',
'plrn':'.true.',
'restart_plrn':'.true.',
'type_plrn':1,
'cal_psir_plrn':'.true.',
'step_wf_grid_plrn':4,
'interp_Ank_plrn':'.true.',
'interp_Bqu_plrn':'.true.',
'filkf':'\'../path.kpt\'',
'filqf':'\'../path.kpt\'',
'nkf1':str(k),
'nkf2':str(k),
'nkf3':str(k),
'nqf1':str(k),
'nqf2':str(k),
'nqf3':str(k),
'eig_read':'.false.'},
name='epw3.'+str(k))
#
lif.prepare(1,type_run='epw3',name = 'epw',infile='epw3.'+str(k)+'.in')
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: None (related to epsil in ph)
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: None (related to epsil in ph)
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- Info -- -- -- -- -- -- -- -- -- -- -- -- -- --
Based on previous pw and ph calculations some parameters are set below
lpolar: None (related to epsil in ph)
cp: cannot stat '../path.kpt': No such file or directory
cp: cannot stat '../path.kpt': No such file or directory
cp: cannot stat '../path.kpt': No such file or directory
[22]:
# Run calculations
lif.run_serial=True
for k in klist:
lif.run(16, type_run='epw3', infile='epw3.'+str(k))
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw3.8 -- -- -- -- -- -- -- -- -- -- --
Running epw3.8 |████████████████████████████████████████| in 1:01.6 (0.02/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw3.10 -- -- -- -- -- -- -- -- -- -- --
Running epw3.10 |████████████████████████████████████████| in 1:01.2 (0.02/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw3.12 -- -- -- -- -- -- -- -- -- -- --
Running epw3.12 |████████████████████████████████████████| in 1:01.5 (0.02/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[23]:
## Linear fit and extrapolate to Np->\infty
# Set unit cell volume in Bohr^3 (to convert Nk to inverse supercell size)
ucell_volume = 112.8044
# Read formation energies from output files
Eform = []
for k in klist:
lif.epw_file = 'epw.k'+str(k)
Eform.append(lif.Eform)
Nk = np.array(klist)
Eform = np.array(Eform)
# Plot
plot_polaron.plot_EvsNk(ucell_volume, Nk, Eform)
Extrapolation to isolated polaron formation energy = -2.054 eV
[ ]: