This tutorial explains how to compute the phonon linewidths induced by the electron-phonon (e-ph) interactionin metallic systems and how to use the Eliashberg spectral function \alpha^2F and the McMillan equationto estimate the superconducting critical temperature T_c within the isotropic Eliashberg formalism.
We start by presenting the basic equations implemented in the code and their connection with the ABINIT variables.Then we discuss how to run isotropic T_c-calculations and how to perform typical convergence studiesfor MgB_2, a well-known phonon-mediated superconductor with T_c = 39 K.For a more complete theoretical introduction, see [Giustino2017] and references therein.
It is assumed the user has already completed the two tutorials RF1 and RF2,and that he/she is familiar with the calculation of ground state and vibrational properties in metals.The user should have read the fourth lesson on Al as wellas the introduction page for the EPH code before running these examples.
This lesson should take about 1.5 hour.
Formalism and connection with the implementation¶
Due to the interaction with electrons, vibrational energies get renormalized and phononsacquire a finite lifetime.These many-body effects are described by the phonon-electron self-energy:
\Pi_{\qq\nu}(\ww, T) = \dfrac{2}{N_\kk} \sum_{mn\kk} |g_{mn\nu}(\kk, \qq)|^2\dfrac{f_\nk - f_{m\kq}}{\ee_{\nk} - \ee_{m\kq} -\ww - i\eta^+}
where only contributions up to the second order in the e-ph vertex g have been included.The sum over the electron wavevector \kk is performed over the first BZ, \eta is a positive real infinitesimaland \gkkp are the e-ph matrix element discussed in the EPH introduction.The self-energy depends on the temperature via the Fermi-Dirac distribution function f(\ee, T)and the factor two accounts for spin degeneracy (henceforth we assume a non-magnetic systemwith scalar wavefunctions i.e. nsppol == 1 and nspinor == 1thus the electron energies \ee and the e-ph matrix elements g do not depend on the spin).
The real part of \Pi_\qnu gives the correction to the vibrational energies due to e-ph interactionwhile the phonon linewidth \gamma_{\qq\nu} (full width at half maximum) is twicethe imaginary part \Pi^{''}_\qnu evaluated at the “bare” phonon frequency \ww_\qnu as computed with DFPT:
\gamma_{\qq\nu}(T) = 2\, \Pi^{''}_\qnu(\ww=\ww_\qnu, T).
Using the Sokhotski–Plemelj theorem:
\lim_{\eta \rightarrow 0^+} \dfrac{1}{x \pm i\eta} = \mp i\pi\delta(x) + \mathcal{P}(\dfrac{1}{x})
the imaginary part of \Pi can we rewritten as:
\Pi^{''}_{\qq\nu}(\ww_\qnu, T) = \dfrac{2\pi}{N_\kk} \sum_{mn\kk} |g_{mn\nu}(\kk, \qq)|^2(f_\nk - f_{m\kq})\,\delta(\ee_{\kk n} - \ee_{m\kpq} -\ww_\qnu)
where the delta function enforces energy conservation in the scattering process.Because phonon energies are typically small compared to electronic energies,the energy difference \ee_{\kk n} - \ee_{m\kpq} is also small hence it is possible toapproximate the difference in the occupation factors with:
f_\nk - f_{m\kq} \approx f'(\ee_\nk) (\ee_{\kk n} - \ee_{m\kpq}) = -f'(\ee_\nk) \ww_\qnu
At low T, the derivative of the Fermi-Dirac occupation function is strongly peaked aroundthe Fermi level thus we can approximate it with:
f'(\ee_\nk) \approx -\delta(\ee_\nk - \ee_F)
By neglecting \ww_\qnu in the argument of the delta, we obtain theso-called double-delta approximation (DDA) for the phonon linewidth [Allen1972]:
\begin{equation}\label{eq:DDA_gamma}\gamma_{\qq\nu} = \dfrac{4\pi \ww_{\qq\nu}}{N_\kk} \sum_{mn\kk}|g_{mn\nu}(\kk, \qq)|^2 \delta(\ee_{\kpq m} -\ee_F) \delta(\ee_{\kk n} -\ee_F)\end{equation}
For a given phonon wavevector \qq, the DDA restricts the BZ integration toelectronic transitions between \kk and \kq states on the Fermi surface (FS)hence very dense \kk-meshes are needed to converge \gamma_{\qq\nu}.The convergence rate is expected to be slower than that required by the electron DOS per spin at \ee_F:
g(\ee_F) = \dfrac{1}{N_\kk} \sum_{n\kk} \delta(\ee_{n\kk} - \ee_F)
in which a single Dirac delta is involved.
Note
The DDA breaks down at small \qq where the phonon frequency cannot be neglected.Note also that the DDA predicts the linewidths to not depend on T.
It is also worth stressing that there is another important contribution to the phonon lifetimes induced bynon-harmonic terms in the Taylor expansion of the Born-Oppenheimer energy surface around the equilibrium point.Within the framework of many-body perturbation theory, these non-harmonic terms leadto phonon-phonon scattering processes that can give a substantial contribution to the lifetimes,especially at “high” T.In the rest of the tutorial, however, non-harmonic terms will be ignored and we will be mainly focusingon the computation of the imaginary part of the phonon self-energy \Pi in the harmonic approximation.
At the level of the implementation, the eph_intmeth input variableselects the technique for integrating the double delta over the FS:eph_intmeth == 2 (default) activates the optimized tetrahedron scheme [Kawamura2014]as implemented in the libtetrabz librarywhile eph_intmeth == 1 replaces the Dirac distribution with a Gaussian function of finite width.In the later case, one can choose between a constant broadening \sigmaspecified in Hartree by eph_fsmear or an adaptive Gaussian scheme (activated if eph_fsmear < 0)in which a state-dependent \sigma_\nk is automatically computed fromthe electron group velocities v_\nk [Li2015].Note that the adaptive Gaussian scheme is internally used by the code also when the optimized tetrahedron methodis employed since the DDA is ill-defined for \qq = \Gamma.
The value of the Fermi level \ee_F is automatically computed from the KS eigenvalues storedin the input WFK file according to the two input variables occopt and tsmear.These variables are usually equal to the ones used for the GS/DFPT calculationsalthough it is possible to change the value of \ee_F during an EPH calculation usingthree (mutually exclusive) input variables: eph_fermie, eph_extrael and eph_doping.This may be useful if one wants to study the effect of doping within the rigid band approximation.
The sum over bands in Eq.\ref{eq:DDA_gamma} is restricted to the Bloch states within asymmetric energy window of thickness eph_fsewin around \ee_F.This value should be chosen according to the FS integration scheme eph_intmeth and the value ofeph_fsmear (if standard Gaussian).Additional details are given in the tutorial.The \kk-mesh for electrons is defined by the input variables ngkpt, nshiftk and shiftk.These parameters must be equal to the ones used to generate the input WFK file passed to the EPH code.
The code computes \gamma_{\qq\nu} for all the \qq-point in the IBZ associatedto the eph_ngqpt_fine \qq-mesh and the DFPT potentials are interpolated startingfrom the ab-initio ddb_ngqpt \qq-mesh associated to the input DVDB/DDB files.If eph_ngqpt_fine is not specified, eph_ngqpt_fine == ddb_ngqpt is assumedand no interpolation of the scattering potentials in \qq-space is performed.
Once the phonon linewidths \gamma_{\qq\nu} are known in the IBZ, ABINIT computesthe isotropic Eliashberg function defined by:
\alpha^2F(\ww) =\dfrac{1}{N_\qq\, N_F} \sum_{\qq\nu} \dfrac{\gamma_{\qq\nu}}{\ww_{\qq\nu}} \delta(\ww - \ww_{\qq \nu})
where N_F is the density of states (DOS) per spin at the Fermi level.The Eliashberg function can be equivalently expressed as:
\alpha^2F(\ww) =\dfrac{1}{N_\qq N_\kk\, N_F} \sum_{\kk\qq mn \nu}|g_{mn\nu}(\kk, \qq)|^2 \delta(\ee_{\kpq m} -\ee_F) \delta(\ee_{\kk n} -\ee_F) \delta(\ww - \ww_{\qq \nu}).
From a physical perspective, \alpha^2F(\ww) gives the strength by which a phononof energy \hbar\ww scatters electronic states on the FS (remember that ABINIT uses atomic units hence \hbar = 1).This quantity is accessible in experiments and experience has shown that\alpha^2F(\ww) is qualitatively similar to the phonon DOS F(\ww):
F(\ww) = \dfrac{1}{N_\qq} \sum_{\qq\nu} \delta(\ww - \ww_{\qq \nu}).
This is not surprising as the equation for \alpha^2F(\ww) resembles the one for F(\ww),except for the weighting factor \frac{\gamma_{\qq\nu}}{\ww_{\qq\nu}}.This also explains the \alpha^2F notation: the Eliashberg function can be seenas the phonon DOS times the positive frequency-dependent prefactor \alpha^2(\ww).
Warning
Converging \alpha^2F(\ww) usually requires very fine \kk- and \qq-meshes.This is especially true if the FS significantly deviates from the free electron picture (spherical FS)and/or the e-ph coupling is strongly anisotropic as in MgB_2.The \kk-mesh affects the quality of the individual \gamma_\qnu linewidthswhereas dense \qq-meshes may be needed to resolve the fine detauls and/orsample regions in \qq-space where the coupling is strong.
The technique used to compute \alpha^2F(\ww) is defined by ph_intmeth (note the ph_ prefix instead of eph_).Both the Gaussian (ph_intmeth = 1 with ph_smear smearing) andthe linear tetrahedron method by [Bloechl1994] (ph_intmeth = 2, default) are implemented.The \alpha^2F(\ww) function is evaluated on a linear mesh of step ph_wstep covering the entirerange of phonon frequencies.The total e-ph coupling strength \lambda (a dimensionless measure of the average strength of the e-ph coupling)is defined as the first inverse moment of \alpha^2F(\ww):
\lambda = \int \dfrac{\alpha^2F(\ww)}{\ww}\dd\ww = \sum_{\qq\nu} \lambda_{\qq\nu}
where we have introduced the mode-dependent contributions:
\lambda_{\qq\nu} = \dfrac{\gamma_{\qq\nu}}{\pi N_F \ww_{\qq\nu}^2}
Warning
Due to the 1/\omega factor, “low-energy” modes are expected to contribute more to \lambdathan “high-energy modes”.On the other hand, \alpha^2F(\ww) should go to zero as \ww^2 for \ww \rightarrow 0so that the integrand function is finite at the acoustic limit.Obviously, if the system is not dynamically stable (“negative” frequencies in the BZ)\lambda will explode but this does not mean you have found a room-temperature superconductor!
In principle, T_c can be obtained by solving the isotropic Eliashberg equations for the superconductinggap (see e.g. [Margine2013]) but many applications prefer to bypass the explicit solutionand estimate T_c using the semi-empirical McMillan equation [McMillan1968]in the improved version proposed by [Allen1975]:
T_c = \dfrac{\ww_{log}}{1.2} \exp \Biggl [ \dfrac{-1.04 (1 + \lambda)}{\lambda ( 1 - 0.62 \mu^*) - \mu^*} \Biggr ]
where \mu^* describes the screened electron-electron interaction and\ww_{\text{log}} is the logarithmic average of the phonon frequencies defined by:
\ww_{\text{log}} = \exp \Biggl [ \dfrac{2}{\lambda} \int \dfrac{\alpha^2F(\ww)}{\ww}\log(\ww)\dd\ww \Biggr ]
In pratical applications, \mu^* is treated as an external parameter,typically between 0.1 and 0.2, that is adjusted to reproduce experimental results.The default value of eph_mustar is 0.1.
Before concluding this brief theoretical introduction, we would like tostress that the present formalism assumes a single superconducting gapwith a weak dependence on \kk so that it is possible to average the anisotropic equations over the FS.This approximation becomes valid in the so-called dirty-limit (metals with impurities)as impurities tend to smear out the anisotropy of the superconducting gap.A more detailed discussion about isotropic/anisotropic formulations can be found in [Margine2013].
Getting started¶
Note
Supposing you made your own installation of ABINIT, the input files to run the examplesare 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 definesome 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 lineexport PATH=$ABI_HOME/src/98_main/:$PATH # Do not change this line: path to executableexport ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests direxport ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/ # Do not change this line: path to pseudos dirExamples in this tutorial use these shell variables: copy and pastethe 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.shThe ‘export PATH’ line adds the directory containing the executables to your PATHso 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*) andcopy 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_eph4isotc in $ABI_TESTS/tutorespfn/Input?
cd $ABI_TESTS/tutorespfn/Inputmkdir Work_eph4isotccd Work_eph4isotcIn this lesson we prefer to focus on e-ph calculations and the associated convergence studies.For this reason, we rely on pre-computed DEN.nc, DDB and DFPT potentialsto bypass both the GS and the DFPT part.The DEN.nc file will be used to perform NSCF computations on arbitrarily dense \kk-meshes while theDFPT POT.nc files will be merged with the mrgdv utility to produce the DVDB database required by the EPH code.
Note that these files are not shipped with the official ABINIT tarball as they are relatively large.In order to run the examples of this tutorial, you need to download the files fromthis github repository.If git is installed on your machine, you can easily fetch the entire repository with:
git clone https://github.com/abinit/MgB2_eph4isotc.gitAlternatively, use wget:
wget https://github.com/abinit/MgB2_eph4isotc/archive/master.zipor curl:
curl -L https://github.com/abinit/MgB2_eph4isotc/archive/master.zip -o master.zipor simply copy the tarball by clicking the “download button” in the github web page,then unzip the archive and rename the directory with:
unzip master.zipmv MgB2_eph4isotc-master MgB2_eph4isotcWarning
The directory with the precomputed files must be located inside Work_eph4isotcand must be named MgB2_eph4isotc.The reason is that all the input files and examples of this tutorial read data from external filesspecified in terms of relative paths.
The AbiPy script used to perform the GS + DFPT steps is availablehere.In order to facilitate comparison with previous studies,we use norm-conserving pseudopotentials with a cutoff energy ecut of 38 Haand the experimental parameters for hexagonal MgB_2 (a = 5.8317 and c/a= 1.1416).All the calculations are performed with a 12x12x12 ngkpt Gamma-centered \kk-grid for electrons (too coarse),and the Marzari smearing (occopt = 4) with tsmear = 0.02 Ha.The DFPT computations is done for 12 irreducible \qq-points correspondingto a \Gamma-centered 4x4x4 \qq-mesh (again, too coarse as we will see in the next sections).
Please keep in mind that several parameters have been tuned in order to reach a reasonable compromise between accuracyand computational cost so do not expect the results obtained at the end of the lesson to be fully converged.It is clear that, in real life, one should start from convergence studies forlattice parameters and vibrational properties as a function of the \kk-mesh and tsmearbefore embarking on EPH calculations.
Merging the DFPT potentials¶
To merge the DFPT potential files, copy the first input file in Work_eph4isotc with:
cp $ABI_TESTS/tutorespfn/Input/teph4isotc_1.abi .and execute the mrgdv tool using:
mrgdv < teph4isotc_1.abiThe first line in teph4isotc_1.abi specifies the name of the output DVDB file, followed by thenumber of partial DFPT POT files and the full list of files we want to merge:
This step produces the teph4isotc_1_DVDB file that will be used in the next examples.Executing:
mrgdv info teph4isotc_1_DVDBshows that all the independent atomic perturbations are available:
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=independent, found=Yes 4) idir= 1, ipert= 2, type=independent, found=Yes 5) idir= 2, ipert= 2, type=symmetric, found=No 6) idir= 3, ipert= 2, type=independent, found=Yes 7) idir= 1, ipert= 3, type=symmetric, found=No 8) idir= 2, ipert= 3, type=symmetric, found=No 9) idir= 3, ipert= 3, type=symmetric, found=No All the independent perturbations are available DoneAnalyzing electronic and vibrational properties¶
Before proceeding with e-ph calculations, it is worth spending some time to analyzethe structural, electronic and vibrational properties of MgB_2 in more detail.We will be using the AbiPy scripts to post-process the data stored in the precomputed netcdf files.
To visualize the crystalline structure with e.g. vesta ,use the abiview.py script and the structure command.
abiview.py structure MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/run.abi -a vestaFull Formula (Mg1 B2)Reduced Formula: MgB2abc : 3.086000 3.086000 3.523000angles: 90.000000 90.000000 120.000000Sites (3) # SP a b c--- ---- -------- -------- --- 0 Mg 0 0 0 1 B 0.333333 0.666667 0.5 2 B 0.666667 0.333333 0.5Visualizing structure with: vestaWriting data to: .xsf with fmt: xsfExecuting MacOSx open command: open -a vesta --args /Users/gmatteo/git_repos/abinit_rmms/_build/tests/.xsfTip
Other graphical applications are supported, see abiview.py structure --help.Note that AbiPy assumes that these tools are already installed and the binariescan be found in $PATH.The same AbiPy command can be used with ABINIT input files as well as output filesin netcdf format e.g. GSR.nc.
MgB_2 crystallizes in the so-calledAlB_2 prototype structurewith Boron atoms forming graphite-like (honeycomb) layers separated by layers of Mg atoms.This structure may be regarded as that of completely intercalated graphite with C replaced by B.Since MgB_2 is formally isoelectronic to graphite,its band dispersion is expected to show some similarity to thatof graphite and graphite intercalation compounds.
Now use the abiopen.py script to plot the electronic bands stored in the GSR file producedby the NSCF calculation (the second task of the first Work i.e. w0/t1):
abiopen.py MgB2_eph4isotc/flow_mgb2_phonons/w0/t1/outdata/out_GSR.nc -eThis command produces the following figures:
The electronic properties of MgB_2 are intensively discussed in literature, see e.g. [Mazin2003].The dispersion is quite similar to that of graphite with three bonding \sigmabands corresponding to in-plane sp^2 hybridization in the boron layer and two \pi bands(bonding \pi and anti-bonding \pi^*) formed by B p_z orbitals.The lowest band in the band plot is a fully occupied \sigma band.The other two \sigma bands are incompletely filled and correspond tothe relatively flat states located slightly above \ee_F along the \Gamma-A segment.Ab-initio calculations showed that these states contribute the most to the e-ph coupling.
The other two bands crossing at the K point at around 1eV are \pi-\pi^* bands.
To compute the so-called fatbands, one can use perforn a NSCF calculation with a \kk-path and prtdos 3.This is left as optional exercise.Fatbands plots produced with AbiPy starting from the FATBANDS.nc file are availablehere.
Tip
The high-symmetry \kk-path has been automatically computed by AbiPy when generating the Flow.To get the explicit list of high-symmetry \kk-point, use the abistruct.py scriptwith the kpath command:
abistruct.py kpath mgb2_DEN.nc...# K-path in reduced coordinates: ndivsm 10 kptopt -11 kptbounds +0.00000 +0.00000 +0.00000 # $\Gamma$ +0.50000 +0.00000 +0.00000 # M +0.33333 +0.33333 +0.00000 # K +0.00000 +0.00000 +0.00000 # $\Gamma$ +0.00000 +0.00000 +0.50000 # A +0.50000 +0.00000 +0.50000 # L +0.33333 +0.33333 +0.50000 # H +0.00000 +0.00000 +0.50000 # A +0.50000 +0.00000 +0.50000 # L +0.50000 +0.00000 +0.00000 # M +0.33333 +0.33333 +0.00000 # K +0.33333 +0.33333 +0.50000 # HThe value of the DOS at the Fermi level plays a very important role when discussing superconducting properties.Actually this is one of the quantities that should be subject to convergence studies with respect to the \kk-gridbefore embarking on DFPT/EPH calculations.To compute the DOS with the tetrahedron method, use the ebands_edos command of abitk
abitk ebands_edos MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_GSR.ncto read the eigenvalues in the IBZ from a GSR.n file (WFK.nc files are supported as well).This command produces a text file named out_GSR.nc_EDOS that can be visualized with:
abiopen.py out_GSR.nc_EDOS -e
Note that the phonon linewidths have a geometrical contribution due to Fermi surface since \gamma_\qnuis expected to be large in correspondence of \qq wave vectors connecting two portions of the FS.Strictly speaking this is true only if the e-ph matrix elements are constant.In real materials the amplitude of g_{mn\nu}(\kk, \qq) is not constantand this may enhance/suppress the value \gamma_\qnu for particular modes.Yet visualizing the FS is rather useful when discussing e-ph properties in metals.For this reason, it is useful to have a look at the FS with an external
Graphical tools for the visualization of the FS usually require an external file withelectronic energies in the full BZ whereas ab-initio codes usually take advantage of symmetriesto compute \ee_\nk in the IBZ only.To produce a BXSF that can be used to visualize the FS with e.g. XCrysDen ,use the ebands_bxsf command of the abitk utility located in src/98_mainand provide a GSR.nc (WFK.nc) file with energies computed on a \kk-mesh in the IBZ:
abitk ebands_bxsf MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_GSR.ncThis command reconstructs the KS energy in the full BZ by symmetry and producesthe out_GSR.nc_BXSF file that can be opened with XCrysDen using:
xcrysden --bxsf out_GSR.nc_BXSF
Other tools such as fermisurfer or pyprocar can read data in the BXSF format as well.For a minimalistic matplotlib-based approach, use can use the abiview.py script with the fs command:
abiview.py fs MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_GSR.ncto visualize the FS in the unit cell of the reciprocal lattice.
Tip
The BXSF file can be produced at the end of the GS calculationor inside the EPH code by setting prtfsurf to 1 but abitk is quite handy if you alreadyhave a file and you don’t want to write a full Abinit input file and rerun the calculationfrom scratch.
We now focus on the vibrational properties of MgB_2.In principle, one can compute phonon frequencies either with anaddb or with the EPH code.However, for many applications, it is much easierto automate the entire process by invoking anaddb directly from the abiview.py script.To compute the phonon band structure using the DDB file produced on the 4x4x4 \qq-mesh, use:
abiview.py ddb MgB2_eph4isotc/flow_mgb2_phonons/w1/outdata/out_DDBthat produces the following figures:
The Projected Phonon DOS (right panel) shows that the low-energy vibrations (below 40 meV) are mainly of Mg characterwhereas the higher energy modes mainly involve motions of B atoms.Of particular interest for us is the two-fold degenerate band located at ~ 60 meV along the \Gamma-A segment.These are the so-called E_{2g} vibrational modes that couple strongly with the two \sigma bands crossing \ee_F.
Tip
A more detailed analysis reveals that E_{2g} involve in-plane stretching of the B-B bonds.To visualize the phonon modes, one can use the phononwebsite web app by H. Miranda.
abiview.py ddb MgB2_eph4isotc/flow_mgb2_phonons/w1/outdata/out_DDBOpening URL: http://henriquemiranda.github.io/phononwebsite/phonon.html?json=http://localhost:8000/__abinb_workdir__/Mg1B2136lmnep.jsonUsing default browser, if the webpage is not displayed correctlytry to change browser either via command line options or directly in the shell with e.g: export BROWSER=firefoxPress Ctrl+C to terminate the HTTP serverIf the page does not work, upload the json file directly.The vibrational spectrum just obtained looks reasonable: no vibrational instability is observed andthe acoustic modes tend to zero linearly for |\qq| \rightarrow 0(note that the acoustic sum-rule is enforced by default via asr = 1).Yet this does not mean that our results are fully converged.This is clearly seen if we compare the phonon bands computed with a 4x4x4 and a 6x6x6 ab-initio \qq-mesh.Also in this case, we can automate the process by using the ddb command of the abicomp.py script that takesin input an arbitrary list of DDB files, calls anaddb for all the DDB files and finally compare the results:
abicomp.py ddb \ MgB2_eph4isotc/flow_mgb2_phonons/w1/outdata/out_DDB \ MgB2_eph4isotc/666q_DDB -e
The figure reveals that the vibrational spectrum interpolated from the 4x4x4 \qq-meshunderestimates the maximum phonon frequency.Other differences are visible above 40 meV, especially in the E_{2g} modes.
Note
The 666q_DDB file was produced with the same AbiPy script by just changing the valueof the ngqpt variable.The reason why we do not provide files obtained with a 6x6x6 \qq-sampling is that the size of thegit repository including all the DFPT potentials is around 200 Mb.For this reason, we continue using the 4x4x4 DDB file but we should takethis into account when comparing our results with previous works.
Our first computation of the isotropic Tc¶
For our first example, we use a relatively simple input file that allows us to introducethe most important variables and the organization of the results.Copy teph4isotc_2.abi in the working directory and run the code using:
abinit teph4isotc_2.abi > log 2> errTip
If you prefer, you can run it in parallel with e.g two MPI processes using:
mpirun -n 2 abinit teph4isotc_2.abi > log 2> errwithout having to introduce any input variable for the MPI parallelizationas the EPH code can automatically distribute the workload over k-points and spins.Further details concerning the MPI version are given in thelast section of the tutorial
We now discuss the meaning of the different variables in more detail.
To activate the computation of \gamma_{\qq\nu} in metals, we use optdriver = 7 and eph_task = 1.The location of the external DDB, DVDB and WFK files is specified viagetddb_filepath getdvdb_filepath getwfk_filepath, respectively.
getddb_filepath "MgB2_eph4isotc/flow_mgb2_phonons/w1/outdata/out_DDB"getwfk_filepath "MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_WFK.nc"getdvdb_filepath "teph4isotc_1_DVDB"The DDB and the WFK files are taken from the git repository while for the DVDB we usethe file produced by mrgdv in the previous section (remember what we said aboutthe use of relative paths in the input files).Note also the use of the new input variable structure (added in Abinit v9)with the abifile prefix
structure "abifile:MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_DEN.nc"to read the crystalline structure from an external fileso that we do avoid repeating the unit cell in every input file.
Next, we have the variables defining the coarse and the fine \qq-mesh(ddb_ngqpt and eph_ngqpt_fine, respectively):
ddb_ngqpt 4 4 4 # The ab-initio q-mesh (DDB, DVDB)eph_ngqpt_fine 6 6 6 # Activate interpolation of DFPT potentials # gamma_{q,nu} are computed of this fine grid.dipdip 0 # No treatment of the dipole-dipole part. OK for metalsdipdip is set to zero as we are dealing with a metal and the inter-atomic forceconstants are short-ranged provided we ignore possible Kohn-anomalies [Kohn1959].
The \kk-point integration is performed with the Gaussian smearing(eph_intmeth == 1 and 0.1 eV for eph_fsmear).
eph_intmeth 1 # Gaussian method for double-delta integration.eph_fsmear 0.1 eV # Constant broadening in Gaussian function.eph_fsewin 0.3 eV # Energy window for wavefunctions.Note
In order to accelerate the calculation, we have decreased eph_fsewinfrom its default value of 1.0 eV to 0.3 eV.This is possible when the Gaussian method is used sincestates whose energy is 3-4 standard deviation from the Fermi level givenegligible contribution to the double delta integral.In other words, for the Gaussian technique, an optimal value of eph_fsewincan be deduced from eph_fsmearUnfortunately, this is not possible when the tetrahedron method is used and, to some extent,also when the adaptive broadening is used.Hopefully, the next versions of the code will provide …0.04 Ha ~ 1 eV
Then we activate two tricks to make the calculation faster:
mixprec 1boxcutmin 1.1As explained in the documentation, using mixprec = 1 and 2 > boxcutmin >= 1.1should not have significant effects on the final results yet these are not the default valuesas users are supposed to compare the results with/without these tricks (especially boxcutmin)before running production calculations.
Output files¶
We now discuss in more detail the main output file produced by the EPH run.
After the standard section with info on the unit cell and the pseudopotentials,we find info on the electronic DOS:
Linear tetrahedron method. Mesh step: 10.0 (meV) with npts: 2757 From emin: -5.1 to emax: 22.5 (eV) Number of k-points in the IBZ: 133 Fermi level: 7.63392108E+00 (eV) Total electron DOS at Fermi level in states/eV: 7.36154553E-01 Total number of electrons at eF: 8.0- Writing electron DOS to file: teph4isotc_2o_DS1_EDOSThen the code outputs some basic quantities concerning the Fermi surfaceand the method for the BZ integration of the double delta:
==== Fermi surface info ==== FS integration done with adaptive gaussian method Total number of k-points in the full mesh: 1728 For spin: 1 Number of BZ k-points close to the Fermi surface: 291 [ 16.8 %] Maximum number of bands crossing the Fermi level: 3 min band: 3 Max band: 5Then, for each \qq-point in the IBZ, the code outputsthe values of \ww_\qnu, \gamma_\qnu and \lambda_\qnu:
q-point = 0.000000E+00 0.000000E+00 0.000000E+00 Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n) 1 0.000000E+00 0.000000E+00 0.000000E+00 2 0.000000E+00 0.000000E+00 0.000000E+00 3 0.000000E+00 0.000000E+00 0.000000E+00 4 1.444402E-03 4.898804E-07 3.731164E-03 5 1.444402E-03 4.898804E-07 3.731164E-03 6 1.674636E-03 2.955247E-15 1.674492E-11 7 2.426616E-03 1.025439E-03 2.767185E+00 8 2.426616E-03 1.025439E-03 2.767185E+00 9 3.127043E-03 3.508418E-08 5.701304E-05The frequencies of the acoustic modes at \Gamma are zero (as they should be)since asr is automatically set to 1 so the acoustic rule is automatically enforced.Also, the linewidths of the acoustic modes are zero.In this case, however, the code uses a rather simple heuristic rule:
abiview.py ddb_asr MgB2_eph4isotc/flow_mgb2_phonons/w1/outdata/out_DDBFinally, we have the value of the isotropic \lambda:
lambda= 0.4286The calculation has produced the following output files:
$ ls teph4isotc_2o_DS1*teph4isotc_2o_DS1_A2F.nc teph4isotc_2o_DS1_NOINTP_A2FW teph4isotc_2o_DS1_PHDOS.ncteph4isotc_2o_DS1_A2FW teph4isotc_2o_DS1_NOINTP_PH_A2FW teph4isotc_2o_DS1_PHGAMMAteph4isotc_2o_DS1_EBANDS.agr teph4isotc_2o_DS1_PHBANDS.agr teph4isotc_2o_DS1_PH_A2FWteph4isotc_2o_DS1_EDOS teph4isotc_2o_DS1_PHBST.ncThe files without the .nc extension are text files that can be plotted with e.g. gnuplot or xmgrace .More specifically,
- Apple
- Pomaceous fruit of plants of the genus Malus in the family Rosaceae.
The A2F.nc netcdf file stores all the results of the calculation in a format that cancan be visualized with abiopen.py :
abiopen.py teph4isotc_2o_DS1_A2F.nc -ethat produces the following plot:
Our first results should be interpreted with a critical eye due to coarse sampling used,yet our \alpha^2F(\ww) already shows features that are observed in other (more accurate) calculationsreported in the literature: a broad peak at ~70 meV that gives the most important contribution to the total \lambdaand a second smaller peak at ~70 meV.The sharp delta-like peaks are an artifact of the linear interpolation and will hopefully disappearif denser \qq-meshes are used.Our initial values for \lambda and \ww_{log} are:
a2f(w) on the [4 4 4] q-mesh (ddb_ngqpt|eph_ngqpt)Isotropic lambda: 0.39, omega_log: 0.067 (eV), 780.042 (K)kIf we compare our results for \lambda with those reported in [Margine2013],
we see that our first calculation underestimates \lambda by almost a factor two.Note, however, that all the previous studies used much denser samplings and all of themuse at least a 6x6x6 \qq-mesh for phonons.It is clear that we need to densify the BZ sampling to get more reliable results yetour results are not that bad considering that the calculation took less than XXX minutes in sequential.
Using the tetrahedron method¶
In this section, we repeat the calculation done in teph4isotc_2.abibut now with the optimized tetrahedron scheme by [Kawamura2014].This test will show that (i) phonon linewidths are very sensitive to the \kk-meshand the integration scheme and that (ii) one has to test different values of eph_fsewinto find the optimizal one.
The reason is that the optimized tetrahedron scheme constructs a third-order interpolation functionwith 20 \kk-points so states that are relatively far from the Fermi level contributes to the integrand.If eph_fsewin is too small, part of the weigt is lost.To find an optimal value, we can use:
ndtset 4eph_intmeth 2 # Tetraeph_fsewin: 0.4 eVeph_fsewin+ 0.4 eVWe do not provide an input file to perform such test as one can easily change teph4isotc_2.abi.and the run the calculation in parallel.To compare the results, we use the a2f command of the abicomp.py script and pass the list of A2F.c files:
abicomp.py a2f teph4isotc_2o_DS*_A2F.nc -e
Preparing the convergence study wrt the k-mesh¶
The NSCF computation of the WFK becomes quite CPU-consuming and memory-demanding if dense \kk-meshes are needed.Fortunately, we can optimize this part since the computation of \gamma_\qnurequires the knowledge of Bloch states inside a relatively small energy window around \ee_F.Similarly to what is done in the eph4mob tutorial, we can therefore take advantage ofthe star-function SKW interpolation to find the \kk wavevectors whose energy is insidethe sigma_erange energy window around the Fermi level.The choice of an optimal window is discussed afterwards.
First of all, we strongly recommend to test whether the SKW interpolation can reproducethe ab-initio results with reasonable accuracyby comparing the ab-initio band structure with the interpolated one byusing abitk with the skw_compare command:
abitk skw_compare \ MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_GSR.nc \ MgB2_eph4isotc/flow_mgb2_phonons/w0/t1/outdata/out_GSR.nc --is-metalThe fist argument is a GSR.nc (WFK.nc) file with the energies in the IBZthat will be used to build the star-function interpolant.The quality of the interpolant obviously depends on the density of this \kk-mesh.The second file contains the ab-initio eigenvalues along a \kk-path computed with a standard NSCF run.In this example, we are using precomputed GSR.nc files with a 12x12x12 \kk-mesh for the SKW interpolation.
Note
Note the use of the --is-metal option else abitk will complain that it cannot compute the gapsas the code by default assumes a semiconductor.
The skw_compare command has produced two netcdf files (abinitio_EBANDS.nc and skw_EBANDS.nc)that we can be compared with the abicomp.py script and the ebands command:
abicomp.py ebands abinitio_EBANDS.nc skw_EBANDS.nc -p combiplotto produce the following plot:
The figure shows that the SKW interpolant nicely reproduces the ab-initio dispersion around \ee_F.Discrepancies between the ab-initio results and the interpolated values are visiblein correspondence of band-crossings.This is expected since SKW is a Fourier-based approach and band-crossing leads to a non-analytic behaviourof the signal that cannot be reproduced with a finite number of Fourier components.Fortunately, these band crossings are relatively far from the Fermi level hence they do not enter into playin the computation of superconducting properties.
Tip
Should the fit be problematic, use the command line options:
abitk skw_compare IBZ_WFK KPATH_WFK [--lpratio 5] [--rcut 0.0] [--rsigma 0]to improve the interpolation and take note of the SKW parametersas the same values should be used when calling ABINIT with the einterp input variable.The worst case scenario of band crossings and oscillations close to the Fermi levelcan be handled by enlarging the sigma_erange energy window.
At this point, we are confident that the SKW interpolation is OK and we can use it to locatethe \kk-wavevectors of a much denser \kk-mesh whose energy is inside the sigma_erange windowaround the Fermi level.This is done in the teph4isotc_4.abi input file:
The most important section of the input file is reproduced below:
optdriver 8wfk_task "wfk_kpts_erange"getwfk_filepath "MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_WFK.nc"# Define fine k-mesh for the SKW interpolationsigma_ngkpt 24 24 24sigma_nshiftk 1sigma_shiftk 0 0 0sigma_erange -0.3 -0.3 eV # Select kpts in the fine mesh within this energy window.einterp 1 5 0 0 # Parameters for star-function interpolation (default values)The first part (optdriver and wfk_task) activates the computation of the KERANGE.nc filewith ab-initio energies in the IBZ taken from getwfk_filepath.The three variables sigma_ngkpt, sigma_nshiftk and sigma_shiftk define the final dense \kk-mesh.Finally, sigma_erange defines the energy window around the Fermi levelwhile einterp lists the parameters passed to the SKW routinesIn this example, we are using the default values.You may need to change einterp if you had to use different options for the SKW interpolationwhen using abitk skw_compare.
Important
When dealing with metals, both entries in sigma_erange must be negativeso that the energy window is refered to \ee_F and not to the CBM/VBM used in semiconductors.Remember to use a value for the window that is reasonably large in order to accountfor possible oscillations and/or inaccuracies of the SKW interpolation around \ee_F.
Once we have the KERANGE.nc file, we can use it to perform a NSCF calculation to generatea customized WFK file on the dense \kk-mesh:This is done in the teph4isotc_4.abi input file:
to perform a NSCF calculation with kptopt 0 to produce a new WFK file on the dense \kk-mesh.
iscf -2tolwfr 1e-18kptopt 0 # Important!with the getkerange_filepath:
Note
Note how we use getden_filepath and the syntax:
getkerange_filepath "teph4isotc_3o_DS1_KERANGE.nc"# Read DEN file to initialize the NSCF run.getden_filepath "MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_DEN.nc"# Init GS wavefunctions from this file (optional).getwfk_filepath "MgB2_eph4isotc/flow_mgb2_phonons/w0/t0/outdata/out_WFK.nc"to start the NSCF run from the precomputed GS DEN file, initialize the trial wafefunctionsfrom getwfk_filepath to accelerate a bit the calculation.
TODO: Discussion about energy range and integration scheme.
Convergence study wrt to the k/q-mesh¶
At this point, we can use the WFK file to perform EPH calculations with denser \kk-meshes.We will be using settings similar to the ones used in teph4isotc_2.abi except forthe use of eph_ngqpt_fine, getwfk_filepath and ngkpt:
Notes on the MPI parallelism¶
EPH calculations performed with eph_task = 1 support 4 different levels of MPI parallelism andthe number of MPI processes for each level can be specified via the eph_np_pqbks input variable.Note that we wrote 4 MPI levels instead of 5 because, for isotropic T_c calculations,the band parallelism is not supported, only atomic perturbations, \qq-points, \kk-points and spins (if any).
The eph_np_pqbks variable is optional in the sense that whatever number of MPI processes you use,the EPH code will be able to parallelize the calculation by distributing over \kk-points and spins (if any).This distribution, however, may not be optimal, especially if you have lot of CPUs availableand/or you need to decrease the memory requirements per CPU.
To decrease memory, we suggest to activate the parallelism over perturbations(up to 3 x natom) to make the memory for W(\rr, \RR, \text{3 x natom}) scale.If memory is not of concern, activate the \qq-point parallelism for better efficiency.
Important
Note that the code is not able to distribute the memory for the wavefunctions althoughonly the states in the eph_fsewin energy window are read and stored in memory.
