1\documentclass[a4paper,11pt,twoside]{article} 2\usepackage[top=3cm,bottom=3cm,left=2cm,right=2cm]{geometry} 3%\usepackage{multirow} 4\usepackage{footnote} 5\usepackage{amsbsy} 6\usepackage{amsmath} 7\usepackage{graphicx} 8\usepackage{fancyhdr} 9\usepackage{subfig} 10\usepackage[T1]{fontenc} % important for having searchable underscores 11\usepackage{bm}% bold math 12\usepackage[sort&compress,numbers]{natbib} 13%\setlength{\parindent}{0in} 14%\setlength{\parskip}{0.05in} 15\setlength{\parskip}{0.1in} 16 17%%%%%%%%%%%%%%%%%%%%%%%%%%% REMOVE AT THE END %%%%%%%%%%%%%%%%%%%%%%%%%%% 18\usepackage[usenames]{color} % colors 19\usepackage{ulem} % allows \sout (strikeout). Can remove at the end (Ivo) 20%\usepackage{soul} % allows \hl (highlight). Can remove at the end (Ivo) 21\newcommand\tent[1]{\textcolor{red}{#1}} % for tentative changes 22\newcommand\tentbis{\color{red}} % for tentative changes 23%%%%%%%%%%%%%%%%%%%%%%%%%%% REMOVE AT THE END %%%%%%%%%%%%%%%%%%%%%%%%%%% 24%\let\oldundercore=\_ 25%\def\_{\oldunderscore} 26 27\usepackage{textcomp} 28 29% set fancy headings 30 31\pagestyle{fancy} 32\lhead[{\it \thepage}]{{\bf\it {\tt wannier90}: Tutorial}} 33\chead{} 34\rhead[{\bf\it {\tt wannier90}: Tutorial}]{{\it \thepage}} 35\renewcommand{\headrulewidth}{0.2pt} 36\lfoot{} 37\cfoot{} 38\rfoot{} 39\renewcommand{\footrulewidth}{0pt} 40\setlength{\footskip}{0.25in} 41\setlength{\parindent}{0in} 42 43\title{\wannier: Tutorial} 44 45\author{Version 3.1} 46 47\date{} 48%\date{\today} 49 50%%% THIS SHOULD BE THE LAST PACKAGE TO BE LOADED! 51%hidelinks to remove colored borders from links 52%plainpages=false needed for the unnumbered pages 53%breaklinks=true to allow to break long links (e.g. long titles) on 54%more than one line 55%bookmarksopen open the bookmarks in the adobe reader plugin 56%bookmarksopenlevel decide the max level at which the bookmarks should 57%be open 58% pdfdisplaydoctitle=true to show the document title instead of the 59% filename in the titlebar of adobereader 60\usepackage[plainpages=false,breaklinks=true,pdfborder=0 0 0,pdfdisplaydoctitle=true,bookmarksopen=true,bookmarksopenlevel=0,pdftex,% 61% Comment the following line if it makes problems (it requires a recent 62% LaTeX distribution) 63 hidelinks, 64 pdftitle={Wannier90 Tutorial}, 65 pdfkeywords={wannier90;postw90;pwscf;tutorial;examples}]{hyperref} 66 67 68\begin{document} 69\newcommand{\wannier}{{\rm\texttt{wannier90}}} 70\newcommand{\postw}{{\rm\texttt{postw90}}} 71\newcommand{\bw}{{\rm\texttt{BoltzWann}}} 72\newcommand{\pwscf}{\textsc{pwscf}} 73\newcommand{\QE}{\textsc{quantum-espresso}} 74\newcommand{\Mkb}{\mathbf{M}^{(\mathbf{k},\mathbf{b})}} 75\newcommand{\Ak}{\mathbf{A}^{(\mathbf{k})}} 76\newcommand{\Uk}{\mathbf{U}^{(\mathbf{k})}} 77 78\newcommand\sectiontitle[1]{\section*{#1}\addcontentsline{toc}{section}{#1}} 79 80\maketitle 81 82\tableofcontents 83 84\newpage 85 86 87\sectiontitle{Preliminaries} 88\label{sec:preliminaries} 89 90Welcome to \wannier! The examples contained in this tutorial are 91designed to help you become familiar with the procedure of generating, 92analysing and using maximally-localised Wannier functions (MLWFs). As 93a first step, install \wannier\ following the instructions in the {\tt 94 README} file of the \wannier\ distribution. For an introduction to 95the theory underlying MLWFs, you are encouraged to refer to the brief 96overview given in the \wannier\ User Guide~\cite{UserGuide}, to the 97two seminal papers of Refs.~\cite{marzari-prb97,souza-prb01}, a recent 98review article~\cite{marzari-rmp12} and to a 99paper~\cite{mostofi-cpc08} describing \wannier. 100 101The following additional programs may be installed in order to 102visualise the output of \wannier\ (they are optional, not all of them 103are necessary) 104\begin{itemize} 105\item {\tt gnuplot} is used to plot bandstructures. It is 106available for many operating systems and is often installed by default on 107 Unix/Linux distributions\\ 108\url{http://www.gnuplot.info} 109\item {\tt xmgrace} may also be used to plot bandstructures.\\ 110\url{http://plasma-gate.weizmann.ac.il/Grace} 111\item {\tt XCrySDen} is used to visualise crystal structures, MLWFs, 112 and Fermi surfaces. It is available for Unix/Linux, 113 Windows (using cygwin), and OSX. To correctly display 114files from \wannier, version 1.4 or later must be used.\\ 115\url{http://www.xcrysden.org} 116\item {\tt vmd} can also be used to visualise crystal structures and 117 MLWFs.\\ 118\url{http://www.ks.uiuc.edu/Research/vmd} 119\item{\tt python} with the {\tt numpy} and {\tt matplotlib} modules 120 is used in examples 17--19\\ 121 \url{http://www.python.org}\\ 122 \url{http://www.numpy.org}\\ 123 \url{http://matplotlib.org} 124\end{itemize} 125 126\sectiontitle{Parallel execution} 127\label{sec:parallel} 128{\tt postw90.x} and {\tt wannier90.x} can be run in parallel to speed up 129the calculations, using the MPI libraries. 130 131To enable the parallel version to be built, you must specify some 132flags in the {\tt make.inc} file of {\tt wannier90} and {\tt postw90}; 133for further information, please refer to the {\tt README.install} file 134in the top directory of the {\tt wannier90} distribution. 135 136Then, to run e.g. with 8 processors, you typically need to run a 137command similar to {\tt postw90} as follows: 138\begin{verbatim} 139mpirun -np 8 postw90.x seedname 140\end{verbatim} 141(the {\tt mpirun} command and its flags may differ depending on the 142MPI libraries installed on your system: refer to your MPI manual and/or to 143your system administrator for further information). 144 145 146 147\sectiontitle{About this tutorial} 148 149The first part of this tutorial comprises four examples taken from 150Refs.~\cite{marzari-prb97,souza-prb01}: gallium arsenide, lead, 151silicon and copper. All of the \wannier\ input files have been 152provided. 153 154The second part of the tutorial covers the generation of \wannier\ 155input files starting from a full electronic structure calculation. We 156have provided input files for the \pwscf\ interface (\url{http://www.quantum-espresso.org}) to \wannier. 157Therefore, you 158will need to install and compile elements of the {\tt 159 quantum-espresso} package, namely {\tt pw.x} and {\tt 160 pw2wannier90.x}, in order to run these 161examples. Please visit \url{http://www.quantum-espresso.org} to download the 162package, and for installation instructions. 163The tutorial examples work with \pwscf\ v5.1.x and v6.0.x. The exception are the examples on symmetry adapted Wannier functions which require v6.0.x together with the 164very latest version of {\tt pw2wannier90.f90}. This can be found in the directory {\tt pwscf/v6.0} in the wannier distribution. It should be moved to {\tt PP/src} in the \pwscf\ distribution and compiled using {\tt make pp}. Later versions v6.x.x should have the most up-to-date version 165 of pw2wannier90.f90 already included in the Quantum ESPRESSO 166 distribution. 167 168 169There are interfaces to a number of other electronic structure codes 170including {\sc abinit} (\url{http://www.abinit.org}), {\sc fleur} 171(\url{http://www.flapw.de}), {\sc OpenMX} (\url{http://www.openmx-square.org/}), 172{\sc GPAW} (\url{https://wiki.fysik.dtu.dk/gpaw/}), {\sc VASP} (\url{http://www.vasp.at}), and 173{\sc Wien2k} (\url{http://www.wien2k.at}) 174 175\sectiontitle{Contact us} 176 177If you have any suggestions regarding ways in which this tutorial may 178be improved, then send us an email. 179% at {\tt developers@wannier.org}. 180 181For other questions, email the \wannier\ forum at {\tt 182 wannier@quantum-espresso.org}. Note that first you will need to 183register in order to post emails. Emails from non-registered users are 184deleted automatically. You can register by following the links at\\ 185\url{http://www.wannier.org/forum.html}. 186 187 188 189\cleardoublepage 190 191\sectiontitle{1: Gallium Arsenide -- MLWFs for the valence bands} 192 193\begin{itemize} 194\item{Outline: \it{Obtain and plot MLWFs for the four valence 195 bands of GaAs.}} 196\item{Generation details: \it{From \pwscf, using norm-conserving 197 pseudopotentials and a 2$\times$2$\times$2 k-point grid. Starting 198 guess: four bond-centred Gaussians.}} 199\item{Directory: {\tt examples/example1/}} 200\item{Input Files} 201\begin{itemize} 202\item{ {\tt gaas.win} {\it The master input file}} 203\item{ {\tt gaas.mmn} {\it The overlap matrices $\Mkb$}} 204\item{ {\tt gaas.amn} {\it Projection $\Ak$ of the Bloch states onto a set 205 of trial localised orbitals}} 206\item{ {\tt UNK00001.1} {\it The Bloch states in the real space unit 207 cell. For plotting only.}} 208\end{itemize} 209\end{itemize} 210 211\begin{enumerate} 212\item Run \wannier\ to minimise the MLWFs spread 213{\tt 214\begin{quote} 215wannier90.x gaas 216\end{quote} } 217Inspect the output file {\tt gaas.wout}. The total spread converges to its 218minimum value after just a few iterations. Note that the geometric centre of 219each MLWF lies along a Ga-As bond, slightly closer to As 220than Ga. Note also that the memory requirement for the minimisation of 221the spread is very low as the MLWFs are defined at each 222k-point by just the 4$\times$4 unitary matrices $\Uk$. 223\item Plot the MLWFs by adding the following keywords to 224 the input file {\tt gaas.win} 225{\tt 226\begin{quote} 227wannier\_\-plot = true 228\end{quote} } and re-running \wannier. To visualise the MLWFs we must 229represent them explicitly on a real space grid (see 230Ref.~\cite{UserGuide}). As a consequence, plotting the MLWFs is slower 231and uses more memory than the minimisation of the spread. The four 232files that are created ({\tt gaas\_00001.xsf}, etc.) can be viewed 233using {\tt XCrySDen},\footnote{Once {\tt XCrySDen} starts, click on 234 {\tt Tools} $\rightarrow$ {\tt Data Grid} in order to specify an 235 isosurface value to plot.} e.g., {\tt 236\begin{quote} 237xcrysden \texttt{-{}-}xsf gaas\_00001.xsf 238\end{quote} } 239 240For large systems, plotting the MLWFs may be time consuming 241and require a lot of memory. Use the keyword {\tt wannier\_plot\_list} 242to plot a subset of the MLWFs. E.g., to plot the 2431st and 3rd MLWFs use 244{\tt 245\begin{quote} 246wannier\_plot\_list = 1 3 247\end{quote} } 248The MLWFs are plotted in a supercell of the unit cell. The 249size of this supercell is set through the keyword {\tt 250 wannier\_plot\_supercell}. The default value is 2 (corresponding to a 251supercell with eight times the unit cell volume). We recommend not using 252values great than 3 as the memory and computational cost scales 253cubically with supercell size. 254 255Plot the 3rd MLWFs in a supercell of size 3. Choose a low 256value for the isosurface (say 0.5). Can you explain what you see? 257 258{\it Hint:} For a finite k-point mesh, the MLWFs are in fact 259periodic and the period is related to the spacing of the k-point mesh. For 260mesh with $n$ divisions in the $i^{\mathrm{th}}$ direction in the 261Brillouin zone, the MLWFs ``live'' in a supercell $n$ times the 262unit cell. 263\end{enumerate} 264 265 266%\cleardoublepage 267 268 269 270\sectiontitle{2: Lead -- Wannier-interpolated Fermi surface} 271 272\begin{itemize} 273\item{Outline: \it{Obtain MLWFs for the four lowest states 274 in lead. Use Wannier interpolation to plot the Fermi surface.}} 275\item{Generation Details: \it{From \pwscf, using norm-conserving 276 pseudopotentials and a 4$\times$4$\times$4 k-point grid. Starting 277 guess: atom-centred sp$^3$ hybrid orbitals}} 278\item{Directory: {\tt examples/example2/}} 279\item{Input Files} 280\begin{itemize} 281\item{ {\tt lead.win} {\it The master input file}} 282\item{ {\tt lead.mmn} {\it The overlap matrices $\Mkb$}} 283\item{ {\tt lead.amn} {\it Projection $\Ak$ of the Bloch states onto a set 284 of trial localised orbitals}} 285\item{ {\tt lead.eig} {\it The Bloch eigenvalues at each k-point. For 286 interpolation only}} 287\end{itemize} 288 289\end{itemize} 290 291The four lowest valence bands in lead are separated in energy from the 292higher conduction states (see Fig.~\ref{fig:pb-bnd}). The MLWFs of 293these states have partial occupancy. MLWFs describing only the occupied 294states would be poorly localised. 295 296\begin{enumerate} 297\item Run \wannier\ to minimise the MLWFs spread 298{\tt 299\begin{quote} 300wannier90.x lead 301\end{quote} } 302Inspect the output file {\tt lead.wout}. 303\item Use Wannier interpolation to generate the Fermi surface of 304 lead. Rather than re-running the whole calculation we can use the 305 unitary transformations obtained in the first calculation and 306 restart from the plotting routine. Add the following keywords to the 307 {\tt lead.win} file: {\tt 308\begin{quote} 309restart = plot 310 311fermi\_energy = 5.2676 312 313fermi\_surface\_plot = true 314\end{quote} } and re-run \wannier. The value of the Fermi energy 315(5.2676\,eV) was obtained from the initial first principles 316calculation. \wannier\ calculates the band energies, through 317 318interpolation, on a dense mesh of k-points in the Brillouin zone. The 319density of this grid is controlled by the keyword {\tt 320 fermi\_surface\_num\_points}. The default value is 50 (i.e., 50$^3$ 321points). The Fermi surface file {\tt lead.bxsf} can be viewed using 322{\tt XCrySDen}, e.g., 323% 324{\tt 325\begin{quote} 326xcrysden \texttt{-{}-}bxsf lead.bxsf 327\end{quote} } 328\end{enumerate} 329 330\begin{figure}[h] 331\begin{center} 332\includegraphics{lead} 333\caption{Bandstructure of lead showing the position of the Fermi 334 level. Only the lowest four bands are included in the calculation.} 335\label{fig:pb-bnd} 336\end{center} 337\end{figure} 338 339%\cleardoublepage 340 341 342\sectiontitle{3: Silicon -- Disentangled MLWFs} 343 344\begin{itemize} 345\item{Outline: \it{Obtain disentangled MLWFs for the valence and 346 low-lying conduction states of Si. Plot the interpolated 347 bandstructure}} 348\item{Generation Details: \it{From \pwscf, using norm-conserving 349 pseudopotentials and a 4$\times$4$\times$4 k-point grid. Starting 350 guess: atom-centred sp$^3$ hybrid orbitals}} 351\item{Directory: {\tt examples/example3/}} 352\item{Input Files} 353\begin{itemize} 354\item{ {\tt silicon.win} {\it The master input file}} 355\item{ {\tt silicon.mmn} {\it The overlap matrices $\Mkb$}} 356\item{ {\tt silicon.amn} {\it Projection $\Ak$ of the Bloch states onto a 357 set of trial localised orbitals}} 358\item{ {\tt silicon.eig} {\it The Bloch eigenvalues at each k-point}} 359\end{itemize} 360\end{itemize} 361The valence and lower conduction states can be represented by MLWFs 362with $sp^3$-like symmetry. The lower conduction states are not 363separated from the higher states by an energy gap. In order to form 364localised WF, we use the disentanglement procedure 365introduced in Ref.~\cite{souza-prb01}. The position of the inner and outer 366energy windows are shown in Fig.~\ref{fig:si.bnd}. 367\begin{enumerate} 368\item Run \wannier. 369{\tt 370\begin{quote} 371wannier90.x silicon 372\end{quote} } 373Inspect the output file {\tt silicon.wout}. The minimisation of the 374spread occurs in a two-step procedure~\cite{souza-prb01}. First, we minimise 375$\Omega_{\rm I}$ -- this is the extraction of the optimal subspace in 376the disentanglement procedure. Then, we minimise $\Omega_{\rm D} + 377\Omega_{{\rm OD}}$. 378 379\item Plot the energy bands by adding the following 380 commands to the input file {\tt silicon.win} {\tt 381\begin{quote} 382restart = plot 383 384bands\_plot = true 385\end{quote} } 386and re-running \wannier. The files {\tt silicon\_band.dat} and {\tt 387 silicon\_band.gnu} are created. 388To plot the bandstructure using gnuplot 389\smallskip 390{\tt 391\begin{quote} 392myshell> gnuplot 393 394gnuplot> load `silicon\_band.gnu' 395\end{quote} } 396The k-point path for the bandstructure interpolation is set in the {\tt 397 kpoint\_path} block. Try plotting along different paths. 398\end{enumerate} 399 400\begin{figure}[h] 401\begin{center} 402\includegraphics{si} 403\caption{Bandstructure of silicon showing the position of the outer 404 and inner energy windows.} 405\label{fig:si.bnd} 406\end{center} 407\end{figure} 408 409%\cleardoublepage 410 411 412\sectiontitle{4: Copper -- Fermi surface, orbital character of energy 413 bands} 414 415\begin{itemize} 416\item{Outline: \it{Obtain MLWFs to describe the states around the 417 Fermi-level in copper}} 418\item{Generation Details: \it{From \pwscf, using ultrasoft 419 pseudopotentials~\cite{vanderbilt-prb90} and a 420 4$\times$4$\times$4 k-point grid. Starting guess: five 421 atom-centred d orbitals, and two s orbitals centred on one of each 422 of the two tetrahedral interstices.}} 423\item{Directory: {\tt examples/example4/}} 424\item{Input Files} 425\begin{itemize} 426\item{ {\tt copper.win} {\it The master input file}} 427\item{ {\tt copper.mmn} {\it The overlap matrices $\Mkb$}} 428\item{ {\tt copper.amn} {\it Projection $\Ak$ of the Bloch states onto a 429 set of trial localised orbitals}} 430\item{ {\tt copper.eig} {\it The Bloch eigenvalues at each k-point}} 431\end{itemize} 432 433\end{itemize} 434 435\begin{enumerate} 436\item Run \wannier\ to minimise the MLWFs spread 437{\tt 438\begin{quote} 439wannier90.x copper 440\end{quote} } 441Inspect the output file {\tt copper.wout}. 442 443\item Plot the Fermi surface, it should look familiar! The Fermi 444 energy is at 12.2103\,eV. 445 446\item Plot the interpolated bandstructure. A suitable path in k-space is 447\smallskip 448{\tt 449\begin{quote} 450begin kpoint\_path 451 452G 0.00 0.00 0.00 X 0.50 0.50 0.00 453 454X 0.50 0.50 0.00 W 0.50 0.75 0.25 455 456W 0.50 0.75 0.25 L 0.00 0.50 0.00 457 458L 0.00 0.50 0.00 G 0.00 0.00 0.00 459 460G 0.00 0.00 0.00 K 0.00 0.50 -0.50 461 462end kpoint\_path 463\end{quote} } 464\item Plot the contribution of the interstitial WF to the 465 bandstructure. Add the following keyword to {\tt copper.win} 466\smallskip 467{\tt 468\begin{quote} 469bands\_plot\_project = 6,7 470\end{quote} } The resulting file {\tt copper\_band\_proj.gnu} can be 471opened with gnuplot. Red lines correspond to a large contribution from 472the interstitial WF (blue line are a small contribution; ie a large 473$d$ contribution). 474 475 476\end{enumerate} 477 478 479 480 481Investigate the effect of the outer and inner energy window on the 482interpolated bands. 483 484 485 486\begin{figure}[h] 487\begin{center} 488\includegraphics{cu} 489\caption{Bandstructure of copper showing the position of the outer 490 and inner energy windows.} 491\label{fig:cu-bnd} 492\end{center} 493\end{figure} 494 495\clearpage 496 497\sectiontitle{Examples Using the {\sc pwscf} Interface} 498\label{sec:using pwscf} 499 500The \pwscf\ plane-wave, density-functional theory code, which is 501available as part of the \QE\ distribution (\url{http://www.quantum-espresso.org}), is fully interfaced to \wannier\ via the 502{\tt pw2wannier90} post-processing code that is also available as part 503of \QE. The latest version of {\tt pw2wannier90} is included as part of 504the \wannier\ distribution. Please see the {\tt pwscf} directory for 505instructions on how to incorporate it into \pwscf. 506 507Note that both the {\tt PWSCF} executable {\tt pw.x} {\it and} {\tt 508 pw2wannier90.x} can be run in parallel, which for large calculations 509can reduce the computation time very significantly. This requires 510compiling the code in its parallel version, using the MPI 511libraries. Refer to the \QE\ package for the documentation on how to 512do so. Note that, unless you specify {\tt wf\_collect=.true.} in your 513{\tt pw.x} input file, you must run {\tt pw2wannier90} with the same 514number of processors as {\tt pw.x}. 515 516Moreover we remind here that both the \wannier\ executable and {\tt postw90.x} can be run in 517parallel. In this case any number of processors can be used, 518independently of the number used for {\tt pw.x} and {\tt 519 pw2wannier90.x}. 520 521%\cleardoublepage 522 523\sectiontitle{5: Diamond -- MLWFs for the valence bands} 524\begin{itemize} 525\item{Outline: \it{Obtain MLWFs for the valence bands of diamond}} 526\item{Directory: {\tt examples/example5/}} 527\item{Input Files} 528\begin{itemize} 529\item{ {\tt diamond.scf} {\it The \pwscf\ input file for ground state 530 calculation}} 531\item{ {\tt diamond.nscf} {\it The \pwscf\ input file to obtain Bloch 532 states on a uniform grid}} 533\item{ {\tt diamond.pw2wan} {\it The input file for {\tt pw2wannier90}}} 534\item{ {\tt diamond.win} {\it The {\tt wannier90} input file}} 535\end{itemize} 536\end{itemize} 537 538\begin{enumerate} 539\item Run \pwscf\ to obtain the ground state of diamond\\ 540{\tt pw.x < diamond.scf > scf.out} 541 542\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 543{\tt pw.x < diamond.nscf > nscf.out} 544 545\item Run \wannier\ to generate a list of the required overlaps (written 546 into the {\tt diamond.nnkp} file).\\ 547{\tt wannier90.x -pp diamond} 548 549\item Run {\tt pw2wannier90} to compute the overlap between Bloch 550 states and the projections for the starting guess (written in the 551 {\tt diamond.mmn} and {\tt diamond.amn} files).\\ 552{\tt pw2wannier90.x < diamond.pw2wan > pw2wan.out} 553 554\item Run \wannier\ to compute the MLWFs.\\ 555{\tt wannier90.x diamond} 556\end{enumerate} 557 558%\cleardoublepage 559 560\sectiontitle{6: Copper -- Fermi surface} 561 562\begin{itemize} 563\item{Outline: \it{Obtain MLWFs to describe the states around the 564 Fermi-level in copper}} 565\item{Directory: {\tt examples/example6/}} 566\item{Input Files} 567\begin{itemize} 568\item{ {\tt copper.scf} {\it The \pwscf\ input file for ground state 569 calculation}} 570\item{ {\tt copper.nscf} {\it The \pwscf\ input file to obtain Bloch states 571 on a uniform grid}} 572\item{ {\tt copper.pw2wan} {\it Input file for {\tt pw2wannier90}}} 573\item{ {\tt copper.win} {\it The {\tt wannier90} input file}} 574\end{itemize} 575 576\end{itemize} 577 578\begin{enumerate} 579\item Run \pwscf\ to obtain the ground state of copper\\ 580{\tt pw.x < copper.scf > scf.out} 581 582\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 583{\tt pw.x < copper.nscf > nscf.out} 584 585\item Run \wannier\ to generate a list of the required overlaps (written 586 into the {\tt copper.nnkp} file).\\ 587{\tt wannier90.x -pp copper} 588 589\item Run {\tt pw2wannier90} to compute the overlap between Bloch 590 states and the projections for the starting guess (written in the 591 {\tt copper.mmn} and {\tt copper.amn} files).\\ 592{\tt pw2wannier90.x < copper.pw2wan > pw2wan.out} 593 594\item Run \wannier\ to compute the MLWFs.\\ 595{\tt wannier90.x copper} 596\end{enumerate} 597 598Inspect the output file {\tt copper.wout}. 599 600\begin{enumerate} 601\item Use Wannier interpolation to obtain the Fermi surface of 602 copper. Rather than re-running the whole calculation we can use the 603 unitary transformations obtained in the first calculation and restart 604 from the plotting routine. Add the following keywords to the {\tt 605 copper.win} file: 606{\tt 607\begin{quote} 608restart = plot 609 610fermi\_energy = [insert your value here] 611 612fermi\_surface\_plot = true 613\end{quote} } and re-run \wannier. The value of the Fermi energy can 614be obtained from the initial first principles calculation. \wannier\ 615calculates the band energies, through Wannier interpolation, on a 616dense mesh of k-points in the Brillouin zone. The density of this grid 617is controlled by the keyword {\tt fermi\_surface\_num\_points}. The 618default value is 50 (i.e., 50$^3$ points). The Fermi surface file 619{\tt copper.bxsf} can be viewed using {\tt XCrySDen}, e.g., 620% 621{\tt 622\begin{quote} 623xcrysden \texttt{-{}-}bxsf copper.bxsf 624\end{quote} } 625 626 627\item Plot the interpolated bandstructure. A suitable path in k-space is 628\smallskip 629{\tt 630\begin{quote} 631begin kpoint\_path 632 633G 0.00 0.00 0.00 X 0.50 0.50 0.00 634 635X 0.50 0.50 0.00 W 0.50 0.75 0.25 636 637W 0.50 0.75 0.25 L 0.00 0.50 0.00 638 639L 0.00 0.50 0.00 G 0.00 0.00 0.00 640 641G 0.00 0.00 0.00 K 0.00 0.50 -0.50 642 643end kpoint\_path 644\end{quote} } 645\end{enumerate} 646 647\subsection*{Further ideas} 648\begin{itemize} 649\item Compare the Wannier interpolated bandstructure with the full 650 \pwscf\ bandstructure. Obtain MLWFs using a denser k-point grid. 651To plot the bandstructure you can use the \pwscf\ tool {\tt bands.x} or the small FORTRAN program available at \url{http://www.tcm.phy.cam.ac.uk/~jry20/bands.html}. 652\item Investigate the effects of the outer and inner energy windows on 653 the interpolated bands. 654\item Instead of extracting a subspace of seven states, we could 655 extract a nine dimensional space (i.e., with $s$, $p$ and $d$ 656 character). Examine this case and compare the interpolated 657 bandstructures. 658\end{itemize} 659 660%\begin{figure}[h] 661%\begin{center} 662%\includegraphics{cu} 663%\caption{Band Structure of Copper showing the position of the outer 664% and inner energy windows.} 665%\label{fig:cu-bnd} 666%\end{center} 667%\end{figure} 668 669%\cleardoublepage 670 671\sectiontitle{7: Silane (SiH$_4$) -- Molecular MLWFs using 672 $\Gamma$-point sampling} 673\begin{itemize} 674\item{Outline: \it{Obtain MLWFs for the occupied states of molecular 675 silane. $\Gamma$-point sampling}} 676\item{Directory: {\tt examples/example7/}} 677\item{Input Files} 678\begin{itemize} 679\item{ {\tt silane.scf} {\it The \pwscf\ input file for ground state 680 calculation}} 681\item{ {\tt silane.nscf} {\it The \pwscf\ input file to obtain Bloch states 682 on a uniform grid}} 683\item{ {\tt silane.pw2wan} {\it Input file for {\tt pw2wannier90}}} 684\item{ {\tt silane.win} {\it The {\tt wannier90} input file}} 685\end{itemize} 686\end{itemize} 687 688\begin{enumerate} 689\item Run \pwscf\ to obtain the ground state of silane\\ 690{\tt pw.x < silane.scf > scf.out} 691 692\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 693{\tt pw.x < silane.nscf > nscf.out} 694 695\item Run \wannier\ to generate a list of the required overlaps (written 696 into the {\tt silane.nnkp} file).\\ 697{\tt wannier90.x -pp silane} 698 699\item Run {\tt pw2wannier90} to compute the overlap between Bloch 700 states and the projections for the starting guess (written in the 701 {\tt silane.mmn} and {\tt silane.amn} files).\\ 702{\tt pw2wannier90.x < silane.pw2wan > pw2wan.out} 703 704\item Run \wannier\ to compute the MLWFs.\\ 705{\tt wannier90.x silane} 706\end{enumerate} 707 708%\cleardoublepage 709 710\sectiontitle{8: Iron -- Spin-polarized WFs, DOS, projected WFs versus MLWFs} 711\begin{itemize} 712\item{Outline: \it{Generate both maximally-localized and projected 713 Wannier functions for ferromagnetic bcc Fe. Calculate the total 714 and orbital-projected density of states by Wannier 715 interpolation.}} 716\item{Directory: {\tt examples/example8/}} 717\item{Input Files} 718\begin{itemize} 719\item{ {\tt iron.scf} {\it The \pwscf\ input file for the 720 spin-polarized ground state calculation}} 721\item{ {\tt iron.nscf} {\it The \pwscf\ input file to obtain Bloch states 722 on a uniform grid}} 723\item{ {\tt iron\_\{up,down\}.pw2wan} {\it Input files for {\tt pw2wannier90}}} 724\item{ {\tt iron\_\{up,down\}.win} {\it Input files for {\tt wannier90} and {\tt 725 postw90}}} 726\end{itemize} 727\item{Note that in a spin-polarized calculation the spin-up and 728 spin-down MLWFs are computed separately. (The more general case of 729 spinor WFs will be treated in Example~17.)} 730\end{itemize} 731 732\begin{enumerate} 733\item Run \pwscf\ to obtain the ferromagnetic ground state of bcc Fe\\ 734{\tt pw.x < iron.scf > scf.out} 735 736\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 737{\tt pw.x < iron.nscf > nscf.out} 738 739\item Run \wannier\ to generate a list of the required overlaps (written 740 into the {\tt .nnkp} files).\\ 741{\tt wannier90.x -pp iron\_up}\\ 742{\tt wannier90.x -pp iron\_dn} 743 744\item Run {\tt pw2wannier90} to compute the overlap between Bloch 745 states and the projections for the starting guess (written in the 746 {\tt .mmn} and {\tt .amn} files).\\ 747 {\tt pw2wannier90.x < iron\_up.pw2wan > pw2wan\_up.out}\\ 748 {\tt pw2wannier90.x < iron\_dn.pw2wan > pw2wan\_dn.out} 749 750\item Run \wannier\ to compute the MLWFs.\\ 751{\tt wannier90.x iron\_up}\\ 752{\tt wannier90.x iron\_dn} 753 754\end{enumerate} 755 756\subsection*{Density of states} 757 758To compute the DOS using a $25\times 25 \times 25$ $k$-point grid add 759to the two {\tt .win} files 760% 761{\tt 762\begin{quote} 763dos = true 764 765dos\_kmesh = 25 766\end{quote} 767} 768% 769run {\tt postw90}, 770% 771{\tt 772\begin{quote} 773postw90.x iron\_up 774 775postw90.x iron\_dn 776\end{quote} 777} 778% 779and plot the DOS with {\tt gnuplot}, 780% 781{\tt 782\begin{quote} 783myshell> gnuplot 784 785gnuplot> plot `iron\_up\_dos.dat' u (-\$2):(\$1-12.6256) w l,`iron\_dn\_dos.dat' u 2:(\$1-12.6256) w l 786 787\end{quote} 788} 789% 790Energies are referred to the Fermi level (12.6256~eV, from {\tt 791 scf.out}). Note the exchange splitting between the up-spin and 792down-spin DOS. Check the convergence by repeating the DOS calculations 793with more $k$-points. 794 795\subsection*{Projected versus maximally-localized Wannier functions} 796 797In the calculations above we chose $s$, $p$, and $d$-type trial 798orbitals in the {\tt .win} files, 799% 800{\tt 801\begin{quote} 802Fe:s;p;d 803\end{quote} 804} 805% 806Let us analyze the evolution of the WFs during the gauge-selection 807step. Open one of the {\tt .wout} files and search for ``{\tt Initial 808 state}''; those are the {\it projected} WFs. As expected they are 809atom-centred, with spreads organized in three groups, 1+3+5: one $s$, 810three $p$, and five $d$. Now look at the final state towards the end 811of the file. The Wannier spreads have re-organized in two groups, 8126+3; moreover, the six more diffuse WFs are off-centred: the initial 813atomic-like orbitals hybridized with one another, becoming more 814localized in the process. It is instructive to visualize the 815final-state MLWFs using {\tt XCrySDen}, following Example~1. For more 816details, see Sec.~IV.B of Ref.~\cite{wang-prb06}. 817 818Let us plot the evolution of the spread functional~$\Omega$, 819% 820{\tt 821\begin{quote} 822myshell> grep SPRD iron\_up.wout > sprd\_up 823 824myshell> gnuplot 825 826gnuplot> plot `sprd\_up' u 6 w l 827\end{quote} 828} 829 830\begin{figure}[h] 831\begin{center} 832\scalebox{0.35}{\includegraphics{Fe-spread}} 833\caption{Evolution of the Wannier spread $\Omega$ of the minority 834 (spin-up) bands of bcc Fe during the iterative minimization of 835 $\widetilde{\Omega}$, starting from $s$, $p$, and $d$-type trial 836 orbitals.} 837\label{fig:Fe-sprd} 838\end{center} 839\end{figure} 840 841 842The first plateau corresponds to atom-centred WFs of separate $s$, 843$p$, and $d$ character, and the sharp drop signals the onset of the 844hybridization. With hindsight, we can redo steps~4 and~5 more 845efficiently using trial orbitals with the same character as the final 846MLWFs, 847% 848{\tt 849\begin{quote} 850Fe:sp3d2;dxy;dxz,dyz 851\end{quote} 852} 853% 854With this choice the minimization converges much more rapidly. 855 856Any reasonable set of localized WFs spanning the states of interest 857can be used to compute physical quantities (they are 858``gauge-invariant''). Let us recompute the DOS using, instead of 859MLWFs, the WFs obtained by projecting onto $s$, $p$, and $d$-type 860trial orbitals, without further iterative minimization of the spread 861functional. This can be done by setting 862% 863{\tt 864\begin{quote} 865num\_iter = 0 866\end{quote} 867} 868% 869But note that we still need to do disentanglement! 870Recalculate the DOS to confirm that it is almost identical to the one 871obtained earlier using the hybridized set of MLWFs. Visualize the 872projected WFs using {\tt XCrySDen}, to see that they retain the pure 873orbital character of the individual trial orbitals. 874 875 876\subsection*{Orbital-projected DOS and exchange splitting} 877 878With projected WFs the total DOS can be separated into $s$, $p$ and 879$d$ contributions, in a similar way to the orbital decomposition of 880the energy bands in Example~4. 881 882In order to obtain the partial DOS projected onto the $p$-type WFs, 883add to the {\tt .win} files 884% 885{\tt 886\begin{quote} 887dos\_project = 2,3,4 888\end{quote} 889} 890% 891and re-run {\tt postw90}. Plot the projected DOS for both up- and 892down-spin bands. Repeat for the $s$ and $d$ projections. 893 894Projected WFs can also be used to quantify more precisely the exchange 895splitting between majority and minority states. Re-run {\tt wannier90} 896after setting {\tt dos=false} and adding to the {\tt .win} files 897% 898{\tt 899\begin{quote} 900write\_hr\_diag = true 901\end{quote} 902} 903% 904This instructs {\tt wannier90} to print in the output file the on-site 905energies $\langle {\bf 0}n\vert H\vert {\bf 0}n\rangle$. The 906difference between corresponding values in {\tt iron\_up.wout} and in 907{\tt iron\_dn.wout} gives the exchange splittings for the individual 908orbitals. Compare their magnitudes with the splittings displayed by 909the orbital-projected DOS plots. In agreement with the Stoner 910criterion, the largest exchange splittings occur for the localized 911$d$-states, which contribute most of the density of states at the 912Fermi level. 913 914 915%\cleardoublepage 916 917 918\sectiontitle{9: Cubic BaTiO$_3$} 919\begin{itemize} 920\item{Outline: \it{Obtain MLWFs for a perovskite}} 921\item{Directory: {\tt examples/example9/}} 922\item{Input Files} 923\begin{itemize} 924\item{ {\tt batio3.scf} {\it The \pwscf\ input file for ground state 925 calculation}} 926\item{ {\tt batio3.nscf} {\it The \pwscf\ input file to obtain Bloch states 927 on a uniform grid}} 928\item{ {\tt batio3.pw2wan} {\it Input file for {\tt pw2wannier90}}} 929\item{ {\tt batio3.win} {\it The {\tt wannier90} input file}} 930\end{itemize} 931\end{itemize} 932 933 To start with, we are going to obtain MLWFs for the oxygen 2p 934 states. From the bandstructure~\cite{marzari-arxiv98}, these form an isolated 935 group of bands. We use the \wannier\ keyword {\tt exclude\_bands} to 936 remove all but the 2p bands from the calculation of the overlap 937 and projection matrices (we don't have to do this, but it saves time). 938 939\begin{enumerate} 940\item Run \pwscf\ to obtain the ground state of BaTiO$_3$\\ 941{\tt pw.x < BaTiO3.scf > scf.out} 942 943\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 944{\tt pw.x < BaTiO3.nscf > nscf.out} 945 946\item Run \wannier\ to generate a list of the required overlaps (written 947 into the {\tt BaTiO3.nnkp} file).\\ 948{\tt wannier90.x -pp BaTiO3} 949 950\item Run {\tt pw2wannier90} to compute the overlap between Bloch 951 states and the projections for the starting guess (written in the 952 {\tt BaTiO3.mmn} and {\tt BaTiO3.amn} files).\\ 953{\tt pw2wannier90.x < BaTiO3.pw2wan > pw2wan.out} 954 955\item Run \wannier\ to compute the MLWFs.\\ 956{\tt wannier90.x BaTiO3} 957\end{enumerate} 958 959Inspect the output file {\tt BaTiO3.wout}. 960 961Plot the second MLWF, as described in Section~1, by adding the 962following keywords to the input file {\tt BaTiO3.win} 963{\tt 964\begin{quote} 965wannier\_plot = true\\ 966restart = plot\\ 967wannier\_plot\_list = 2\\ 968wannier\_plot\_supercell = 3 969\end{quote} } 970and re-running \wannier. Visualise it using {\tt XCrySDen}, 971{\tt 972\begin{quote} 973xcrysden \texttt{-{}-}xsf BaTiO3\_00002.xsf 974\end{quote} } 975 976We can now simulate the ferroelectric phase by displacing the Ti 977 atom. Change its position to 978{\tt 979\begin{quote} 980Ti 0.505 0.5 0.5 981\end{quote} 982} 983and regenerate the MLWFs (i.e., compute the ground-state charge 984density and Bloch states using \pwscf, etc.) and look at the change in 985the second MLWF. 986 987\subsection*{Further ideas} 988\begin{itemize} 989\item Look at MLWFs for other groups of bands. What happens if you form 990 MLWFs for the whole valence manifold? 991 992\item Following Ref.~\cite{marzari-arxiv98}, compute the Born effective charges from the 993 change in Wannier centres under an atomic displacement. 994\end{itemize} 995 996%\cleardoublepage 997 998\sectiontitle{10: Graphite} 999\begin{itemize} 1000\item{Outline: \it{Obtain MLWFs for graphite (AB, Bernal)}} 1001\item{Directory: {\tt examples/example10/}} 1002\item{Input Files} 1003\begin{itemize} 1004\item{ {\tt graphite.scf} {\it The \pwscf\ input file for ground 1005 state calculation}} 1006\item{ {\tt graphite.nscf} {\it The \pwscf\ input file to obtain Bloch 1007 states on a uniform grid}} 1008\item{ {\tt graphite.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1009\item{ {\tt graphite.win} {\it The {\tt wannier90} input file}} 1010\end{itemize} 1011\end{itemize} 1012 1013\begin{enumerate} 1014\item Run \pwscf\ to obtain the ground state of graphite\\ 1015{\tt pw.x < graphite.scf > scf.out} 1016 1017\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 1018{\tt pw.x < graphite.nscf > nscf.out} 1019 1020\item Run \wannier\ to generate a list of the required overlaps (written 1021 into the {\tt graphite.nnkp} file).\\ 1022{\tt wannier90.x -pp graphite} 1023 1024\item Run {\tt pw2wannier90} to compute the overlap between Bloch 1025 states and the projections for the starting guess (written in the 1026 {\tt graphite.mmn} and {\tt graphite.amn} files).\\ 1027{\tt pw2wannier90.x < graphite.pw2wan > pw2wan.out} 1028 1029\item Run \wannier\ to compute the MLWFs.\\ 1030{\tt wannier90.x graphite} 1031\end{enumerate} 1032 1033 1034%\cleardoublepage 1035 1036 1037\sectiontitle{11: Silicon -- Valence and low-lying conduction states} 1038\subsection*{Valence States} 1039\begin{itemize} 1040\item{Outline: \it{Obtain MLWFs for the valence bands of silicon.}} 1041\item{Directory: {\tt examples/example11/}} 1042\item{Input Files} 1043\begin{itemize} 1044\item{ {\tt silicon.scf} {\it The \pwscf\ input file for ground state 1045 calculation}} 1046\item{ {\tt silicon.nscf} {\it The \pwscf\ input file to obtain Bloch 1047 states on a uniform grid}} 1048\item{ {\tt silicon.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1049\item{ {\tt silicon.win} {\it The {\tt wannier90} input file}} 1050\end{itemize} 1051 1052\end{itemize} 1053 1054\begin{enumerate} 1055\item Run \pwscf\ to obtain the ground state of silicon\\ 1056{\tt pw.x < silicon.scf > scf.out} 1057 1058\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 1059 grid. Note that we request the lower 4 (valence) bands\\ 1060{\tt pw.x < silicon.nscf > nscf.out} 1061 1062\item Run \wannier\ to generate a list of the required overlaps (written 1063 into the {\tt silicon.nnkp} file).\\ 1064{\tt wannier90.x -pp silicon} 1065 1066\item Run {\tt pw2wannier90} to compute the overlap between Bloch 1067 states and the projections for the starting guess (written in the 1068 {\tt silicon.mmn} and {\tt silicon.amn} files).\\ 1069{\tt pw2wannier90.x < silicon.pw2wan > pw2wan.out} 1070 1071\item Run \wannier\ to compute the MLWFs.\\ 1072{\tt wannier90.x silicon} 1073 1074\end{enumerate} 1075 1076Inspect the output file {\tt silicon.wout}. The total spread converges to its 1077minimum value after just a few iterations. Note that the geometric centre of 1078each MLWF lies at the centre of the Si-Si bond. 1079Note also that the memory requirement for the minimisation of 1080the spread is very low as the MLWFs are defined 1081by just the 4$\times$4 unitary matrices $\Uk$. 1082 1083Plot the MLWFs by adding the following keywords to the input file {\tt 1084 silicon.win} 1085{\tt 1086\begin{quote} 1087wannier\_plot = true 1088\end{quote} } 1089and re-running \wannier. Visualise them using {\tt XCrySDen}, e.g., 1090{\tt 1091\begin{quote} 1092xcrysden \texttt{-{}-}xsf silicon\_00001.xsf 1093\end{quote} } 1094 1095\subsection*{Valence + Conduction States} 1096 1097\begin{itemize} 1098\item{Outline: \it{Obtain MLWFs for the valence and low-lying 1099 conduction-band states of Si. Plot the interpolated 1100 bandstructure. Apply a scissors correction to the conduction 1101 bands.}} 1102\item{Input Files} 1103\begin{itemize} 1104\item{ {\tt silicon.scf} {\it The \pwscf\ input file for ground state 1105 calculation}} 1106\item{ {\tt silicon.nscf} {\it The \pwscf\ input file to obtain Bloch 1107 states on a uniform grid}} 1108\item{ {\tt silicon.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1109\item{ {\tt silicon.win} {\it The {\tt wannier90} input file}} 1110\end{itemize} 1111\end{itemize} 1112The valence and lower conduction states can be represented by MLWFs 1113with $sp^3$-like symmetry. The lower conduction states are not 1114separated by an energy gap from the higher states. In order to form 1115localised WF we use the disentanglement procedure introduced in 1116Ref.~\cite{souza-prb01}. The position of the inner and outer energy 1117windows are shown in Fig.~\ref{fig:si.bnd}. 1118\begin{enumerate} 1119\item Modify the input file and run \pwscf\ and \wannier.\\ 1120Inspect the output file {\tt silicon.wout}. The minimisation of the 1121spread occurs in a two-step procedure. First, we minimise $\Omega_{\rm 1122 I}$ -- this is the extraction of the optimal subspace in the 1123disentanglement procedure. Then, we minimise $\Omega_{\rm 1124 O}+\Omega_{{\rm OD}}$. 1125 1126\item Plot the bandstructure by adding the following commands to the 1127 input file {\tt silicon.win} 1128{\tt 1129\begin{quote} 1130restart = plot 1131 1132bands\_plot = true 1133\end{quote} } 1134and re-running \wannier. The files {\tt silicon\_band.dat} and {\tt 1135 silicon\_band.gnu} are created. To plot the bandstructure using 1136 gnuplot \smallskip 1137{\tt 1138\begin{quote} 1139myshell> gnuplot 1140 1141gnuplot> load `silicon\_band.gnu' 1142\end{quote} } 1143The k-point path for the bandstructure interpolation is set in the {\tt 1144 kpoint\_path} block. Try plotting along different paths. 1145 1146%\item Shift the conduction bands upwards in energy by 1147% 0.65~eV. This {\it ad-hoc} ``scissors correction'' is applied to 1148% the Hamiltonian in the Wannier basis using the instructions 1149% 1150%{\tt 1151%\begin{quote} 1152%num\_valence\_bands = 4 1153% 1154%scissors\_shift = 0.65 1155%\end{quote} } 1156% 1157%Plot together the original and scissors-corrected bands. 1158 1159\end{enumerate} 1160 1161\subsection*{Further ideas} 1162 1163\begin{itemize} 1164\item Compare the Wannier-interpolated bandstructure with the full 1165 \pwscf\ bandstructure. Recompute the MLWFs using a finer $k$-point 1166 grid (e.g., 6$\times$6$\times$6 or 8$\times$8$\times$8) and note how 1167 the accuracy of the interpolation increases~\cite{yates-prb07}. 1168\item Compute four MLWFs spanning the low-lying conduction states (see 1169 Ref.~\cite{souza-prb01}). 1170\end{itemize} 1171 1172%\begin{figure}[h] 1173%\begin{center} 1174%\includegraphics{si} 1175%\caption{Band Structure of Silicon showing the position of the outer 1176%and inner energy windows.} 1177%\label{fig:si.bnd} 1178%\end{center} 1179%\end{figure} 1180 1181%\cleardoublepage 1182 1183 1184\sectiontitle{12: Benzene -- Valence and low-lying conduction states} 1185\subsection*{Valence States} 1186\begin{itemize} 1187\item{Outline: \it{Obtain MLWFs for the valence states of benzene}} 1188\item{Directory: {\tt examples/example12/}} 1189\item{Input Files} 1190\begin{itemize} 1191\item{ {\tt benzene.scf} {\it The \pwscf\ input file for ground state 1192 calculation}} 1193\item{ {\tt benzene.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1194\item{ {\tt benzene.win} {\it The {\tt wannier90} input file}} 1195\end{itemize} 1196 1197\end{itemize} 1198 1199\begin{enumerate} 1200\item Run \pwscf\ to obtain the ground state of benzene\\ 1201{\tt pw.x < benzene.scf > scf.out} 1202 1203\item Run \wannier\ to generate a list of the required overlaps (written 1204 into the {\tt benzene.nnkp} file).\\ 1205{\tt wannier90.x -pp benzene} 1206 1207\item Run {\tt pw2wannier90} to compute the overlap between Bloch 1208 states and the projections for the starting guess (written in the 1209 {\tt benzene.mmn} and {\tt benzene.amn} files).\\ 1210{\tt pw2wannier90.x < benzene.pw2wan > pw2wan.out} 1211 1212\item Run \wannier\ to compute the MLWFs.\\ 1213{\tt wannier90.x benzene} 1214 1215\end{enumerate} 1216 1217Inspect the output file {\tt benzene.wout}. The total spread converges 1218to its minimum value after just a few iterations. 1219 1220Plot the MLWFs by adding the following keywords to the input file {\tt 1221 benzene.win} 1222{\tt 1223\begin{quote} 1224restart = plot\\ 1225wannier\_plot = true\\ 1226wannier\_plot\_format = cube\\ 1227wannier\_plot\_list = 2-4 1228\end{quote} } 1229and re-running \wannier. Visualise them using, e.g., {\tt XCrySDen}. 1230 1231\subsection*{Valence + Conduction States} 1232 1233\begin{itemize} 1234\item{Outline: \it{Obtain MLWFs for the valence and low-lying 1235 conduction states of benzene.}} 1236\item{Input Files} 1237\begin{itemize} 1238\item{ {\tt benzene.scf} {\it The \pwscf\ input file for ground state 1239 calculation}} 1240\item{ {\tt benzene.nscf} {\it The \pwscf\ input file to obtain Bloch 1241 states for the conduction states}} 1242\item{ {\tt benzene.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1243\item{ {\tt benzene.win} {\it The {\tt wannier90} input file}} 1244\end{itemize} 1245\end{itemize} 1246In order to form localised WF we use the disentanglement 1247procedure. The position of the inner energy window is set to lie in 1248the energy gap; the outer energy window is set to 4.0\,eV. Modify the 1249input file appropriately. 1250\begin{enumerate} 1251\item Run \pwscf\ and \wannier.\\ 1252Inspect the output file {\tt benzene.wout}. The minimisation of the 1253spread occurs in a two-step procedure. First, we minimise $\Omega_{\rm 1254 I}$. Then, we minimise $\Omega_{\rm O}+\Omega_{{\rm OD}}$. 1255 1256\item Plot the MLWFs by adding the following commands to the 1257 input file {\tt benzene.win} 1258{\tt 1259\begin{quote} 1260restart = plot\\ 1261wannier\_plot = true\\ 1262wannier\_plot\_format = cube\\ 1263wannier\_plot\_list = 1,7,13 1264\end{quote} } 1265and re-running \wannier. Visualise them using, e.g., {\tt XCrySDen}. 1266\end{enumerate} 1267 1268%\cleardoublepage 1269 1270\sectiontitle{13: (5,5) Carbon Nanotube -- Transport properties} 1271%\subsection*{Transport properties} 1272 1273\begin{itemize} 1274 \item{Outline: \it{Obtain the bandstructure, quantum conductance and 1275 density of states of a metallic (5,5) carbon nanotube}} 1276 \item{Directory: {\tt examples/example13/}} 1277 \item{Input Files} 1278 \begin{itemize} 1279 \item{ {\tt cnt55.scf} {\it The \pwscf\ input file for ground state 1280 calculation}} 1281 \item{ {\tt cnt55.nscf} {\it The \pwscf\ input file to obtain Bloch 1282 states for the conduction states}} 1283 \item{ {\tt cnt55.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1284 \item{ {\tt cnt55.win} {\it The {\tt wannier90} input file}} 1285 \end{itemize} 1286\end{itemize} 1287 1288In order to form localised WF that describe both the occupied and 1289unoccupied $\pi$ and $\pi^{\ast}$ manifolds, we use the 1290disentanglement procedure to extract a smooth manifold of states that 1291has dimension equal to 2.5 times the number of carbon atoms per unit 1292cell~\cite{lee-prl05}. The positions of the energy windows are shown in 1293Fig.~\ref{fig:cnt.win}. 1294 1295The part of the \wannier\ input file that controls the transport part 1296of the calculation looks like: 1297 1298{\tt 1299\begin{quote} 1300transport = true\\ 1301transport\_mode = bulk\\ 1302one\_dim\_axis = z\\ 1303dist\_cutoff = 5.5\\ 1304fermi\_energy = -1.06\\ 1305tran\_win\_min = -6.5\\ 1306tran\_win\_max = 6.5\\ 1307tran\_energy\_step = 0.01\\ 1308dist\_cutoff\_mode = one\_dim\\ 1309translation\_centre\_frac = 0.0 0.0 0.0 1310\end{quote} } 1311 1312Descriptions of these and other keywords related to the calculation of 1313transport properties can be found in the User Guide. 1314 1315\begin{enumerate} 1316\item Run \pwscf\ and \wannier.\\ 1317Inspect the output file {\tt cnt55.wout}. The minimisation of the 1318spread occurs in a two-step procedure. First, we minimise $\Omega_{\rm 1319 I}$. Then, we minimise $\Omega_{\rm O}+\Omega_{{\rm OD}}$. 1320\item Note that the initial $p_{z}$ projections on the carbon atoms 1321are oriented in the radial direction with respect to the nanotube 1322axis. 1323\item The interpolated bandstructure is written to {\tt 1324cnt55\_band.agr} (since {\tt bands\_plot\_format = xmgr} in the input 1325file). 1326\item The quantum conductance and density of states are written to the 1327files {\tt cnt55\_qc.dat} and {\tt cnt55\_dos.dat}, respectively. 1328Note that this part of the calculation may take some time. You can 1329follow its progress by monitoring the output to these files. 1330Use a package such as {\tt gnuplot} or {\tt xmgrace} in order to visualise 1331the data. You should get something that looks like Fig.~\ref{fig:cnt.tran}. 1332\end{enumerate} 1333 1334\begin{figure}[h] 1335\begin{center} 1336\includegraphics[width=10cm]{cnt_win} 1337\caption{Bandstructure of (5,5) carbon nanotube showing the position 1338 of the outer and inner energy windows.} 1339\label{fig:cnt.win} 1340\end{center} 1341\end{figure} 1342 1343\begin{figure}[h] 1344\begin{center} 1345\scalebox{0.8}{\includegraphics{cnt_tran}} 1346\caption{Wannier interpolated bandstructure, quantum conductance and 1347density of states of (5,5) carbon nanotube. Note that the Fermi level has been shifted 1348by 1.06eV with respect to Fig.~\ref{fig:cnt.win}.} 1349\label{fig:cnt.tran} 1350\end{center} 1351\end{figure} 1352 1353%\cleardoublepage 1354 1355\sectiontitle{14: Linear Sodium Chain -- Transport properties} 1356 1357\begin{itemize} 1358 \item{Outline: \it{Compare the quantum conductance of a periodic 1359 linear chain of Sodium atoms with that of a defected chain}} 1360 \item{\begin{tabbing} 1361 Directories: \= {\tt examples/example14/periodic}\\ 1362 \> {\tt examples/example14/defected} 1363 \end{tabbing}} 1364 \item{Input Files} 1365 \begin{itemize} 1366 \item{ {\tt Na\_chain.scf} {\it The \pwscf\ input file for ground state 1367 calculation}} 1368 \item{ {\tt Na\_chain.nscf} {\it The \pwscf\ input file to obtain Bloch 1369 states for the conduction states}} 1370 \item{ {\tt Na\_chain.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1371 \item{ {\tt Na\_chain.win} {\it The {\tt wannier90} input file}} 1372 \end{itemize} 1373\end{itemize} 1374 1375The periodic system contains two unit cells evenly distributed along 1376the supercell. Transport calculations are performed using {\tt 1377 transport\_mode = bulk} and so the resulting quantum conductance 1378represents that of an infinite periodic chain. 1379 1380The part of the \wannier\ input file that controls the transport part 1381of the calculation looks like: 1382 1383{\tt 1384\begin{quote} 1385transport = true\\ 1386transport\_mode = bulk\\ 1387tran\_read\_ht = false\\ 1388one\_dim\_axis = x\\ 1389fermi\_energy = -2.7401\\ 1390tran\_win\_min = -5.0\\ 1391tran\_win\_max = 5.0\\ 1392tran\_energy\_step = 0.01\\ 1393translation\_centre\_frac = 0.5 0.5 0.5\\ 1394tran\_num\_bb = 2 1395 1396\end{quote} } 1397 1398The defected system uses a 13 atom supercell with the central atom 1399position altered to break symmetry. Setting {\tt transport\_mode = lcr} with tell 1400\wannier\ to treat the system as an infinite system with the defect at its centre. 1401The supercell is chosen so that is conforms to the 2c2 geometry (see User Guide 1402for details). Each principal layer is 2 atoms long so that the conductor 1403region contains the defected atom plus a single atom on either side. 1404 1405The transport section of the input file contains these key differences: 1406 1407{\tt 1408\begin{quote} 1409transport\_mode = lcr\\ 1410tran\_num\_ll = 2\\ 1411tran\_num\_cell\_ll = 2\\ 1412 1413\end{quote} } 1414 1415Descriptions of these and other keywords related to the calculation of 1416transport properties can be found in the User Guide. 1417 1418\begin{enumerate} 1419\item Run \pwscf\ and \wannier\ for the periodic system. 1420\item Run \pwscf\ and \wannier\ for the defected system. 1421\item The quantum conductance is written to the files {\tt periodic/Na\_chain\_qc.dat} 1422and \linebreak{\tt defected/Na\_chain\_dos.dat}, respectively. 1423Compare the quantum conductance of the periodic (bulk) calculation with the 1424defected (LCR) calculation. Your plot should look like Fig.~\ref{fig:Na_qc}. 1425\end{enumerate} 1426 1427 1428\begin{figure}[h] 1429\begin{center} 1430\scalebox{0.3}{\includegraphics{Na_qc}} 1431\caption{Quantum conductance of periodic Sodium chain (black) compared to that of the defected Sodium chain (red).} 1432\label{fig:Na_qc} 1433\end{center} 1434\end{figure} 1435 1436%\cleardoublepage 1437 1438\sectiontitle{15: (5,0) Carbon Nanotube -- Transport properties} 1439 1440\emph{Note that these systems require reasonably large-scale electronic 1441structure calculations.} 1442 1443\subsection*{Bulk Transport properties} 1444 1445\begin{itemize} 1446 \item{Outline: \it{Obtain the quantum conductance of a pristine single-walled carbon nanotube}} 1447 \item{Directory: {\tt examples/example14/periodic}} 1448 \item{Input Files} 1449 \begin{itemize} 1450 \item{ {\tt cnt.scf} {\it The \pwscf\ input file for ground state 1451 calculation}} 1452 \item{ {\tt cnt.nscf} {\it The \pwscf\ input file to obtain Bloch 1453 states for the conduction states}} 1454 \item{ {\tt cnt.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1455 \item{ {\tt cnt.win} {\it The {\tt wannier90} input file}} 1456 \end{itemize} 1457\end{itemize} 1458 1459 1460First we consider a single unit cell, with 10 k-points. With 1461{\tt transport\_mode = bulk} we compute the transport properties 1462of a pristine, infinite, periodic (5,0) carbon nanotube. Later, we will 1463compare the quantum conductance of this system with a defected 1464nanotube. 1465 1466\begin{enumerate} 1467\item Run \pwscf\ and \wannier.\\ 1468\item The quantum conductance and density of states are written to the 1469files {\tt cnt\_qc.dat} and {\tt cnt\_dos.dat}, respectively. 1470\end{enumerate} 1471 1472\subsection*{LCR transport properties -- Defected nanotube} 1473 1474\begin{itemize} 1475 \item{Outline: \it{Use the automated LCR routine to investigate the effect of 1476 a single silicon atom in a infinite (5,0) carbon nanotube.}} 1477 \item{Directory: {\tt examples/example15/defected}} 1478 \item{Input Files} 1479 \begin{itemize} 1480 \item{ {\tt cnt+si.scf} {\it The \pwscf\ input file for ground state 1481 calculation}} 1482 \item{ {\tt cnt+si.nscf} {\it The \pwscf\ input file to obtain Bloch 1483 states for the conduction states}} 1484 \item{ {\tt cnt+si.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1485 \item{ {\tt cnt+si.win} {\it The {\tt wannier90} input file}} 1486 \end{itemize} 1487\end{itemize} 1488 1489In this calculation an 11-atom supercell is used with a single silicon 1490substitutional defect in the central unit cell. The supercell is 1491chosen so that is conforms to the 2c2 geometry (see User Guide for 1492details) with principal layers set to be two unit cells long. 1493 1494\begin{enumerate} 1495\item Run \pwscf\ and \wannier. Again these are large calculations, progress 1496can be monitored by viewing respective output files.\\ 1497\item The quantum conductance is written to {\tt cnt+si\_qc.dat}. 1498Compare the quantum conductance with the periodic (bulk) calculation. 1499Your plot should look like Fig.~\ref{fig:cnt_qc}. 1500 1501\end{enumerate} 1502 1503\begin{figure}[h] 1504\begin{center} 1505\scalebox{0.3}{\includegraphics{cnt_qc}} 1506\caption{Quantum conductance of infinite pristine nanotube (black) 1507compared to that of the infinite nanotube with the substitutional silicon 1508defect (red).} 1509\label{fig:cnt_qc} 1510\end{center} 1511\end{figure} 1512 1513\subsection*{Further ideas} 1514\begin{itemize} 1515\item Set {\tt write\_hr = true} in the bulk case. Consider the magnitude of Hamiltonian 1516 elements between Wannier functions in increasingly distant unit cells. Are two 1517 unit cell principal layers really large enough, or are significant errors introduced? 1518\item Does one unit cell either side of the defected unit cell shield the disorder 1519 so that the leads are ideal? Does the quantum conductance change if these 1520 `buffer' regions are increased? 1521\end{itemize} 1522 1523%\cleardoublepage 1524 1525\sectiontitle{16: Silicon -- Boltzmann transport} 1526\begin{itemize} 1527\item{Outline: \it{Obtain MLWFs for the valence and low-lying 1528 conduction states of Si. Calculate the electrical conductivity, the 1529 Seebeck coefficient and the thermal conductivity in the constant 1530 relaxation time approximation using the \bw\ module.}} 1531\end{itemize} 1532\subsection*{If you want to use Quantum ESPRESSO} 1533\begin{itemize} 1534\item{Directory: {\tt examples/example16-withqe/}} 1535\item{Input Files} 1536\begin{itemize} 1537\item{ {\tt Si.scf} {\it The \pwscf\ input file for ground state 1538 calculation}} 1539\item{ {\tt Si.nscf} {\it The \pwscf\ input file to obtain Bloch 1540 states on a uniform grid}} 1541\item{ {\tt Si.pw2wan} {\it Input file for {\tt pw2wannier90}}} 1542\item{ {\tt Si.win} {\it The \wannier\ and \postw\ input file}}\ 1543\end{itemize} 1544\end{itemize} 1545\subsection*{If you do not want to use Quantum ESPRESSO} 1546\begin{itemize} 1547\item{Directory: {\tt examples/example16-noqe/}} 1548\item{Input Files} 1549\begin{itemize} 1550\item{ {\tt Si.win} {\it The \wannier\ and \postw\ input file}} 1551\item{ {\tt Si.mmn} {\it The overlap matrices $\Mkb$}} 1552\item{ {\tt Si.amn} {\it Projection $\Ak$ of the Bloch states onto a set 1553 of trial localised orbitals}} 1554\item{ {\tt Si.eig} {\it The Bloch eigenvalues at each k-point. For 1555 interpolation only}} 1556\end{itemize} 1557\end{itemize} 1558 1559Note the first five steps in the following are the same of Example~11, and 1560are needed only if you want to use the \texttt{PWscf} code of Quantum ESPRESSO. 1561Otherwise, if you have already run Example~11 with Quantum ESPRESSSO 1562(in particular, the section ``\emph{Valence + Conduction States}'') 1563you can start from those files and continue from point 6, after having 1564added the \bw\ flags to the input file. 1565 1566If instead you do not have Quantum ESPRESSO installed, or you do not want to 1567use it, you can start from step 5 using the files in the 1568{\tt examples/example16-noqe/} folder. 1569 1570\begin{enumerate} 1571\item Run \pwscf\ to obtain the ground state of silicon\\ 1572{\tt pw.x < Si.scf > scf.out} 1573 1574\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 1575 grid. Details on the disentanglement procedure are discussed in Example~11.\\ 1576{\tt pw.x < Si.nscf > nscf.out} 1577 1578\item Run \wannier\ to generate a list of the required overlaps (written 1579 into the {\tt Si.nnkp} file).\\ 1580{\tt wannier90.x -pp Si} 1581 1582\item Run {\tt pw2wannier90} to compute the overlap between Bloch 1583 states and the projections for the starting guess (written in the 1584 {\tt Si.mmn} and {\tt Si.amn} files).\\ 1585{\tt pw2wannier90.x < Si.pw2wan > pw2wan.out} 1586 1587\item Run \wannier\ to compute the MLWFs.\\ 1588{\tt wannier90.x Si} 1589 1590Inspect the output file {\tt Si.wout} and check if the convergence was reached both in the 1591disentanglement and in the wannierisation steps (as discussed in further detail in Example~11). 1592You may also want to plot the Wannier functions and the interpolated band structure. 1593 1594\item Run \postw\ to calculate the transport coefficients.\\ 1595{\tt postw90.x Si} (serial execution) \\ 1596{\tt mpirun -np 8 postw90.x Si} (example of parallel execution with 8 1597MPI processes) 1598\end{enumerate} 1599 1600Inspect the output file {\tt Si.wpout}. It summarizes the main details of the calculation (more details can be obtained by setting a larger value of the \verb#iprint# flag). 1601Check if no warnings are issued. Note that if no special flags are passed to \bw, it assumes that 1602the ab-initio calculation did not include magnetization effects, and thus it sets to 2 the 1603number of electrons per state. 1604 1605Note also that the value of the relaxation time $\tau=10$~fs in the 1606example is set only as a representative 1607value; note also that only the electrical and thermal conductivity 1608depend on $\tau$, while the Seebeck coefficient is independent of 1609$\tau$. 1610 1611Using your favourite plotting program, plot the {\tt Si\_boltzdos.dat} file to inspect the DOS. 1612 1613Using your favourite plotting program, plot columns 1 and 3 of the {\tt Si\_seebeck.dat} file to inspect the $S_{xx}$ component of the Seebeck coefficient as a function of the chemical potential $\mu$, at $T=300$~K. 1614 1615\subsection*{Further ideas} 1616 1617\begin{itemize} 1618\item Change the interpolation to a $60\times 60\times 60$ mesh and run again \postw\ to check if the results for the transport properties are converged. 1619 1620\item Change the {\tt Si.win} input file so that it calculates the transport coefficients for temperatures from 300 to 700~K, with steps of 200~K. Rerun \postw\ and verify that the increase in execution time is neglibile (in fact, most of the time is spent to interpolate the band structure on the $k$ mesh). 1621 1622Plot the Seebeck coefficient for the three temperatures $T=300$~K, $T=500$~K and $T=700$~K. To do this, you have to filter the {\tt Si\_seebeck.dat} to select only those lines where the second column is equal to the required temperature. A possible script to select the $S_{xx}$ component of the Seebeck coefficient for $T=500$~K using the {\tt awk/gawk} command line program is the following: 1623\begin{verbatim} 1624awk `{if ($2 == 500) {print $1, $3;}}' < Si_seebeck.dat \ 1625 > Si_seebeck_xx_500K.dat 1626\end{verbatim} 1627Then, you can plot columns 1 and 2 of the output file \verb#Si_seebeck_xx_500K.dat#. 1628\item Try to calculate the Seebeck coefficient as a function of the temperature, for a $n-$doped sample with, e.g., $n=10^{18}$ cm$^{-3}$. Note that to this aim, you need to calculate consistently the value $\mu(T)$ of the chemical potential as a function of the temperature, so as to reproduce the given value of $n$. Then, you have to write a small program/script to interpolate the output of \bw, that you should have run on a suitable grid of $(\mu,T)$ points. 1629\end{itemize} 1630 1631%\begin{figure}[h] 1632%\begin{center} 1633%\includegraphics{si} 1634%\caption{Band Structure of Silicon showing the position of the outer 1635%and inner energy windows.} 1636%\label{fig:si.bnd} 1637%\end{center} 1638%\end{figure} 1639 1640%\cleardoublepage 1641 1642 1643\sectiontitle{17: Iron -- Spin-orbit-coupled bands and 1644Fermi-surface contours} 1645 1646Note: It is recommended that you go through Example 8 first (bcc Fe 1647without spin-orbit). 1648 1649Note: This example requires a recent version of the {\tt pw2wannier90} interface. 1650 1651\begin{itemize} 1652\item{Outline: \it{Plot the spin-orbit-coupled bands of ferromagnetic 1653 bcc Fe. Plot the Fermi-surface contours on a plane in the 1654 Brillouin zone.}} 1655\item{Directory: {\tt examples/example17/}} 1656\item{Input files} 1657\begin{itemize} 1658\item{ {\tt Fe.scf} {\it The \pwscf\ input file for ground state 1659 calculation}} 1660\item{ {\tt Fe.nscf} {\it The \pwscf\ input file to obtain Bloch 1661 states on a uniform grid}} 1662\item{ {\tt Fe.pw2wan} {\it The input file for {\tt pw2wannier90}}} 1663 1664\item{ {\tt Fe.win} {\it The {\tt wannier90} and {\tt postw90} input file}} 1665\end{itemize} 1666\end{itemize} 1667 1668Note that {\tt num\_wann =18} in {\tt Fe.win}, but only nine trial 1669orbitals are provided. The line 1670 {\tt 1671\begin{quote} 1672spinors = true 1673\end{quote} 1674}tells {\tt wannier90} to use in step~3 below the specified trial 1675orbitals on both the up- and down-spin channels, effectively doubling 1676their number. 1677 1678\begin{enumerate} 1679\item Run \pwscf\ to obtain the ferromagnetic ground state of 1680 iron\footnote{Please note the following counterintuitive feature in 1681 {\tt pwscf}: in order to obtain a ground state with magnetization 1682 along the {\it positive} $z$-axis, one should use a {\it negative} 1683 value for the variable {\tt starting\_magnetization}.}\\ 1684 {\tt pw.x < Fe.scf > scf.out} 1685 1686\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 1687 grid\\ 1688{\tt pw.x < Fe.nscf > nscf.out} 1689 1690\item Run \wannier\ to generate a list of the required overlaps (written 1691 into the {\tt Fe.nnkp} file)\\ 1692{\tt wannier90.x -pp Fe} 1693 1694\item Run {\tt pw2wannier90} to compute: 1695 \begin{itemize} 1696 1697 \item[{\bf --}] The overlaps $\langle u_{n{\bf k}}\vert u_{m{\bf 1698 k}+{\bf b}}\rangle$ between {\it spinor} Bloch states (written 1699 in the {\tt Fe.mmn} file) 1700 1701 \item[{\bf --}] The projections for the starting guess (written in 1702 the {\tt Fe.amn} file) 1703 1704 \item[{\bf --}] The spin matrix elements $\langle \psi_{n{\bf 1705 k}}\vert \sigma_i\vert \psi_{m{\bf k}}\rangle$, $i=x,y,z$ 1706 (written in the {\tt Fe.spn} file) 1707%\item{The matrix elements $\langle u_{n{\bf k}+{\bf b}_1}\vert H_{\bf k}\vert 1708%u_{m{\bf k}+{\bf 1709% b}_2}\rangle$ (written in the {\tt 1710% Fe.uHu} file)} 1711 \end{itemize} 1712{\tt pw2wannier90.x < Fe.pw2wan > pw2wan.out} 1713 1714\item Run \wannier\ to compute the MLWFs.\\ 1715{\tt wannier90.x Fe} 1716 1717\item Run \postw\ to compute the energy eigenvalues and spin 1718 expectation values.\\ 1719 {\tt postw90.x Fe} (serial execution) \\ 1720 {\tt mpirun -np 8 postw90.x Fe} (example of parallel execution with 1721 8 MPI processes) 1722 1723\end{enumerate} 1724 1725 In this example we use the module {\tt kpath} to plot the energy 1726 bands coloured by the expectation value of the spin along [001]: 1727 {\tt 1728\begin{quote} 1729kpath = true 1730 1731kpath\_task = bands 1732 1733kpath\_bands\_colour = spin 1734 1735kpath\_num\_points=500 1736\end{quote} } 1737 1738 1739To plot the bands using {\tt gnuplot} (version 4.2 or higher) issue 1740% 1741{\tt 1742\begin{quote} 1743myshell> gnuplot 1744 1745gnuplot> load `Fe-bands.gnu' 1746\end{quote} } 1747% 1748or, using {\tt python}, 1749% 1750{\tt 1751\begin{quote} 1752myshell> python Fe-bands.py 1753\end{quote} } 1754 1755Next we plot the Fermi-surface contours on the (010) plane $k_y=0$, 1756using the {\tt kslice} module. Set {\tt kpath = false} and uncomment 1757the following instructions in {\tt Fe.win}, 1758% 1759{\tt 1760\begin{quote} 1761kslice = true 1762 1763kslice\_task = fermi\_lines 1764 1765fermi\_energy = [insert your value here] 1766 1767kslice\_corner = 0.0 0.0 0.0 1768 1769kslice\_b1 = 0.5 -0.5 -0.5 1770 1771kslice\_b2 = 0.5 0.5 0.5 1772 1773kslice\_2dkmesh = 200 200 1774\end{quote} } 1775 1776taking the Fermi level value from {\tt scf.out}. The energy 1777eigenvalues are computed on a $200\times 200$ $k$-point grid covering 1778the BZ slice. The lines of intersection between the Fermi surface and 1779the (010) plane can be visualized with the {\tt gnuplot} or {\tt 1780 python} scripts generated at runtime, 1781% 1782{\tt 1783\begin{quote} 1784myshell> gnuplot 1785 1786gnuplot> load `Fe-kslice-fermi\_lines.gnu' 1787\end{quote} } 1788% 1789%(do not be concerned by the warning messages, they are caused by the 1790%fact that not all bands cross the Fermi energy) 1791or 1792% 1793{\tt 1794\begin{quote} 1795myshell> python Fe-kslice-fermi\_lines.py 1796\end{quote} } 1797% 1798The Fermi lines can be colour-coded by the spin expectation value 1799$\langle S_z\rangle$ of the states on the Fermi surface. Add to {\tt 1800 Fe.win} the line {\tt 1801\begin{quote} 1802kslice\_fermi\_lines\_colour = spin 1803\end{quote} } 1804% 1805and re-run {\tt postw90}. The names of the {\tt gnuplot} and {\tt 1806 python} scripts generated at runtime are unchanged. (However, the 1807plotting algorithm is different in this case, and the lines are not as 1808smooth as before. You may want to increase {\tt kslice\_2dkmesh}.) 1809 1810 1811\subsection*{Further ideas} 1812 1813\begin{itemize} 1814 1815\item Redraw the Fermi surface contours on the (010) plane starting 1816 from a calculation without spin-orbit coupling, by adding to the 1817 input files {\tt iron\_\{up,down\}.win} in Example~8 the lines {\tt 1818\begin{quote} 1819kslice = true 1820 1821kslice\_task = fermi\_lines 1822 1823fermi\_energy = [insert your value here] 1824 1825kslice\_corner = 0.0 0.0 0.0 1826 1827kslice\_b1 = 0.5 -0.5 -0.5 1828 1829kslice\_b2 = 0.5 0.5 0.5 1830 1831kslice\_2dkmesh = 200 200 1832\end{quote} } 1833% 1834before running {\tt postw90}, 1835% 1836{\tt 1837\begin{quote} 1838postw90.x iron\_up 1839 1840postw90.x iron\_dn 1841\end{quote} 1842} 1843% 1844The {\tt python} scripts generated at runtime draw the up- and 1845down-spin Fermi lines on separate figures. To draw them together, use 1846the script {\tt iron\_updn-kslice-fermi\_lines.py} provided with 1847Example~17 (or merge the two generated scripts). Compare the Fermi 1848lines with and without spin-orbit, and note the spin-orbit-induced 1849avoided crossings. 1850 1851\item In Example~8 we obtained MLWFs separately for the up- and down-spin 1852channels of bcc Fe without spin-orbit. The Wannier-interpolated DOS 1853was therefore automatically separated into minority and majority 1854contributions. For a spinor calculation we can still spin-decompose 1855the DOS, using 1856% 1857{\tt 1858\begin{quote} 1859dos = true 1860 1861spin\_decomp = true 1862 1863dos\_kmesh = 25 25 25 1864\end{quote} } 1865% 1866The data file {\tt Fe-dos.dat} created by {\tt postw90} contains the 1867up-spin and down-spin contributions in the third and fourth columns, 1868% 1869{\tt 1870\begin{quote} 1871myshell> gnuplot 1872 1873gnuplot> plot 'Fe-dos.dat' u (-\$3):(\$1-12.6285) w l,'Fe-dos.dat' u (\$4):(\$1-12.6285) w l 1874\end{quote} } 1875% 1876(You should replace 12.6285 with your value of the Fermi energy). An 1877alternative approach is to project the DOS onto the up-spin and 1878down-spin WFs separately. To find the DOS projected onto the up-spin 1879(odd-numbered) WFs replace {\tt spin\_decomp = true} with 1880% 1881{\tt 1882\begin{quote} 1883 dos\_project = 1,3,5,7,9,11,13,15,17 1884\end{quote} } 1885% 1886and re-run {\tt postw90}. This approach has the advantage that it does 1887not require the {\tt Fe.spn} file. 1888 1889 1890\end{itemize} 1891 1892%\cleardoublepage 1893 1894\sectiontitle{18: Iron -- Berry curvature, anomalous Hall 1895 conductivity and optical conductivity} 1896 1897Note: This example requires a recent version of the {\tt pw2wannier90} interface. 1898 1899\begin{itemize} 1900\item{Outline: \it{Calculate the Berry curvature, anomalous Hall 1901 conductivity, and (magneto)optical conductivity of ferromagnetic 1902 bcc Fe with spin-orbit coupling. In preparation for this example 1903 it may be useful to read Ref.~\cite{yao-prl04} and Ch.~11 of the 1904 User Guide.}} 1905\item{Directory: {\tt examples/example18/}} 1906\item{Input files} 1907\begin{itemize} 1908\item{ {\tt Fe.scf} {\it The \pwscf\ input file for ground state 1909 calculation}} 1910\item{ {\tt Fe.nscf} {\it The \pwscf\ input file to obtain Bloch 1911 states on a uniform grid}} 1912\item{ {\tt Fe.pw2wan} {\it The input file for {\tt pw2wannier90}}} 1913\item{ {\tt Fe.win} {\it The {\tt wannier90} and {\tt postw90} input file}} 1914\end{itemize} 1915\end{itemize} 1916 1917The sequence of steps below is the same of Example~17. If you have 1918already run that example, you can reuse the output files from steps 19191--5, and only step 6 must be carried out again using the new input 1920file {\tt Fe.win}. 1921 1922\begin{enumerate} 1923\item Run \pwscf\ to obtain the ground state of iron\\ 1924{\tt pw.x < Fe.scf > scf.out} 1925 1926\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 1927 grid\\ 1928{\tt pw.x < Fe.nscf > nscf.out} 1929 1930\item Run \wannier\ to generate a list of the required overlaps (written 1931 into the {\tt Fe.nnkp} file)\\ 1932{\tt wannier90.x -pp Fe} 1933 1934\item Run {\tt pw2wannier90} to compute the overlaps between Bloch 1935 states and the projections for the starting guess (written in the 1936 {\tt Si.mmn} and {\tt Si.amn} files)\\ 1937{\tt pw2wannier90.x < Fe.pw2wan > pw2wan.out} 1938 1939\item Run \wannier\ to compute the MLWFs\\ 1940{\tt wannier90.x Fe} 1941 1942\item Run \postw\ \\ 1943 {\tt postw90.x Fe} (serial execution)\\ 1944 {\tt mpirun -np 8 postw90.x Fe} (example of parallel execution with 1945 8 MPI processes) 1946 1947\end{enumerate} 1948 1949\subsection*{Berry curvature plots} 1950 1951The Berry curvature $\Omega_{\alpha\beta}({\bf k})$ of the occupied 1952states is defined in Eq.~(11.18) of the User Guide. The following 1953lines in {\tt Fe.win} are used to calculate the energy bands and the 1954Berry curvature (in bohr$^2$) along high-symmetry lines in $k$-space. 1955{\tt 1956\begin{quote} 1957fermi\_energy = [insert your value here] 1958 1959berry\_curv\_unit = bohr2 1960 1961kpath = true 1962 1963kpath\_task = bands+curv 1964 1965kpath\_bands\_colour = spin 1966 1967kpath\_num\_points = 1000 1968\end{quote} } 1969 1970%The path specification in {\tt Fe.win} is the same as in Example~17, and 1971 1972 1973After executing {\tt postw90}, plot the Berry curvature component 1974$\Omega_z({\bf k})=\Omega_{xy}({\bf k})$ along the magnetization 1975direction using the script generated at runtime, 1976% 1977{\tt 1978\begin{quote} 1979myshell> python Fe-bands+curv\_z.py 1980\end{quote} } 1981 1982and compare with Fig.~2 of Ref.~\cite{yao-prl04}. 1983 1984In Example~17 we plotted the Fermi lines on the (010) plane $k_y=0$. 1985To combine them with a heatmap plot of (minus) the Berry curvature set 1986{\tt kpath = false}, uncomment the following lines in {\tt Fe.win}, 1987\smallskip {\tt 1988\begin{quote} 1989kslice = true 1990 1991kslice\_task = curv+fermi\_lines 1992 1993kslice\_corner = 0.0 0.0 0.0 1994 1995kslice\_b1 = 0.5 -0.5 -0.5 1996 1997kslice\_b2 = 0.5 0.5 0.5 1998 1999kslice\_2dkmesh = 200 200 2000\end{quote} } 2001% 2002re-run {\tt postw90}, and issue 2003% 2004{\tt 2005\begin{quote} 2006myshell> python Fe-kslice-curv\_z+fermi\_lines.py 2007\end{quote} } 2008\medskip 2009% 2010Compare with Fig.~3 in Ref.~\cite{yao-prl04}. Note how the Berry 2011curvature ``hot-spots'' tend to occur near spin-orbit-induced avoided 2012crossings (the Fermi lines with and without spin-orbit were generated 2013in Example~17). 2014 2015\subsection*{Anomalous Hall conductivity} 2016 2017The intrinsic anomalous Hall conductivity (AHC) is proportional to the 2018BZ integral of the Berry curvature. In bcc Fe with the magnetization 2019along $\hat{\bf z}$, the only nonzero components are 2020$\sigma_{xy}=-\sigma_{yx}$. To evaluate the AHC using a $25\times 202125\times 25$ $k$-point mesh, set {\tt kslice = false}, uncomment the 2022following lines in {\tt Fe.win}, {\tt 2023\begin{quote} 2024berry = true 2025 2026berry\_task = ahc 2027 2028berry\_kmesh = 25 25 25 2029 2030\end{quote} } and re-run {\tt postw90}. The AHC is written in the 2031output file {\tt Fe.wpout} in vector form. For bcc Fe with the 2032magnetization along [001], only the $z$-component $\sigma_{xy}$ is 2033nonzero. 2034 2035As a result of the strong and rapid variations of the Berry curvature 2036across the BZ, the AHC converges rather slowly with $k$-point 2037sampling, and a $25\times 25\times 25$ does not yield a well-converged 2038value. 2039 2040 \begin{itemize} 2041 2042 \item[{\bf --}] Increase the BZ mesh density by changing {\tt 2043 berry\_kmesh}. 2044 2045 \item[{\bf --}] To accelerate the convergence, adaptively refine the 2046 mesh around spikes in the Berry curvature, by adding to {\tt 2047 Fe.win} the lines \smallskip {\tt 2048 \begin{quote} 2049 berry\_curv\_adpt\_kmesh = 5 2050 2051 berry\_curv\_adpt\_kmesh\_thresh = 100.0 2052 \end{quote} } 2053 2054 \end{itemize} 2055 2056 This adds a $5\times 5\times 5$ fine mesh around those points 2057 where $\vert{\bm \Omega}({\bf k})\vert$ exceeds 100~bohr$^2$. The 2058 percentage of points triggering adaptive refinement is reported in 2059 {\tt Fe.wpout}. 2060 2061 Compare the converged AHC value with those obtained in 2062 Refs.~\cite{wang-prb06} and~\cite{yao-prl04}. 2063 2064 The Wannier-interpolation formula for the Berry curvature 2065 comprises three terms, denoted $D$-$D$, $D$-$\overline{A}$, and 2066 $\overline{\Omega}$ in Ref.~\cite{wang-prb06}, and $J2$, $J1$, and 2067 $J0$ in Ref.~\cite{lopez-prb12}. To report in {\tt Fe.wpout} the 2068 decomposition of the total AHC into these three terms, set {\tt 2069 iprint} (verbosity level) to a value larger than one in {\tt 2070 Fe.win}. 2071 2072\subsection*{Optical conductivity} 2073 2074The optical conductivity tensor of bcc Fe with magnetization along 2075$\hat{\bf z}$ has the form 2076% 2077$$ 2078\bm{\sigma}=\bm{\sigma}^{\rm S}+\bm{\sigma}^{\rm A}= 2079\left( 2080\begin{array}{ccc} 2081\sigma_{xx} & 0 & 0\\ 20820 & \sigma_{xx} & 0\\ 20830 & 0 & \sigma_{zz} 2084\end{array} 2085\right)+ 2086\left( 2087\begin{array}{ccc} 20880 & \sigma_{xy} & 0 \\ 2089-\sigma_{xy} & 0 & 0\\ 20900 & 0 & 0 2091\end{array} 2092\right)$$ 2093% 2094where ``S'' and ``A'' stand for the symmetric and antisymmetric parts 2095and $\sigma_{xx}=\sigma_{yy}\not=\sigma_{zz}$. The dc AHC calculated 2096earlier corresponds to $\sigma_{xy}$ in the limit $\omega\rightarrow 20970$. At finite frequency $\sigma_{xy}=-\sigma_{yx}$ acquires an 2098imaginary part which describes magnetic circular dichoism (MCD). 2099 2100To compute the complex optical conductivity for $\hbar\omega$ 2101up to 7~eV, replace 2102{\tt 2103\begin{quote} 2104berry\_task = ahc 2105\end{quote} } 2106% 2107with 2108% 2109{\tt 2110\begin{quote} 2111berry\_task = kubo 2112\end{quote} } 2113% 2114add the line 2115% 2116{\tt 2117\begin{quote} 2118kubo\_freq\_max = 7.0 2119\end{quote} } 2120% 2121and re-run {\tt postw90}. Reasonably converged spectra can be 2122obtained with a $125\times 125\times 125$ $k$-point mesh. Let us first 2123plot the ac AHC in S/cm, as in the lower panel of Fig.~5 in 2124Ref.~\cite{yao-prl04}, {\tt 2125\begin{quote} 2126myshell> gnuplot 2127 2128gnuplot> plot `Fe-kubo\_A\_xy.dat' u 1:2 w l 2129\end{quote} } 2130 2131Comapare the $\omega\rightarrow 0$ limit with the result obtained 2132earlier by integrating the Berry curvature.\footnote{The calculation 2133 of the AHC using {\tt berry\_task = kubo} involves a truncation of 2134 the sum over empty states in the Kubo-Greenwood formula: see 2135 description of the keyword {\tt kubo\_eigval\_max} in the User 2136 Guide. As discussed around Eq.~(11.17) of the User Guide, no 2137 truncation is done with {\tt berry\_task = ahc}.} 2138 2139 2140Next we plot the MCD spectrum. Following Ref.~\cite{yao-prl04}, we 2141plot ${\rm Im}[\omega\sigma_{xy}(\hbar\omega)]$, in units of 2142$10^{29}$~sec$^{-2}$. The needed conversion factor is $9\times 214310^{-18}\times e/\hbar\simeq 0.0137$ ($e$ and $\hbar$ in SI units), 2144{\tt 2145\begin{quote} 2146gnuplot> set yrange[-5:15] 2147 2148gnuplot> plot `Fe-kubo\_A\_xy.dat' u 1:(\$1)*(\$3)*0.0137 w l 2149\end{quote} } 2150 2151\subsection*{Further ideas} 2152 2153\begin{itemize} 2154 2155\item Recompute the AHC and optical spectra of bcc Fe using projected 2156 $s$, $p$, and $d$-type Wannier functions instead of the hybridrized 2157 MLWFs (see Example~8), and compare the results. 2158 2159\item A crude way to model the influence of heterovalent alloying on 2160 the AHC is to assume that its only effect is to donate or deplete 2161 electrons, i.e., to shift the Fermi level of the pure 2162 crystal~\cite{yao-prb07}. 2163 2164Recalculate the AHC of bcc Fe for a range of Fermi energies within 2165$\pm 0.5$~eV of the true Fermi level. This calculation can be 2166streamlined by replacing in {\tt Fe.win} {\tt 2167\begin{quote} 2168fermi\_energy = [insert your value here] 2169\end{quote} } 2170% 2171with 2172% 2173{\tt 2174\begin{quote} 2175fermi\_energy\_min = [insert here your value minus 0.5] 2176 2177fermi\_energy\_max = [insert here your value plus 0.5] 2178\end{quote} } 2179% 2180Use a sufficiently dense BZ mesh with adaptive refinement. To plot 2181$\sigma_{xy}$ versus $\varepsilon_F$, issue 2182% 2183{\tt 2184\begin{quote} 2185myshell> gnuplot 2186 2187gnuplot> plot `Fe-ahc-fermiscan.dat' u 1:4 w lp 2188\end{quote} } 2189 2190\end{itemize} 2191 2192 2193%\cleardoublepage 2194 2195\sectiontitle{19: Iron -- Orbital magnetization} 2196 2197Note: This example requires a recent version of the {\tt pw2wannier90} interface. 2198 2199\begin{itemize} 2200\item{Outline: \it{Calculate the orbital magnetization of 2201 ferromagnetic bcc Fe by Wannier interpolation.}} 2202\item{Directory: {\tt examples/example19/}} 2203\item{Input files} 2204\begin{itemize} 2205\item{ {\tt Fe.scf} {\it The \pwscf\ input file for ground state 2206 calculation}} 2207\item{ {\tt Fe.nscf} {\it The \pwscf\ input file to obtain Bloch 2208 states on a uniform grid}} 2209\item{ {\tt Fe.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2210\item{ {\tt Fe.win} {\it The {\tt wannier90} and {\tt postw90} input 2211 file}} 2212\end{itemize} 2213\end{itemize} 2214 2215The sequence of steps below is the same of Examples~17 and 18. If you 2216have already run one of those examples, you can reuse the output files 2217from steps 1--3 and 5. Steps~4 and 6 should be carried out again using 2218the new input files {\tt Fe.pw2wan} and {\tt Fe.win}. 2219 2220\begin{enumerate} 2221\item Run \pwscf\ to obtain the ground state of iron\\ 2222{\tt pw.x < Fe.scf > scf.out} 2223 2224\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 2225 grid\\ 2226{\tt pw.x < Fe.nscf > nscf.out} 2227 2228\item Run \wannier\ to generate a list of the required overlaps (written 2229 into the {\tt Fe.nnkp} file).\\ 2230{\tt wannier90.x -pp Fe} 2231 2232 2233\item Run {\tt pw2wannier90} to compute: 2234 \begin{itemize} 2235 2236 \item[{\bf --}] The overlaps $\langle u_{n{\bf k}}\vert u_{m{\bf k}+{\bf 2237 b}}\rangle$ (written in the {\tt Fe.mmn} file) 2238 2239 \item[{\bf --}] The projections for the starting guess (written in the {\tt 2240 Fe.amn} file) 2241 2242 \item[{\bf --}] The matrix elements $\langle u_{n{\bf k}+{\bf b}_1}\vert 2243 H_{\bf k}\vert u_{m{\bf k}+{\bf b}_2}\rangle$ (written in the 2244 {\tt Fe.uHu} file) 2245 2246% \item[{\bf --}] The spin matrix elements $\langle \psi_{n{\bf 2247% k}}\vert \sigma_i\vert \psi_{m{\bf k}}\rangle$ (written in the 2248% {\tt Fe.spn} file) 2249 2250 \end{itemize} 2251{\tt pw2wannier90.x < Fe.pw2wan > pw2wan.out} 2252 2253\item Run \wannier\ to compute the MLWFs.\\ 2254{\tt wannier90.x Fe} 2255 2256\item Run \postw\ to compute the orbital magnetization.\\ 2257 {\tt postw90.x Fe} (serial execution)\\ 2258 {\tt mpirun -np 8 postw90.x Fe} (example of parallel execution with 2259 8 MPI processes) 2260 2261 2262\end{enumerate} 2263 2264The orbital magnetization is computed as the BZ integral of the 2265quantity ${\bf M}^{\rm orb}({\bf k})$ defined in Eq.~(11.20) of the 2266User Guide. The relevant lines in {\tt Fe.win} are 2267% 2268{\tt 2269\begin{quote} 2270berry = true 2271 2272berry\_task = morb 2273 2274berry\_kmesh = 25 25 25 2275 2276fermi\_energy = [insert your value here] 2277\end{quote} } 2278% 2279After running {\tt postw90}, compare the value of the orbital 2280magnetization reported in {\tt Fe.wpout} with the spin magnetization 2281in {\tt scf.out}. Set {\tt iprint = 2} to report the decomposition of 2282${\bf M}^{\rm orb}$ into the terms $J0$, $J1$, and $J2$ defined in 2283Ref.~\cite{lopez-prb12}. 2284%\footnote{Alternatively, the 2285%spin magnetization can be computed by {\tt postw90}. This requires adding the 2286%instruction {\tt spin\_moment = T} to {\tt Fe.win}, assuming that {\tt pw2wannier90} 2287%was executed with {\tt write\_spn = .true.} in {\tt Fe.pw2wan}.} 2288%Both point along $\hat{\bf 2289% z}$, but the orbital contribution is much smaller (it is strongly 2290%quenched by the crystal field). 2291 2292To plot $M_z^{\rm orb}({\bf k})$ along high-symmetry lines set {\tt 2293 berry = false} and uncomment in {\tt Fe.win} the block of 2294instructions containing {\tt 2295\begin{quote} 2296kpath = true 2297 2298kpath\_task = bands+morb 2299\end{quote} 2300} 2301After running {\tt postw90}, issue 2302% 2303{\tt 2304\begin{quote} 2305myshell> python Fe-bands+morb\_z.py 2306\end{quote} } 2307% 2308Compare with Fig.~2 of Ref.~\cite{lopez-prb12}, bearing in mind the 2309factor of $-1/2$ difference in the definition of ${\bf M}^{\rm 2310 orb}({\bf k})$ (see Ch.~11 in the User Guide). 2311 2312To plot $M_z^{\rm orb}({\bf k})$ together with the Fermi contours on the 2313(010) BZ plane set {\tt kpath = false}, uncomment in {\tt 2314 Fe.win} the block of instructions containing 2315{\tt 2316\begin{quote} 2317kslice = true 2318 2319kslice\_task = morb+fermi\_lines 2320\end{quote} 2321} 2322re-run {\tt postw90}, and issue 2323{\tt 2324\begin{quote} 2325myshell> python Fe-kslice-morb\_z+fermi\_lines.py 2326\end{quote} } 2327% 2328$M_z^{\rm orb}({\bf k})$ is much more evenly distributed in $k$-space 2329than the Berry curvature (see Example 18). As a result, the integrated 2330orbital magnetization converges more rapidly with the BZ sampling. 2331 2332 2333%%\sectiontitle{20: Disentanglement only in small (spherical) regions of $k$ space} 2334%% 2335\sectiontitle{20: Disentanglement restricted inside spherical regions of $k$ space} 2336\subsection*{LaVO$_3$} 2337\begin{itemize} 2338\item{Outline: \it{Obtain disentangled MLWFs for strained LaVO$_3$.}} 2339\item{Directory: {\tt examples/example20/}} 2340\item{Input Files} 2341\begin{itemize} 2342\item{ {\tt LaVO3.scf} {\it The \pwscf\ input file for ground state 2343 calculation}} 2344\item{ {\tt LaVO3.nscf} {\it The \pwscf\ input file to obtain Bloch 2345 states on a uniform grid}} 2346\item{ {\tt LaV03.pw2wan} {\it Input file for {\tt pw2wannier90}}} 2347\item{ {\tt LaVO3.win} {\it The {\tt wannier90} input file}} 2348\end{itemize} 2349 2350\end{itemize} 2351 2352\begin{enumerate} 2353\item Run \pwscf\ to obtain the ground state of LaVO$_3$.\\ 2354{\tt pw.x < LaVO3.scf > scf.out} 2355 2356\item Run \pwscf\ to obtain the Bloch states on a uniform k-point 2357 grid.\\ 2358{\tt pw.x < LaVO3.nscf > nscf.out} 2359 2360\item Run \wannier\ to generate a list of the required overlaps (written 2361 into the {\tt LaVO3.nnkp} file).\\ 2362{\tt wannier90.x -pp LaVO3} 2363 2364\item Run {\tt pw2wannier90} to compute the overlap between Bloch 2365 states and the projections for the starting guess (written in the 2366 {\tt LaVO3.mmn} and {\tt LaVO3.amn} files).\\ 2367{\tt pw2wannier90.x < LaVO3.pw2wan > pw2wan.out} 2368 2369\item Run \wannier\ to compute the MLWFs.\\ 2370{\tt wannier90.x LaVO3} 2371 2372\end{enumerate} 2373 2374Inspect the output file {\tt LaVO3.wout}. In the initial summary, you 2375will see that the disentanglement was performed only within one sphere 2376of radius 0.2 arount the point $A=(0.5, 0.5, 0.5)$ in reciprocal space: 2377\begin{verbatim} 2378 | Number of spheres in k-space : 1 | 2379 | center n. 1 : 0.500 0.500 0.500, radius = 0.200 | 2380\end{verbatim} 2381 2382Compare the band structure that Wannier90 produced with the one 2383obtained using Quantum ESPRESSO. You should get something similar to Fig.~\ref{fig:lavo3}. 2384%% 2385Notice how the $t_{2g}$-bands are entangled with other bands at $A$ and the Wannier-interpolated band structure deviates from the Bloch bands only in a small region around that $k$-point. 2386It is important to keep in mind that all symmetry equivalent $k$-points within the first Brillouin zone must be written explicitly in the list of sphere centers. 2387For instance, the $A$ point in the simple tetragonal lattice of this example is non-degenerate, while the $X$ point has degeneracy two, hence one must specify both $(1/2,0,0)$ and $(0,1/2,0)$ (see the SrMnO$_3$ example here below). 2388 2389\begin{figure}[h] 2390\begin{center} 2391\includegraphics[width=10cm]{LaVO3} 2392\caption{Band structure of epitaxially-strained (tetragonal) LaVO$_3$. Black: Bloch bands; red circles: Wannier-interpolated band structure. The disentanglement was performed only for $k$-points within a sphere of radius 0.2 \AA$^{-1}$ centered in $A$.} 2393\label{fig:lavo3} 2394\end{center} 2395\end{figure} 2396 2397\subsection*{Further ideas} 2398 2399\begin{itemize} 2400%\item Try to disentangle the Wannier functions without disentanglement, and with standard disentanglement without using the spheres. You will notice that without disentanglement Wannier functions will not converge. With standard disentanglement (without spheres), also bands away from the A point in $k-$space are perturbed, even in regions of $k-$space where there is no entanglement. 2401\item Try to obtain the Wannier functions using the standard disentanglement procedure (without spheres, \verb+dis_spheres_num = 0+). You will notice that the Wannier-interpolated band structure now shows deviations also in regions of $k$-space far away from $A$, where disentanglement is actually not necessary. If you disable the disentanglement completely, instead, the Wannierisation procedure does not converge. 2402 2403\item In order to illustrate all possible cases, it is instructive to apply this method to SrMnO$_3$, where the $t_{2g}$ bands are entangled with the above-lying $e_g$ bands, and also with the deeper O-$2p$ states. 2404 In the SrMnO$_3$ subfolder, you can find input files for building three different sets of Wannier functions: only $t_{2g}$ states, only $e_g$ states, or all V-$3d$-derived states ($t_{2g} + e_g$). In each case one needs to specify different disentanglement spheres, according to which region(s) in $k$-space show entanglement of the targeted bands. 2405 Also the index \verb+dis_sphere_first_wan+ needs to be adapted to the new disentanglement window, which here contains also states below the lowest-lying Wannier function (at variance with the LaVO$_3$ case). 2406\end{itemize} 2407 2408%\cleardoublepage 2409 2410\sectiontitle{21: Gallium Arsenide -- Symmetry-adapted Wannier functions} 2411 2412Note: This example requires a recent version of the {\tt pw2wannier90} interface. 2413 2414\begin{itemize} 2415\item{Outline: \it{Obtain symmetry-adapted Wannier functions out of four valence bands of GaAs. 2416For the theoretical background of the symmetry-adapted Wannier functions, see R. Sakuma, Phys. Rev. B {\bf 87}, 235109 (2013).}} 2417\item{Directory: {\tt examples/example21/atom\_centered\_As\_sp/} \\ 2418\phantom{Directory: }{\tt examples/example21/atom\_centered\_Ga\_p/} \\ 2419\phantom{Directory: }{\tt examples/example21/atom\_centered\_Ga\_s/} \\ 2420\phantom{Directory: }{\tt examples/example21/atom\_centered\_Ga\_sp/} \\ 2421\phantom{Directory: }{\tt examples/example21/bond\_centered/} 2422} 2423\item{Input Files} 2424\begin{itemize} 2425\item{ {\tt GaAs.scf} {\it The \pwscf\ input file for ground state 2426 calculation}} 2427\item{ {\tt GaAs.nscf} {\it The \pwscf\ input file to obtain Bloch 2428 states on a uniform grid}} 2429\item{ {\tt GaAs.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2430\item{ {\tt GaAs.win} {\it The {\tt wannier90} input file}} 2431\end{itemize} 2432\end{itemize} 2433 2434\begin{enumerate} 2435\item Run \pwscf\ to obtain the ground state of GaAs\\ 2436{\tt pw.x < GaAs.scf > scf.out} 2437 2438\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 2439{\tt pw.x < GaAs.nscf > nscf.out} 2440 2441\item Run \wannier\ to generate a list of the required overlaps (written 2442 into the {\tt GaAs.nnkp} file).\\ 2443{\tt wannier90.x -pp GaAs} 2444 2445\item Run {\tt pw2wannier90} to compute the overlap between Bloch 2446 states, the projections for the starting guess, and the symmetry information needed for symmetry-adapted mode (written in the 2447 {\tt GaAs.mmn}, {\tt GaAs.amn}, and {\tt GaAs.dmn} files, respectively).\\ 2448{\tt pw2wannier90.x < GaAs.pw2wan > pw2wan.out} 2449 2450\item Run \wannier\ to compute the MLWFs.\\ 2451{\tt wannier90.x GaAs} 2452\end{enumerate} 2453 2454Each directory creates different kind of symmetry-adapted Wannier function. 2455See more detail in {\tt examples/example21/README}. 2456%Compare the results with those of Sec. III.A in R. Sakuma, Phys. Rev. B {\bf 87}, 235109 (2013). 2457 2458 2459%\cleardoublepage 2460 2461\sectiontitle{22: Copper -- Symmetry-adapted Wannier functions} 2462 2463Note: This example requires a recent version of the {\tt pw2wannier90} interface. 2464 2465\begin{itemize} 2466\item{Outline: \it{Obtain symmetry-adapted Wannier functions for Cu. By symmetry-adapted mode, for example, we can make atomic centered $s$-like Wannier function, which is not possible in the usual procedure to create maximally localized Wannier functions. 2467For the theoretical background of the symmetry-adapted Wannier functions, see R. Sakuma, Phys. Rev. B {\bf 87}, 235109 (2013).}} 2468\item{Directory: {\tt examples/example22/s\_at\_0.00/} \\ 2469\phantom{Directory: }{\tt examples/example22/s\_at\_0.25/} \\ 2470\phantom{Directory: }{\tt examples/example22/s\_at\_0.50/} 2471} 2472\item{Input Files} 2473\begin{itemize} 2474\item{ {\tt Cu.scf} {\it The \pwscf\ input file for ground state 2475 calculation}} 2476\item{ {\tt Cu.nscf} {\it The \pwscf\ input file to obtain Bloch 2477 states on a uniform grid}} 2478\item{ {\tt Cu.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2479\item{ {\tt Cu.sym} {\it Used only in {\tt examples/example22/s\_at\_0.25/}. {\tt pw2wannier90} reads this file when {\tt ``read\_sym = .true.''} in {\tt Cu.pw2wan}. By default, {\tt ``read\_sym = .false.'' and {\tt Cu.sym} is the output of {\tt pw2wannier90}}, in which the symmetry operations employed in the calculation are written for reference. } } 2480\item{ {\tt Cu.win} {\it The {\tt wannier90} input file}} 2481\end{itemize} 2482\end{itemize} 2483 2484\begin{enumerate} 2485\item Run \pwscf\ to obtain the ground state of Cu\\ 2486{\tt pw.x < Cu.scf > scf.out} 2487 2488\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 2489{\tt pw.x < Cu.nscf > nscf.out} 2490 2491\item Run \wannier\ to generate a list of the required overlaps (written 2492 into the {\tt Cu.nnkp} file).\\ 2493{\tt wannier90.x -pp Cu} 2494 2495\item Run {\tt pw2wannier90} to compute the overlap between Bloch 2496 states, the projections for the starting guess, and the symmetry information needed for symmetry-adapted mode (written in the 2497 {\tt Cu.mmn}, {\tt Cu.amn}, and {\tt Cu.dmn} files, respectively).\\ 2498{\tt pw2wannier90.x < Cu.pw2wan > pw2wan.out} 2499 2500\item Run \wannier\ to compute the MLWFs.\\ 2501{\tt wannier90.x Cu} 2502\end{enumerate} 2503 2504Each directory creates $s$-like symmetry-adapted Wannier function centered at different position on top of atomic centered $d$-like Wannier functions. 2505See more detail in {\tt examples/example22/README}. 2506%Compare the results with those of Sec. III.B in R. Sakuma, Phys. Rev. B {\bf 87}, 235109 (2013). 2507 2508 2509\sectiontitle{23: Silicon -- $G_0W_0$ bands structure interpolation} 2510 2511Note: This example requires a recent version of the {\tt ypp} post-processing code of {\tt yambo}. 2512 2513\begin{itemize} 2514\item{Outline: \it{Interpolate the bands structure of silicon obtained from many-body perturbation theory at the $G_0W_0$ level. Using the {\tt yambo} code, the quasi-particle corrections (QP) are summed to Kohn-Sham eigenvalues, while the wavefunctions remain the same. }} 2515\item{Directory: {\tt examples/example23/}} 2516\item{Input Files} 2517\begin{itemize} 2518\item{ {\tt silicon.scf} {\it The \pwscf\ input file for the ground state 2519 calculation}} 2520\item{ {\tt silicon.nscf } {\it The \pwscf\ input file to obtain Bloch 2521 states on a uniform grid}} 2522\item{ {\tt silicon.gw.nscf } {\it The \pwscf\ input file to obtain Bloch 2523 states on a reduced grid with many empty bands}} 2524\item{ {\tt silicon.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2525\item{ {\tt silicon.win} {\it The {\tt wannier90} input file}} 2526\item{ {\tt silicon.gw.win} {\it The {\tt wannier90} input file} (for the $G_0W_0$ step)} 2527\item{ {\tt yambo.in} {\it The {\tt yambo} input file}} 2528\item{ {\tt ypp.in} {\it The {\tt ypp} input file}} 2529\end{itemize} 2530\end{itemize} 2531 2532\begin{enumerate} 2533\item Copy the input files from the {\tt INPUT directory} into a working directory (e.g. {\tt WORK}) 2534\item Run \pwscf\ to obtain the ground state charge of silicon \\ 2535{\tt pw.x < silicon.scf > scf.out} 2536\item Run \pwscf\ to obtain the Bloch states reduced grid. We use a 8x8x8 with many bands (many empty bands are needed to perform a $G_0W_0$ with {\tt yambo})\\ 2537{\tt pw.x < silicon.gw.nscf > nscf.gw.out} 2538\item Use the {\tt k\_mapper.py} utility to find the indexes of a 4x4x4 uniform grid into the 8x8x8 reduced grid \\ 2539{\tt ./k\_mapper.py 4 4 4 "../examples/example23/WORK/nscf.gw.out"}\\ 2540Use the output to complete the {\tt yambo.in} input file (you also need to specify the on how many bands you want to compute the QP corrections, here you can use all the bands from 1 to 14). Then, you should have obtained something like:\\ 2541 {\tt 25421| 1| 1|14| \\ 25433| 3| 1|14| \\ 25445| 5| 1|14| \\ 254513| 13| 1|14| \\ 2546...\tt} 2547\item Enter the {\tt si.save} directory and run {\tt p2y}. A {\tt SAVE} folder is created, you can move it up in the {\tt /WORK/} directory.\\ 2548\item Run a $G_0W_0$ calculation from the {\tt /WORK/} directory (remember, we are using a 8x8x8 grid but computing QP corrections only on a 4x4x4 grid)\\ 2549{\tt yambo } 2550 2551 2552\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 2553{\tt pw.x < silicon.nscf > nscf.out} 2554 2555\item Run \wannier\ to generate a list of the required overlaps (written 2556 into the {\tt silicon.nnkp} file).\\ 2557{\tt wannier90.x -pp silicon} 2558 2559\item Run {\tt pw2wannier90} to compute the overlap between Bloch 2560 states, the projections for the starting guess (written in the 2561 {\tt silicon.mmn} and {\tt silicon.amn} respectively).\\ 2562{\tt pw2wannier90.x < silicon.pw2wan > pw2wan.out} 2563 2564\item Run \wannier\ to compute the MLWFs.\\ 2565{\tt wannier90.x silicon}\\ 2566At this point, you should have obtained the interpolated valence bands for silicon at the DFT level. 2567\item Run a {\tt ypp} calculation (just type {\tt ypp})\\ 2568You should obtain a file {\tt silicon.gw.unsorted.eig} which contains the QP corrections on a uniform 4x4x4 grid. 2569\item Run the gw2wannier90.py script to reorder, align and correct all matrices and files using the QP corrections\\ 2570{\tt ../../../utility/gw2wannier90.py silicon mmn amn} 2571\item Run \wannier\ to compute the MLWFs.\\ 2572{\tt wannier90.x silicon.gw}\\ 2573At this point, you should have obtained the interpolated valence bands for silicon at the $G_0W_0$ level. 2574\end{enumerate} 2575After you completed the tutorial for the valence bands only, you can repeat the final steps to interpolate also some conduction bands using disentanglement (the code is already present as comments in the input files). 2576%\cleardoublepage 2577 2578 2579\sectiontitle{24: Tellurium -- gyrotropic effects} 2580 2581 2582\begin{itemize} 2583\item{Outline: {\it Calculate the gyrotropic effects in trigonal right-handed Te} Similar to the calculations of \cite{tsirkin-arxiv17}} 2584\item{Directory: {\tt examples/example24/}} 2585\item{Input files} 2586\begin{itemize} 2587\item{ {\tt Te.scf} {\it The \pwscf\ input file for ground state 2588 calculation}} 2589\item{ {\tt Te.nscf} {\it The \pwscf\ input file to obtain Bloch 2590 states on a uniform grid}} 2591\item{ {\tt Te.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2592\item{ {\tt Te.win} {\it The {\tt wannier90} input 2593 file}} 2594\end{itemize} 2595\end{itemize} 2596 2597To make things easy, the example treats Te without spin-orbit 2598 2599 2600 2601\begin{enumerate} 2602\item Run \pwscf\ to obtain the ground state of tellurium\\ 2603{\tt pw.x < Te.scf > scf.out} 2604 2605\item Run \pwscf\ to obtain the Bloch states on a uniform {\tt 3x3x4} k-point 2606 grid\\ 2607{\tt pw.x < Te.nscf > nscf.out} 2608 2609\item Run \wannier\ to generate a list of the required overlaps (written 2610 into the {\tt Te.nnkp} file).\\ 2611{\tt wannier90.x -pp Te} 2612 2613\item Run {\tt pw2wannier90} to compute: 2614 \begin{itemize} 2615 2616 \item[{\bf --}] The overlaps $\langle u_{n{\bf k}}\vert u_{m{\bf k}+{\bf 2617 b}}\rangle$ (written in the {\tt Te.mmn} file) 2618 2619 \item[{\bf --}] The projections for the starting guess (written in the {\tt 2620 Te.amn} file) 2621 2622 \item[{\bf --}] The matrix elements $\langle u_{n{\bf k}+{\bf b}_1}\vert 2623 H_{\bf k}\vert u_{m{\bf k}+{\bf b}_2}\rangle$ (written in the 2624 {\tt Te.uHu} file) 2625 2626 \item[{\bf --}] The spin matrix elements $\langle \psi_{n{\bf 2627 k}}\vert \sigma_i\vert \psi_{m{\bf k}}\rangle$ (would be written in the 2628 {\tt Te.spn} file, but only if spin-orbit is included, which is not the case for the present example) 2629 2630 \end{itemize} 2631{\tt pw2wannier90.x < Te.pw2wan > pw2wan.out} 2632 2633\item Run \wannier\ to compute the MLWFs.\\ 2634{\tt wannier90.x Te} 2635 2636\item Add the following lines to the {\tt wannier90.win} file:\\ 2637{\tt gyrotropic=true \\ 2638gyrotropic\_task=-C-dos-D0-Dw-K \\ 2639fermi\_energy\_step=0.0025\\ 2640fermi\_energy\_min=5.8\\ 2641fermi\_energy\_max=6.2\\ 2642gyrotropic\_freq\_step=0.0025\\ 2643gyrotropic\_freq\_min=0.0\\ 2644gyrotropic\_freq\_max=0.1\\ 2645gyrotropic\_smr\_fixed\_en\_width=0.01\\ 2646gyrotropic\_smr\_max\_arg=5\\ 2647gyrotropic\_degen\_thresh=0.001\\ 2648gyrotropic\_box\_b1=0.2 0.0 0.0\\ 2649gyrotropic\_box\_b2=0.0 0.2 0.0\\ 2650gyrotropic\_box\_b3=0.0 0.0 0.2\\ 2651gyrotropic\_box\_center=0.33333 0.33333 0.5\\ 2652gyrotropic\_kmesh=50 50 50 2653} 2654 2655 2656 2657\item Run \postw\ \\to compute the gyrotropic properties: tensors $D$, $\widetilde{D}$, $K$, $C$ (See the User Guide):.\\ 2658 {\tt postw90.x Te} (serial execution)\\ 2659 {\tt mpirun -np 8 postw90.x Te} (example of parallel execution with 2660 8 MPI processes) \\ 2661 2662 2663The integration in the $k$-space is limited to a small area around the H point. Thus it is valid only for Fermi levels near the band gap. 2664And one needs to multiply the results by 2, to account for the H' point. To integrate over the entire Brillouin zone, one needs to remove the 2665{\tt gyrotropic\_box\_$\ldots$} parameters 2666 2667\item Now change the above lines to \\ {\tt 2668gyrotropic=true\\ 2669gyrotropic\_task=-NOA\\ 2670fermi\_energy=5.95\\ 2671gyrotropic\_freq\_step=0.0025\\ 2672gyrotropic\_freq\_min=0.0\\ 2673gyrotropic\_freq\_max=0.3\\ 2674gyrotropic\_smr\_fixed\_en\_width=0.01\\ 2675gyrotropic\_smr\_max\_arg=5\\ 2676gyrotropic\_band\_list=4-9\\ 2677gyrotropic\_kmesh=50 50 50\\ 2678} 2679 2680and compute the interband natural optical activity\\ 2681 2682 {\tt postw90.x Te} (serial execution)\\ 2683 {\tt mpirun -np 8 postw90.x Te} (example of parallel execution with 2684 8 MPI processes) \\ 2685 2686 2687 2688 2689 2690\end{enumerate} 2691 2692\sectiontitle{25: Gallium Arsenide -- Nonlinear shift current} 2693 2694\begin{itemize} 2695 2696\item Outline: \textit{Calculate the nonlinear shift current of inversion asymmetric fcc Gallium Arsenide. In preparation for this example it may be useful to read Ref. 2697\cite{ibanez-azpiroz_ab_2018} } 2698 2699 2700 2701\item Directory: \verb|examples/example25/| 2702 2703\item Input files: 2704 2705\begin{itemize} 2706 2707\item[--] \verb|GaAs.scf| \textit{The {\tt PWSCF} input file for ground state calculation} 2708\item[--] \verb|GaAs.nscf| \textit{The {\tt PWSCF} input file to obtain Bloch states on a uniform grid} 2709\item[--] \verb|GaAs.pw2wan| \textit{The input file for} \verb|pw2wannier90| 2710\item[--] \verb|GaAs.win| \textit{The} \verb|wannier90| \textit{and} \verb|postw90| \textit{input file} 2711 2712 2713\end{itemize} 2714 2715 2716\begin{enumerate} 2717 2718\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide 2719 2720\verb|pw.x < GaAs.scf > scf.out| 2721 2722 2723\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide 2724 2725\verb|pw.x < GaAs.nscf > nscf.out| 2726 2727\item Run {\tt Wannier90} to generate a list of the required overlaps (written into the \verb|GaAs.nnkp| file) 2728 2729\verb|wannier90.x -pp GaAs| 2730 2731 2732\item Run {\tt pw2wannier90} to compute: 2733 2734\begin{itemize} 2735\item[--] The overlaps $\langle u_{n\bf{k}}|u_{n\bf{k+b}}\rangle$ between spinor 2736Bloch states (written in the \verb|GaAs.mmn| file) 2737\item[--] The projections for the starting guess (written in the \verb|GaAs.amn| file) 2738 2739\end{itemize} 2740 2741 2742\verb|pw2wannier90.x < GaAs.pw2wan > pw2wan.out| 2743 2744\item Run {\tt wannier90} to compute MLWFs 2745 2746\verb|wannier90.x GaAs| 2747 2748\item Run {\tt postw90} to compute nonlinear shift current 2749 2750\verb|postw90.x GaAs| (serial execution) 2751 2752\verb|mpirun -np 8 postw90.x GaAs| (example of parallel execution with 8 MPI processes) 2753 2754\end{enumerate} 2755 2756 2757 2758 2759\end{itemize} 2760 2761 2762\subsection*{Shift current $\sigma^{abc}$} 2763 2764The shift current tensor of GaAs has only one independent component that is finite, namely $\sigma^{xyz}$. 2765For its computation, set 2766\begin{verbatim} 2767berry = true 2768berry_task = sc 2769\end{verbatim} 2770Like the linear optical conductivity, the shift current is a frequency-dependent quantity. 2771The frequency window and step is controlled by \verb|kubo_freq_min|, \verb|kubo_freq_max| and 2772\verb|kubo_freq_step|, as explained in the users guide. 2773 2774The shift current requires an integral over the Brillouin zone. The interpolated k-mesh is controlled by \verb|berry_kmesh|, 2775which has been set to 2776\begin{verbatim} 2777berry_kmesh = 100 100 100 2778 \end{verbatim} 2779We also need to input the value of the Fermi level in eV: 2780\begin{verbatim} 2781fermi_energy = [insert your value here] 2782\end{verbatim} 2783 2784Due to the sum over intermediate states involved in the calculation of the shift current, 2785one needs to consider a small broadening parameter to avoid numerical problems due to possible degeneracies 2786(see parameter $\eta$ in Eq. (36) of Ref. \cite{ibanez-azpiroz_ab_2018} and related discussion). 2787This parameter is controlled by \verb|sc_eta|. It is normally found that values between 0.01 eV and 0.1 eV 2788yield an stable spectrum. The default value is set to $0.04$ eV. 2789 2790Finally, \verb|sc_phase_conv| controls the phase convention used for the Bloch sums. 2791\verb|sc_phase_conv=1| uses the so-called tight-binding convention, whereby the Wannier centres are included 2792into the phase, while \verb|sc_phase_conv=2| leaves the Wannier centres out of the phase. 2793These two possible conventions are explained in Ref. \cite{pythtb}. 2794Note that the overall shift-current spectrum does not depend on the chosen convention, 2795but the individual terms that compose it do. 2796 2797 2798On output, the program generates a set of 18 files named \verb|SEED-sc_***.dat|, 2799which correspond to the different tensor components of the shift current 2800(note that the 9 remaining components until totaling $3\times3\times3=27$ 2801can be obtained from the 18 outputed by taking into account that $\sigma^{abc}$ is 2802symmetric under $b\leftrightarrow c$ index exchange). 2803For plotting the only finite shift-current component of GaAs $\sigma^{xyz}$ (units of A/V$^{2}$) as in the upper 2804panel of Fig. 3 in Ref. \cite{ibanez-azpiroz_ab_2018}, 2805\begin{verbatim} 2806myshell> gnuplot 2807gnuplot> plot 'GaAs-sc_xyz.dat' u 1:2 w l 2808\end{verbatim} 2809 2810 2811\sectiontitle{26: Gallium Arsenide -- Selective localization and constrained centres} 2812 2813\begin{itemize} 2814 2815\item Outline: \textit{Application of the selectively localised Wannier function (SLWF) method to gallium arsenide (GaAs), following the example in Ref. \cite{Marianetti}, which is essential reading for this tutorial example.} 2816 2817 2818\item Directory: \verb|examples/example26/| 2819 2820 2821\item Input files: 2822 2823\begin{itemize} 2824 2825\item[--] \verb|GaAs.scf| \textit{The {\tt PWSCF} input file for ground state calculation} 2826\item[--] \verb|GaAs.nscf| \textit{The {\tt PWSCF} input file to obtain Bloch states on a uniform grid} 2827\item[--] \verb|GaAs.pw2wan| \textit{The input file for} \verb|pw2wannier90| 2828\item[--] \verb|GaAs.win| \textit{The} \verb|wannier90| \textit{and} \verb|postw90| \textit{input file} 2829 2830 2831\end{itemize} 2832 2833\begin{enumerate} 2834 2835\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide 2836 2837\verb|pw.x < GaAs.scf > scf.out| 2838 2839 2840\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide 2841 2842\verb|pw.x < GaAs.nscf > nscf.out| 2843 2844\item Run {\tt Wannier90} to generate a list of the required overlaps (written into the \verb|GaAs.nnkp| file) 2845 2846\verb|wannier90.x -pp GaAs| 2847 2848 2849\item Run {\tt pw2wannier90} to compute: 2850 2851\begin{itemize} 2852\item[--] The overlaps $\langle u_{n\bf{k}}|u_{n\bf{k+b}}\rangle$ between 2853Bloch states (written in the \verb|GaAs.mmn| file) 2854\item[--] The projections for the starting guess (written in the \verb|GaAs.amn| file) 2855 2856\end{itemize} 2857 2858 2859\verb|pw2wannier90.x < GaAs.pw2wan > pw2wan.out| 2860 2861\item Inspect the {\tt .win} file. 2862 2863\begin{itemize} 2864\item[--] Make sure you understand the new keywords corresponding to the selective localisation algorithm. 2865\item[--] Run {\tt wannier90} to compute the SLWFs, in this case using one objective Wannier function. 2866 2867\end{itemize} 2868 2869 2870\verb|wannier90.x GaAs| 2871 2872\end{enumerate} 2873 2874To constrain the centre of the SLWF you need to add \mbox{{\tt slwf\_constrain = true}} and \\ 2875\mbox{{\tt slwf\_lambda = 1}} to the input file and uncomment the \mbox{{\tt slwf\_centres}} block. This will add a penalty functional to the total spread, which will try to constrain the centre of the SLWF to be on the As atom (as explained in Ref.~\cite{Marianetti}, particularly from Eq.~24 to Eq.~35). 2876 2877Look at the value of the penalty functional, is this what you would expect at convergence? 2878Does the chosen value of the Lagrange multiplier {\tt slwf\_lambda} give a SLWF function centred on the As atom? 2879 2880Alternatively, you can modify the {\tt slwf\_centres} block to constrain the centre of the SLWF to be on the Ga atom. 2881Do you need a different value of {\tt slwf\_lambda} in this case to converge? 2882Take a look at the result in Vesta and explain what you see. Do these functions transform like the identity under the action of the $T_d$ group? 2883 2884 2885\end{itemize} 2886 2887\sectiontitle{27: Silicon -- Selected columns of density matrix algorithm for automated MLWFs} 2888 2889Note: This example requires a recent version of the {\tt pw2wannier90.x} post-processing code of {\tt Quantum ESPRESSO} (v6.4 or above). 2890 2891\begin{itemize} 2892\item{Outline: \it{For bulk crystalline Silicon, generate the $A_{mn}$ matrices via the selected columns of density matrix (SCDM) algorithm and the corresponding MLWFs for 1) Valence bands 2) Valence bands and 4 low-lying conduction bands 3) Conduction bands only. To better understand the input files and the results of these calculations, it is crucial that the Reader has familiarized with the concepts and methods explained in Ref.~\cite{LinLin-ArXiv2017}. More info on the keywords related to the SCDM method may be found in the user\_guide.}} 2893\item{Directory: {\tt examples/example27/}} 2894\item{Input Files: {\tt input\_files}, and in the three subfolders {\tt isolated, erfc} and {\tt gaussian}. 2895The {\tt input\_files} folder contains:} 2896\begin{itemize} 2897\item{ {\tt si.scf} {\it The \pwscf\ input file for the ground state 2898 calculation}} 2899\item{ {\tt si\_4bands.nscf } {\it The \pwscf\ input file to obtain Bloch 2900 states on a uniform grid for 4 bands.}} 2901\item{ {\tt si\_12bands.nscf } {\it The \pwscf\ input file to obtain Bloch 2902 states on a uniform grid for 12 bands.}} 2903\end{itemize} 2904\item{ Whereas the three subfolders {\tt isolated, erfc} and {\tt gaussian} contain the {\tt si.win} \wannier~ input files and {\tt si.pw2wan} {\tt pw2wannier90} input files each corresponding to one of the scenarios listed in the outline.} 2905\end{itemize} 2906 2907 2908\begin{itemize} 2909 \item[1]{Valence bands: In this case we will compute 4 localized WFs corresponding to the 4 valence bands of Silicon. These 4 bands constitute a manifold that is separated in energy from other bands. In this case the columns of the density matrix are already localized in real space and no extra parameter is required.} 2910 \begin{enumerate} 2911 \item Copy the input files {\tt si.scf} and {\tt si\_4bands.nscf} from the {\tt input\_files} directory into the {\tt isolated} folder 2912 \item Run \pwscf\ to obtain the ground state charge of bulk Silicon. \\ 2913 {\tt pw.x < si.scf > scf.out} 2914 2915 \item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid of 4x4x4 for 4 bands. \\ 2916 {\tt pw.x < si\_4bands.nscf > nscf.out} 2917 2918 \item Inspect the {\tt si.win} input file and make sure that the {\tt auto\_projections} flag is set to {\tt .true.}. Also, make sure that no projections block is present. 2919 \item Run \wannier\ to generate a list of the required overlaps and also info on the SCDM method (written 2920 into the {\tt si.nnkp} file).\\ 2921 {\tt wannier90.x -pp si} 2922 \item Inspect the {\tt si.nnkp} file and make sure you find the {\tt auto\_projections} block and that no projections have been written in the {\tt projections} block. 2923 2924 \item Inspect the {\tt .pw2wan} input file. You will find two new keywords, i.e. {\tt scdm\_proj} and {\tt scdm\_entanglement}. The former, will instruct {\tt pw2wannier90.x} to use the SCDM method when generating the $A_{mn}$ matrix. The latter, defines which formula to adopt for the function $f(\varepsilon_{n\mathbf{k}})$ (see \cite{LinLin-ArXiv2017} and point below). 2925 \item Run {\tt pw2wannier90} to compute the overlap between Bloch 2926 states and the projections via the SCDM method (written in the 2927 {\tt si.mmn} and {\tt si.amn} respectively).\\ 2928 {\tt pw2wannier90.x < si.pw2wan > pw2wan.out} 2929 2930 \item Run \wannier\ to compute the MLWFs.\\ 2931 {\tt wannier90.x si}\\ 2932 At this point, you should have obtained 4 Wannier functions and the interpolated valence bands for Silicon. Inspect the output file {\tt si.wout}. In particular, look at the geometric centres of each WF, do they lie at the centre of the Si-Si bond as for the MLWFs computed from user-defined initial $s$-like projections (see Example11)? 2933 Plot these WFs using Vesta. Do they show the $\sigma$ character one would expect from chemical arguments? 2934 \end{enumerate} 2935 \item[2]{Valence bands + conduction bands: In this case we will compute 8 localized WFs corresponding to the 4 valence bands and 4 low-lying conduction bands. Here, we don't have a separate manifold, since the conduction bands are entangled with other high-energy bands and the columns of the density matrix are not exponentially localized by construction. A modified density matrix is required in this case\cite{LinLin-ArXiv2017}, and it is defined as: $$P(\mathbf{r},\mathbf{r}') = \sum_{n,\mathbf{k}} \psi_{n\mathbf{k}}(\mathbf{r})f(\varepsilon_{n,\mathbf{k}})\psi_{n\mathbf{k}}^\ast(\mathbf{r}'),$$ 2936 where $\psi_{n\mathbf{k}}$ and $\varepsilon_{n,\mathbf{k}}$ are the energy eigestates and eigenvalues from the first-principle calculation respectively. The function $f(\varepsilon_{n,\mathbf{k}})$ contains two free parameters $\mu$ and $\sigma$ and is defined as a complementary error function: $$f(\varepsilon_{n,\mathbf{k}}) = \frac{1}{2}\mathrm{erfc}\left(\frac{\varepsilon_{n,\mathbf{k}} - \mu}{\sigma}\right).$$ } 2937 \begin{enumerate} 2938 \item Copy the input files {\tt si.scf} and {\tt si\_12bands.nscf} from the {\tt input\_files} folder into the {\tt erfc} folder 2939 \item Run \pwscf\ to obtain the ground state charge of bulk Silicon. \\ 2940 {\tt pw.x < si.scf > scf.out} 2941 2942 \item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid of 4x4x4 for 12 bands this time. \\ 2943 {\tt pw.x < si\_12bands.nscf > nscf.out} 2944 2945 \item Inspect the {\tt si.win} input file and make sure that the {\tt auto\_projections} flag is set to {\tt .true.}. Also, make sure that no projection block is present. 2946 \item Run \wannier\ to generate a list of the required overlaps and also info on the SCDM method (written 2947 into the {\tt si.nnkp} file).\\ 2948 {\tt wannier90.x -pp si} 2949 \item Inspect the {\tt si.nnkp} file and make sure you find the {\tt auto\_projections} block and that no projections have been written in the {\tt projections} block. 2950 2951 \item Inspect the {\tt .pw2wan} input file. You will find other two new keywords, i.e. {\tt scdm\_mu} and {\tt scdm\_sigma}. These are the values in eV of $\mu$ and $\sigma$ in $f(\varepsilon_{n,\mathbf{k}})$, respectively. 2952 \item Run {\tt pw2wannier90} to compute the overlap between Bloch 2953 states and the projections via the SCDM method (written in the 2954 {\tt si.mmn} and {\tt si.amn} respectively).\\ 2955 {\tt pw2wannier90.x < si.pw2wan > pw2wan.out} 2956 2957 \item Run \wannier\ to compute the MLWFs.\\ 2958 {\tt wannier90.x si}\\ 2959 At this point, you should have obtained 8 localized Wannier functions and the interpolated valence and conduction bands for Silicon. Again, compare the results for the geometric centres and the individual spreads with the ones from Example11. Is the final value of total spread bigger or smaller than the one from Example11? Look at the WFs with Vesta. Can you explain what you see? Where do the major lobes of the $sp3$-like WFs point in this case? 2960 \end{enumerate} 2961 \item[3]{Conduction bands only: In this case we will compute 4 localized WFs corresponding to the 4 low-lying conduction bands only. As for the previous point, we need to define a modified density matrix\cite{LinLin-ArXiv2017}. Since we are only interested in a subset of the conduction states, within a bounded energy region, a good choice for $f(\varepsilon_{n,\mathbf{k}})$ is: $$f(\varepsilon_{n,\mathbf{k}}) = \exp\left(-\frac{(\varepsilon_{n,\mathbf{k}} - \mu)^2}{\sigma^2}\right).$$} 2962 \begin{enumerate} 2963 \item Copy the input files {\tt si.scf} and {\tt si\_12bands.nscf} from the {\tt input\_files} directory into the {\tt gaussian} folder 2964 \item Run \pwscf\ to obtain the ground state charge of bulk Silicon. \\ 2965 {\tt pw.x < si.scf > scf.out} 2966 2967 \item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid of 4x4x4 for 12 bands this time. \\ 2968 {\tt pw.x < si\_12bands.nscf > nscf.out} 2969 2970 \item Inspect the {\tt si.win} input file and make sure that the {\tt auto\_projections} flag is set to {\tt .true.}. Also, make sure that no projections block is present. 2971 \item Run \wannier\ to generate a list of the required overlaps and also info on the SCDM method (written 2972 into the {\tt si.nnkp} file).\\ 2973 {\tt wannier90.x -pp si} 2974 \item Inspect the {\tt si.nnkp} file and make sure you find the {\tt auto\_projections} block and that no projections have been written in the {\tt projections} block. 2975 2976 \item Run {\tt pw2wannier90} to compute the overlap between Bloch 2977 states, the projections for the starting guess via the SCDM method (written in the 2978 {\tt si.mmn} and {\tt si.amn} respectively).\\ 2979 {\tt pw2wannier90.x < si.pw2wan > pw2wan.out} 2980 2981 \item Run \wannier\ to compute the MLWFs.\\ 2982 {\tt wannier90.x si}\\ 2983 At this point, you should have obtained 4 localized Wannier functions and the interpolated conduction bands for Silicon. From chemical intuition, we would expect these functions to be similar to anti-bonding orbitals of molecules with tetrahedral symmetry. Plot the WFs and check if this is confirmed. 2984 \end{enumerate} 2985\end{itemize} 2986 2987\sectiontitle{28: Diamond -- plotting of MLWFs using Gaussian cube format and VESTA} 2988\begin{itemize} 2989\item{Outline: \it{Obtain MLWFs for the valence bands of diamond and output them in Gaussian cube format}} 2990\item{Directory: {\tt examples/example28/}} 2991The input files for this examples are the same as the ones in example05 2992\item{Input Files} 2993\begin{itemize} 2994\item{ {\tt diamond.scf} {\it The \pwscf\ input file for ground state 2995 calculation}} 2996\item{ {\tt diamond.nscf} {\it The \pwscf\ input file to obtain Bloch 2997 states on a uniform grid}} 2998\item{ {\tt diamond.pw2wan} {\it The input file for {\tt pw2wannier90}}} 2999\item{ {\tt diamond.win} {\it The {\tt wannier90} input file}} 3000\end{itemize} 3001\end{itemize} 3002 3003\begin{enumerate} 3004\item Run \pwscf\ to obtain the ground state of diamond\\ 3005{\tt pw.x < diamond.scf > scf.out} 3006 3007\item Run \pwscf\ to obtain the Bloch states on a uniform k-point grid\\ 3008{\tt pw.x < diamond.nscf > nscf.out} 3009 3010\item Run \wannier\ to generate a list of the required overlaps (written 3011 into the {\tt diamond.nnkp} file).\\ 3012{\tt wannier90.x -pp diamond} 3013 3014\item Run {\tt pw2wannier90} to compute the overlap between Bloch 3015 states and the projections for the starting guess (written in the 3016 {\tt diamond.mmn} and {\tt diamond.amn} files).\\ 3017{\tt pw2wannier90.x < diamond.pw2wan > pw2wan.out} 3018 3019\item When the lattice vectors are non-orthogonal, not all the visualisation programs are capable to plot volumetric data in the Gaussian cube format. 3020One program that can read volumetric data for these systems is VESTA. 3021To instruct \wannier\ to output the MLWFs data in Gaussian cube format you need to add the following lines to the {\tt .win} file 3022\begin{verbatim} 3023wannier_plot = .true. 3024wannier_plot_supercell = 3 3025wannier_plot_format = cube 3026wannier_plot_mode = crystal 3027wannier_plot_radius = 2.5 3028wannier_plot_scale = 1.0 3029\end{verbatim} 3030Run \wannier\ to compute the MLWFs and output them in the Gaussian cube file.\\ 3031{\tt wannier90.x diamond} 3032 3033\item Plot the first MLWF with VESTA 3034 {\tt vesta diamond\_00001.cube} 3035\end{enumerate} 3036 3037Extra: Instead of using {\tt wannier\_plot\_mode = crystal} try to use the molecule mode as {\tt wannier\_plot\_mode = molecule} (see the user guide for the definition of this keyword). 3038Add the following line to the {\tt .win} file: 3039\begin{verbatim} 3040restart = plot 3041\end{verbatim} 3042and re-run \wannier. Use VESTA to plot the resulting MLWFs, do you see any difference from the {\tt crystal} mode case? Can you explain why? 3043Try to change the size of the supercell from 3 to 5, do you expect the results to be different? ({\it Hint:} When using the Gaussian cube format the code outputs the WF on a grid that is smaller than the super 3044unit-cell. The size of the grid is specified by {\tt wannier\_plot\_scale} and {\tt wannier\_plot\_radius}.) 3045 3046\sectiontitle{29: Platinum -- Spin Hall conductivity} 3047 3048\begin{itemize} 3049 \item{Outline: {\it Calculate spin Hall conductivity (SHC) and 3050 plot Berry curvature-like term 3051 of fcc Pt considering spin-orbit coupling. 3052 To gain a better understanding of this example, 3053 it is suggested to read Ref.~\cite{qiao-prb2018} for a detailed 3054 description of the theory and Ch.~12.5 of the User Guide.}} 3055 \item{Directory: {\tt examples/example29/}} 3056 \item{Input files} 3057 \begin{itemize} 3058 \item{ {\tt Pt.scf} {\it The \pwscf\ input file for ground state 3059 calculation}} 3060 \item{ {\tt Pt.nscf} {\it The \pwscf\ input file to obtain Bloch 3061 states on a uniform grid}} 3062 \item{ {\tt Pt.pw2wan} {\it The input file for {\tt pw2wannier90}}} 3063 \item{ {\tt Pt.win} {\it The {\tt wannier90} and {\tt postw90} input file}} 3064 \end{itemize} 3065\end{itemize} 3066 3067\begin{enumerate} 3068 \item Run \pwscf\ to obtain the ground state of platinum\\ 3069 {\tt pw.x < Pt.scf > scf.out} 3070 3071 \item Run \pwscf\ to obtain the Bloch states on a uniform $k$-point 3072 grid\\ 3073 {\tt pw.x < Pt.nscf > nscf.out} 3074 3075 \item Run \wannier\ to generate a list of the required overlaps (written 3076 into the {\tt Pt.nnkp} file)\\ 3077 {\tt wannier90.x -pp Pt} 3078 3079 \item Run {\tt pw2wannier90} to compute the overlaps between Bloch 3080 states and the projections for the starting guess (written in the 3081 {\tt Pt.mmn} and {\tt Pt.amn} files)\\ 3082 {\tt pw2wannier90.x < Pt.pw2wan > pw2wan.out} 3083 3084 \item Run \wannier\ to compute the MLWFs\\ 3085 {\tt wannier90.x Pt} 3086 3087 \item Run \postw\ \\ 3088 {\tt postw90.x Pt} (serial execution)\\ 3089 {\tt mpirun -np 8 postw90.x Pt} (example of parallel execution with 3090 8 MPI processes) 3091 3092\end{enumerate} 3093 3094\subsection*{Spin Hall conductivity} 3095 3096The intrinsic spin Hall conductivity 3097$\sigma_{\alpha\beta}^{\text{spin}\gamma}$ is proportional to the 3098BZ integral of the Berry curvature-like term. 3099To evaluate the SHC using a $25\times 310025\times 25$ $k$-point mesh, set the 3101following lines in {\tt Pt.win}, {\tt 3102 \begin{quote} 3103 berry = true 3104 3105 berry\_task = shc 3106 3107 berry\_kmesh = 25 25 25 3108 3109\end{quote} } 3110When calculating SHC, adaptive smearing can be used 3111by commenting the following two lines, 3112{\tt 3113 \begin{quote} 3114 3115 \#kubo\_adpt\_smr = false 3116 3117 \#kubo\_smr\_fixed\_en\_width = 1 3118\end{quote} } 3119Then set the Fermi energy $\varepsilon_F$ to a specific value 3120{\tt 3121 \begin{quote} 3122 fermi\_energy = [insert your value here] 3123\end{quote} } 3124or invoke Fermi energy scan by setting 3125{\tt 3126 \begin{quote} 3127 fermi\_energy\_min = [insert here your lower range] 3128 3129 fermi\_energy\_max = [insert here your upper range] 3130 3131 fermi\_energy\_step = [insert here your step] 3132 3133\end{quote} } 3134and re-run {\tt postw90}. 3135The SHC is written in the 3136output file {\tt Pt-shc-fermiscan.dat}. 3137If only {\tt fermi\_energy} is set, the output file will contain 3138SHC at this specific energy; if a list of Fermi energies are set, 3139the output file will contain SHC calculated at each 3140energy point in the list: we call this the ``Fermi energy scan'' of SHC. 3141 3142To plot the Fermi energy scan of SHC 3143$\sigma_{xy}^{\text{spin}z}$ versus $\varepsilon_F$, issue 3144% 3145{\tt 3146 \begin{quote} 3147 myshell> gnuplot 3148 3149 gnuplot> plot `Pt-shc-fermiscan.dat' u 2:3 w lp 3150\end{quote} } 3151 3152As a result of the strong and rapid variations of the 3153Berry curvature-like term 3154across the BZ, the SHC converges rather slowly with $k$-point 3155sampling, and a $25\times 25\times 25$ kmesh does not yield a well-converged value. 3156 3157\begin{itemize} 3158 3159 \item[{\bf --}] Increase the kmesh density by changing {\tt 3160 berry\_kmesh}. 3161 3162 \item[{\bf --}] To accelerate the convergence, adaptively refine the 3163 kmesh around spikes in the Berry curvature-like term, by adding to {\tt 3164 Pt.win} the lines \smallskip {\tt 3165 \begin{quote} 3166 berry\_curv\_adpt\_kmesh = 5 3167 3168 berry\_curv\_adpt\_kmesh\_thresh = 100.0 3169 \end{quote} } 3170 This adds a $5\times 5\times 5$ fine mesh around those points 3171 where $\vert{\Omega_{\alpha\beta}^{\text{spin}\gamma}}({\bm k})\vert$ exceeds 100~{\tt berry\_curv\_unit}. The 3172 percentage of points triggering adaptive refinement is reported in 3173 {\tt Pt.wpout}. 3174 3175\end{itemize} 3176 3177 3178 3179Compare the converged SHC value with those obtained in 3180Refs.~\cite{qiao-prb2018} and~\cite{guo-prl2008}. 3181 3182Note some rough estimations of computation progress and time are reported in {\tt Pt.wpout} 3183(see the SHC part of the Solution Booklet). These may be helpful if the computation time is very long. 3184 3185\subsection*{Notes} 3186% 3187\begin{itemize} 3188 \item Since the Kubo formula of SHC involves unoccupied bands, 3189 we need to include some unoccupied bands and construct more MLWF. 3190 Thus the following parameters should be increased accordingly:\smallskip {\tt 3191 \begin{quote} 3192 dis\_froz\_max 3193 3194 dis\_win\_max 3195 3196 projections 3197 \end{quote} } 3198 \item Normally we calculate the SHC 3199 $\sigma_{xy}^{\text{spin}z}$, 3200 i.e. $\alpha = x, \beta = y, \gamma = z$. To calculate other components, the following parameters can be set as 3201 {\tt 1, 2, 3} \smallskip {\tt 3202 \begin{quote} 3203 shc\_alpha = [insert here the $\alpha$ direction] 3204 3205 shc\_beta = [insert here the $\beta$ direction] 3206 3207 shc\_gamma = [insert here the $\gamma$ direction] 3208 \end{quote} } 3209 with {\tt 1, 2, 3} standing for {\tt x, y, z} respectively. 3210\end{itemize} 3211 3212\subsection*{Berry curvature-like term plots} 3213 3214The band-projected Berry curvature-like term 3215$\Omega_{n,\alpha\beta}^{\text{spin} \gamma}({\bm k})$ 3216is defined in Eq.~(12.22) of the User Guide. The following 3217lines in {\tt Pt.win} are used to calculate the energy bands colored by the 3218band-projected Berry curvature-like term 3219$\Omega_{n,\alpha\beta}^{\text{spin} \gamma}({\bm k})$ (in \AA$^2$), 3220as well as the $k$-resolved Berry curvature-like term 3221$\Omega_{\alpha\beta}^{\text{spin} \gamma}({\bm k})$ 3222along high-symmetry lines in $k$-space, i.e. the {\tt kpath} plot. First comment the line {\tt berry = true} and then set 3223{\tt 3224 \begin{quote} 3225 3226 kpath = true 3227 3228 kpath\_task = bands+shc 3229 3230 kpath\_bands\_colour = shc 3231 3232 kpath\_num\_points = 400 3233 3234 kubo\_adpt\_smr = false 3235 3236 kubo\_smr\_fixed\_en\_width = 1 3237 3238 fermi\_energy = [insert your value here] 3239 3240 berry\_curv\_unit = ang2 3241\end{quote} } 3242 3243After executing {\tt postw90}, four files are generated: {\tt Pt-bands.dat}, {\tt Pt-path.kpt}, {\tt Pt-shc.dat} and {\tt Pt-bands+shc.py}. 3244Then plot the band-projected Berry curvature-like term 3245$\Omega_{n,\alpha\beta}^{\text{spin}\gamma}({\bm k})$ 3246using the script generated at runtime, 3247% 3248{\tt 3249 \begin{quote} 3250 myshell> python Pt-bands+shc.py 3251\end{quote} } 3252 3253and compare with Fig.~2 of Ref.~\cite{qiao-prb2018}. Note a 3254large fixed smearing of 1~eV is used to recover the result in Ref.~\cite{qiao-prb2018}. 3255You can adjust the {\tt kubo\_smr\_fixed\_en\_width} 3256as you like to draw a visually appealing figure. 3257A {\tt kpath} plot of 0.05 eV smearing is shown in the 3258Solution Booklet. 3259 3260Besides, you can set {\tt kpath\_task = shc} to only draw 3261$k$-resolved term 3262$\Omega_{\alpha\beta}^{\text{spin} \gamma}({\bm k})$ 3263(the lower panel of the figure), or set 3264{\tt kpath\_task = bands} and {\tt kpath\_bands\_colour = shc} 3265to only draw energy bands colored by the 3266band-projected term 3267$\Omega_{n,\alpha\beta}^{\text{spin} \gamma}({\bm k})$ 3268(the upper panel of the figure). 3269 3270Similar to that of AHC, we can get a heatmap plot of 3271the $k$-resolved Berry curvature-like term 3272$\Omega_{\alpha\beta}^{\text{spin}\gamma}({\bm k})$, 3273i.e. the {\tt kslice} plot. To move forward, 3274set {\tt kpath = false} and uncomment the following lines in {\tt Pt.win},\smallskip {\tt 3275 \begin{quote} 3276 kslice = true 3277 3278 kslice\_task = shc+fermi\_lines 3279 3280 kslice\_corner = 0.0 0.0 0.0 3281 3282 kslice\_b1 = 1.0 0.0 0.0 3283 3284 kslice\_b2 = 0.3535533905932738 1.0606601717798214 0.00 3285 3286 kslice\_2dkmesh = 200 200 3287\end{quote} } 3288% 3289Note the {\tt kslice\_b2} is actually $(\frac{\sqrt{2}}{4}, \frac{3\sqrt{2}}{4},0.0)$ which leads to a square slice in the BZ, 3290making it easier to plot in the generated {\tt python} script. 3291Re-run {\tt postw90}, and issue% 3292{\tt 3293 \begin{quote} 3294 myshell> python Pt-kslice-shc+fermi\_lines.py 3295\end{quote} } 3296\medskip 3297% 3298Compare the generated figure with Fig.~3 in Ref.~\cite{qiao-prb2018}, or 3299the Solution Booklet. 3300 3301\subsection*{Notes} 3302% 3303\begin{itemize} 3304\item Adaptive smearing depends on a uniform kmesh, so when running 3305{\tt kpath} and {\tt kslice} plots adaptive smearing should not be 3306used. A fixed smearing is needed to avoid 3307near zero number in the denominator of the Kubo formula, 3308Eq.~(12.22) in the User Guide. To add a fixed smearing of 33090.05 eV, add the following keywords in the {\tt Pt.win}, \smallskip {\tt 3310 \begin{quote} 3311 kubo\_adpt\_smr = .false. 3312 3313 kubo\_smr\_fixed\_en\_width = 0.05 3314\end{quote} } 3315\end{itemize} 3316 3317\subsection*{Input parameters for SHC} 3318Finally, we provide a complete list of input parameters that 3319can be used to control the SHC calculation, including 3320the calculation of alternating current (ac) SHC which will be introduced in 3321the next tutorial. 3322\begin{itemize} 3323\item general controls for SHC 3324 \smallskip {\tt 3325 \begin{quote} 3326 shc\_freq\_scan, 3327 shc\_alpha, 3328 shc\_beta, 3329 shc\_gamma, 3330 3331 kubo\_eigval\_max, 3332 exclude\_bands, 3333 berry\_curv\_unit 3334 \end{quote} } 3335\item kmesh 3336 \smallskip {\tt 3337 \begin{quote} 3338 berry\_task, 3339 berry\_kmesh, 3340 3341 berry\_curv\_adpt\_kmesh, 3342 berry\_curv\_adpt\_kmesh\_thresh 3343 \end{quote} } 3344\item ac SHC 3345 \smallskip {\tt 3346 \begin{quote} 3347 kubo\_freq\_min, 3348 kubo\_freq\_max, 3349 kubo\_freq\_step, 3350 3351 shc\_bandshift, 3352 shc\_bandshift\_firstband, 3353 shc\_bandshift\_energyshift, 3354 3355 scissors\_shift, 3356 num\_valence\_bands 3357 \end{quote} } 3358\item smearing 3359 \smallskip {\tt 3360 \begin{quote} 3361 [kubo\_]adpt\_smr, 3362 [kubo\_]adpt\_smr\_fac, 3363 [kubo\_]adpt\_smr\_max, 3364 3365 [kubo\_]smr\_fixed\_en\_width 3366 \end{quote} } 3367\item Fermi energy 3368 \smallskip {\tt 3369 \begin{quote} 3370 fermi\_energy, 3371 fermi\_energy\_min, 3372 fermi\_energy\_max, 3373 fermi\_energy\_step 3374 \end{quote} } 3375\item kpath 3376 \smallskip {\tt 3377 \begin{quote} 3378 kpath, 3379 kpath\_task, 3380 kpath\_num\_points, 3381 kpath\_bands\_colour 3382 \end{quote} } 3383\item kslice 3384 \smallskip {\tt 3385 \begin{quote} 3386 kslice, 3387 kslice\_task, 3388 kslice\_corner, 3389 kslice\_b1, 3390 kslice\_b2, 3391 kslice\_2dkmesh, 3392 3393 kslice\_fermi\_level, 3394 kslice\_fermi\_lines\_colour 3395 \end{quote} } 3396\end{itemize} 3397Their meanings and usages can be found in Ch.~11.5 of the User Guide. 3398 3399\sectiontitle{30: Gallium Arsenide -- Frequency-dependent spin Hall conductivity} 3400 3401\begin{itemize} 3402 \item{Outline: {\it Calculate the alternating current (ac) spin Hall conductivity 3403 of gallium arsenide considering spin-orbit coupling. 3404 To gain a better understanding of this example, 3405 it is suggested to read Ref.~\cite{qiao-prb2018} for a detailed 3406 description of the theory and Ch.~12.5 of the User Guide.}} 3407 \item{Directory: {\tt examples/example30/}} 3408 \item{Input files} 3409 \begin{itemize} 3410 \item{ {\tt GaAs.scf} {\it The \pwscf\ input file for ground state 3411 calculation}} 3412 \item{ {\tt GaAs.nscf} {\it The \pwscf\ input file to obtain Bloch 3413 states on a uniform grid}} 3414 \item{ {\tt GaAs.pw2wan} {\it The input file for {\tt pw2wannier90}}} 3415 \item{ {\tt GaAs.win} {\it The {\tt wannier90} and {\tt postw90} input file}} 3416 \end{itemize} 3417\end{itemize} 3418 3419\begin{enumerate} 3420 \item Run \pwscf\ to obtain the ground state of gallium arsenide\\ 3421 {\tt pw.x < GaAs.scf > scf.out} 3422 3423 \item Run \pwscf\ to obtain the Bloch states on a uniform $k$-point 3424 grid\\ 3425 {\tt pw.x < GaAs.nscf > nscf.out} 3426 3427 \item Run \wannier\ to generate a list of the required overlaps (written 3428 into the {\tt GaAs.nnkp} file)\\ 3429 {\tt wannier90.x -pp GaAs} 3430 3431 \item Run {\tt pw2wannier90} to compute the overlaps between Bloch 3432 states and the projections for the starting guess (written in the 3433 {\tt GaAs.mmn} and {\tt GaAs.amn} files)\\ 3434 {\tt pw2wannier90.x < GaAs.pw2wan > pw2wan.out} 3435 3436 \item Run \wannier\ to compute the MLWFs\\ 3437 {\tt wannier90.x GaAs} 3438 3439 \item Run \postw\ \\ 3440 {\tt postw90.x GaAs} (serial execution)\\ 3441 {\tt mpirun -np 8 postw90.x GaAs} (example of parallel execution with 3442 8 MPI processes) 3443 3444\end{enumerate} 3445 3446\subsection*{ac spin Hall conductivity} 3447 3448The spin Hall conductivity is also dependent on the 3449frequency $\omega$ in the Eq.~(12.22) of the User Guide. 3450The direct current (dc) SHC calculated in the previous example 3451corresponds to $\sigma_{\alpha\beta}^{\text{spin}\gamma}$ in the limit $\omega\rightarrow 34520$ and it is a real number. 3453At finite frequency $\sigma_{\alpha\beta}^{\text{spin}\gamma}$ acquires an imaginary part. 3454 3455To compute the ac spin Hall conductivity for $\hbar\omega$ 3456up to 8~eV, 3457add the lines 3458% 3459{\tt 3460 \begin{quote} 3461 shc\_freq\_scan = true 3462 3463 kubo\_freq\_min = 0.0 3464 3465 kubo\_freq\_max = 8.0 3466 3467 kubo\_freq\_step = 0.01 3468\end{quote} } 3469% 3470and re-run {\tt postw90}. The file {\tt GaAs-shc-freqscan.dat} contains the calculated ac SHC. Reasonably converged spectra can be 3471obtained with a $250\times 250\times 250$ $k$-point mesh. 3472To plot the ac SHC, issue the following commands {\tt 3473 \begin{quote} 3474 myshell> gnuplot 3475 3476 gnuplot> plot `GaAs-shc-freqscan.dat' u 2:3 w l title `Re', `GaAs-shc-freqscan.dat' u 2:4 w l title `Im' 3477\end{quote} } 3478and then compare the result with Fig.~4 in 3479Ref.~\cite{qiao-prb2018} or the Solution Booklet. 3480 3481\subsection*{Notes} 3482% 3483\begin{itemize} 3484 \item When calculating ac SHC, adaptive smearing can be used 3485 by add the following keywords in the {\tt GaAs.win}, \smallskip {\tt 3486 \begin{quote} 3487 kubo\_adpt\_smr = true 3488 3489 kubo\_adpt\_smr\_fac = [insert here your smearing factor] 3490 3491 kubo\_adpt\_smr\_max = [insert here your maximum smearing] 3492 \end{quote} } 3493 3494 \item Adaptive kmesh refinement is not implemented for ac SHC calculation. 3495 3496 \item The first 10 semi-core states are excluded from the calculation by using the following keywords\smallskip {\tt 3497 \begin{quote} 3498 exclude\_bands = 1-10 3499 \end{quote} } 3500 and in the case of GaAs disentanglement is not adopted so\smallskip {\tt 3501 \begin{quote} 3502 num\_bands = 16 3503 3504 num\_wann = 16 3505 \end{quote} } 3506 \item Since the band gap is often under estimated by LDA/GGA calculations, a scissors shift is applied to recover 3507 the experimental band gap by using the 3508 following keywords\smallskip {\tt 3509 \begin{quote} 3510 shc\_bandshift = true 3511 3512 shc\_bandshift\_firstband = 9 3513 3514 shc\_bandshift\_energyshift = 1.117 3515 \end{quote} } 3516or by\smallskip {\tt 3517 \begin{quote} 3518 num\_valence\_bands = 8 3519 3520 scissors\_shift = 1.117 3521\end{quote} } 3522\end{itemize} 3523 3524 3525\sectiontitle{31: Platinum -- Selected columns of density matrix algorithm for spinor wavefunctions} 3526 3527Note: This example requires a recent version of the {\tt pw2wannier90.x} post-processing code of {\tt Quantum ESPRESSO} (v6.3 or above). 3528 3529\begin{itemize} 3530\item Outline: {\it For bulk crystalline platinum with spin-orbit coupling, generate the $A_{mn}$ matrices via the selected columns of density matrix (SCDM) algorithm and the corresponding spinor-MLWFs. To better understand the input files and the results of these calculations, it is crucial that the Reader has familiarized with the concepts and methods explained in Ref.~\cite{LinLin-ArXiv2017}. More info on the keywords related to the SCDM method may be found in the user\_guide.} 3531 3532This example focuses on the use of the SCDM method for spin-noncollinear systems. For the overview of the use of SCDM method to spinless systems, please refer to example27. 3533\item Directory: {\tt examples/example31/} 3534 3535The input files for this examples are similar to the ones in example 29, except that a coarser k-point grid is used and that the keywords related to {\tt postw90.x} are removed. 3536\item Input Files: 3537\begin{itemize} 3538\item {\tt Pt.scf} {\it The \pwscf\ input file for the ground state calculation} 3539\item {\tt Pt.nscf} {\it The \pwscf\ input file to obtain Bloch states on a uniform grid} 3540\item {\tt Pt.pw2wan} {\it The input file for {\tt pw2wannier90} with keywords related to the SCDM method} 3541\item {\tt Pt.win} {\it The {\tt wannier90} input file} 3542\end{itemize} 3543\end{itemize} 3544 3545 3546We will compute 18 localized WFs. Since the band structure of platinum is metallic, the low-lying bands are entangled with other high-energy bands, and the columns of the density matrix are not exponentially localized by construction. 3547Thus, we use a modified density matrix \cite{LinLin-ArXiv2017}, with the function $f(\varepsilon_{n,\mathbf{k}})$ defined as a complementary error function. Refer to example 27 for the definition of the modified density matrix and the functional form of $f(\varepsilon_{n,\mathbf{k}})$. 3548 3549\begin{enumerate} 3550 \item Run \pwscf\ to obtain the ground state of platinum\\ 3551 {\tt pw.x < Pt.scf > scf.out} 3552 3553 \item Run \pwscf\ to obtain the Bloch states on a uniform $7\times 7\times 7$ $k$-point 3554 grid\\ 3555 {\tt pw.x < Pt.nscf > nscf.out} 3556 3557 \item Inspect the {\tt Pt.win} input file and make sure that the {\tt auto\_projections} flag is set to {\tt .true.}. Also, make sure that no projection block is present. 3558 3559 \item Run \wannier\ to generate a list of the required overlaps (written 3560 into the {\tt Pt.nnkp} file)\\ 3561 {\tt wannier90.x -pp Pt} 3562 3563 \item Inspect the {\tt Pt.nnkp} file and make sure you find the {\tt auto\_projections} block and that no projections have been written in the {\tt projections} block. 3564 3565 \item Inspect the {\tt Pt.pw2wan} input file. You will find four SCDM-related keywords: {\tt scdm\_proj}, {\tt scdm\_entanglement}, {\tt scdm\_mu} and {\tt scdm\_sigma}. 3566 In particular, the keyword {\tt scdm\_proj} will instruct {\tt pw2wannier90.x} to use the SCDM method when generating the $A_{mn}$ matrix. 3567 The remaining three keywords defines the formula and parameters to define the function $f(\varepsilon_{n\mathbf{k}})$ (see Ref.~\cite{LinLin-ArXiv2017} and example 27). 3568 3569 \item Run {\tt pw2wannier90} to compute the overlap between Bloch 3570 states and the projections via the SCDM method (written in the 3571 {\tt Pt.mmn} and {\tt Pt.amn} respectively).\\ 3572 {\tt pw2wannier90.x < Pt.pw2wan > pw2wan.out} 3573 3574 \item Inspect the {\tt pw2wan.out} output file. 3575 Compared to the spinless case, you will find the following two additional lines. 3576 \begin{verbatim} 3577 Number of pivot points with spin up : 9 3578 Number of pivot points with spin down: 9 3579 \end{verbatim} 3580 These lines give information on the pivots obtained by the QR decomposition with column pivoting (QRCP) in the SCDM algorithm. Each pivot determines a point in the real-space grid and a spin state. The basis of the spin state is determined by the basis used in the electronic structure code. In \pwscf, the basis states are spin up and down states along the Cartesian $z$-axis. 3581 3582 \item Run \wannier\ to compute the MLWFs\\ 3583 {\tt wannier90.x Pt}\\ 3584 3585\end{enumerate} 3586 3587\sectiontitle{32: Tungsten --- SCDM parameters from projectability} 3588 3589\begin{itemize} 3590 \item{Outline: {\it Compute the Wannier interpolated band structure of tungsten (W) 3591 using the SCDM method to calculate the initial guess (see Example 27 for more details). The free parameters 3592 in the SCDM method, i.e., $\mu$ and $\sigma$, are obtained by fitting 3593 a complementary error 3594 function to the projectabilities. The number of MLWFs is given by the number 3595 of pseudo-atomic orbitals (PAOs) in the pseudopotential, $21$ in this case. All the steps shown in this example have been automated in the AiiDA\cite{Pizzi_AiiDA} workflow that can be downloaded from the MaterialsCloud website\cite{MaterialsCloudArchiveEntry}}.} 3596 \item{Directory: {\tt examples/example31/}} 3597 \item{Input files} 3598 \begin{itemize} 3599 \item{ {\tt W.scf} {\it The \pwscf\ input file for ground state 3600 calculation}} 3601 \item{ {\tt W.nscf} {\it The \pwscf\ input file to obtain Bloch 3602 states on a uniform grid}} 3603 \item{ {\tt W.pw2wan} {\it The input file for {\tt pw2wannier90}}} 3604 \item{ {\tt W.proj} {\it The input file for {\tt projwfc}}} 3605 \item{ {\tt generate\_weights.sh} {\it The bash script to extract the projectabilities from the output of {\tt projwfc} }} 3606 \item{ {\tt W.win} {\it The {\tt wannier90} input file}} 3607 \end{itemize} 3608\end{itemize} 3609 3610\begin{enumerate} 3611 \item Run \pwscf\ to obtain the ground state of tungsten\\ 3612 {\tt pw.x -in W.scf > scf.out} 3613 3614 \item Run \pwscf\ to obtain the Bloch states on a $10\times10\times10$ uniform $k$-point 3615 grid\\ 3616 {\tt pw.x -in W.nscf > nscf.out} 3617 3618 \item Run \wannier\ to generate a list of the required overlaps (written 3619 into the {\tt W.nnkp} file)\\ 3620 {\tt wannier90.x -pp W} 3621 3622 \item Run {\tt projwfc} to compute the projectabilities of the Bloch states 3623 onto the Bloch sums obtained from the PAOs in the pseudopotential \\ 3624 {\tt projwfc.x -in W.proj > proj.out} 3625 3626 \item Run {\tt generate\_weights} to extract the projectabilitites from {\tt proj.out} 3627 in a format suitable to be read by Xmgrace or gnuplot \\ 3628 {\tt ./generate\_weights.sh} 3629 3630 \item Plot the projectabilities and fit the data with the complementary error function 3631 $$f(\epsilon;\mu,\sigma) = \frac{1}{2}\mathrm{erfc}(-\frac{\mu - \epsilon}{\sigma}).$$ 3632 We are going Xmgrace to plot the projectabilities and perform the fitting. Open Xmgrace \\ 3633 {\tt xmgrace } 3634 3635 To Import the {\tt p\_vs\_e.dat} file, click on {\tt Data} from the top bar and then {\tt Import -> ASCII...}. 3636 At this point a new window {\tt Grace: Read sets} should pop up. Select {\tt p\_vs\_e.dat} in the {\tt Files} section, click {\tt Ok} at the bottom and close the window. You should now be able to see a quite noisy function that is bounded between 1 and 0. You can modify the appearence of the plot by clicking on {\tt Plot} in the top bar and then {\tt Set appearance...}. In the {\tt Main} section of the pop-up window change the symbol type from {\tt None} to {\tt Circle}. Change the line type from straight to none, since the lines added by default by Xmgrace are not meaningful. For the fitting, go to {\tt Data -> Transformations -> Non-linear curve fitting}. In this window, select the source from the {\tt Set} box and in the {\tt Formula} box insert the following \\ 3637 3638 {\tt y = 0.5 * erfc( ( x - A0 ) / A1 )} \\ 3639 3640Select 2 as number of parameters, give 40 as initial condition for {\tt A0} and 7 for {\tt A1}. Click {\tt Apply}. A new window should pop up with the stats of the fitting. In particular you should find a {\tt Correlation coefficient} of 0.96 and a value of $39.9756$ for {\tt A0} and $6.6529$ for {\tt A1}. These are the value of $\mu_{fit}$ and $\sigma_{fit}$ we are going to use for the SCDM method. In particular, $\mu_{SCDM} = \mu_{fit} - 3\sigma_{fit} = 20.0169$ eV and $\sigma_{SCDM} = \sigma_{fit} = 6.6529$ eV. The motivation for this specific choice of $\mu_{fit}$ and $\sigma_{fit}$ may be found in Ref.~\cite{Vitale2019automated}, where the authors also show validation of this approach on a dataset of 200 materials. You should now see the fitting function, as well as the projectabilities, in the graph (see Fig. \ref{fig:W_fit}-(a)).\\ 3641 3642 \item Open {\tt W.pw2wan} and append the following lines 3643 {\tt 3644 3645 scdm\_entanglement = \textquotesingle erfc\textquotesingle 3646 3647 scdm\_mu = 20.0169 3648 3649 scdm\_proj = .true. 3650 3651 scdm\_sigma = 6.6529 3652 3653 /} 3654 3655 \item Run {\tt pw2wannier90} to compute the overlaps between Bloch 3656 states and the projections for the starting guess (written in the 3657 {\tt W.mmn} and {\tt W.amn} files)\\ 3658 {\tt pw2wannier90.x -in W.pw2wan > pw2wan.out} 3659 3660 \item Run \wannier\ to obtain the interpolated bandstructure (see Fig. \ref{fig:W_fit}-(b)).\\ 3661 {\tt wannier90.x W} 3662 3663Please cite Ref.~\cite{Vitale2019automated} in any publication employing the procedure outlined in this example to obtain $\mu$ and $\sigma$. 3664\end{enumerate} 3665 3666\begin{figure} 3667\centering 3668\subfloat[]{\includegraphics[width=0.45\columnwidth]{./W_fit.png}} 3669\subfloat[]{\includegraphics[width=0.45\columnwidth]{./W_bs.pdf}} 3670\caption{ a) Each blue dot represents the projectability as defined in Eq. (22) of Ref. \cite{Vitale2019automated} of the state 3671$\vert n\mathbf{k} \rangle$ as a function of the corresponding energy $\epsilon_{n\mathbf{k}}$ for tungsten. The yellow line shows the fitted complementary error function. The vertical red line represents the value of $\sigma_{fit}$ while the vertical green line represents the optimal value of $\mu_{SCDM}$, i.e. $\mu_{SCDM} = \mu_{fit} - 3\sigma_{fit}$. b) Band structure of tungsten on the $\Gamma$-H-N-$\Gamma$ path from DFT calculations (solid black) and Wannier interpolation using the SCDM method to construct the initial guess (red dots).} 3672\label{fig:W_fit} 3673\end{figure} 3674 3675%\cleardoublepage 3676 3677\bibliographystyle{apsrev4-1} 3678\bibliography{../wannier90} 3679 3680\end{document} 3681
