|
| 1 | +# # Wannierization using Wannier.jl or Wannier90 |
| 2 | +# |
| 3 | +# DFTK features an interface with the programs |
| 4 | +# [Wannier.jl](https://wannierjl.org) and [Wannier90](http://www.wannier.org/), |
| 5 | +# in order to compute maximally-localized Wannier functions (MLWFs) |
| 6 | +# from an initial self consistent field calculation. |
| 7 | +# All processes are handled by calling the routine `wannier_model` (for Wannier.jl) or `run_wannier90` (for Wannier90). |
| 8 | +# |
| 9 | +# !!! warning "No guarantees on Wannier interface" |
| 10 | +# This code is at an early stage and has so far not been fully tested. |
| 11 | +# Bugs are likely and we welcome issues in case you find any! |
| 12 | +# |
| 13 | +# This example shows how to obtain the MLWFs corresponding |
| 14 | +# to the first five bands of graphene. Since the bands 2 to 11 are entangled, |
| 15 | +# 15 bands are first computed to obtain 5 MLWFs by a disantanglement procedure. |
| 16 | + |
| 17 | +using DFTK |
| 18 | +using Plots |
| 19 | +using Unitful |
| 20 | +using UnitfulAtomic |
| 21 | + |
| 22 | +d = 10u"Å" |
| 23 | +a = 2.641u"Å" # Graphene Lattice constant |
| 24 | +lattice = [a -a/2 0; |
| 25 | + 0 √3*a/2 0; |
| 26 | + 0 0 d] |
| 27 | + |
| 28 | +C = ElementPsp(:C; psp=load_psp("hgh/pbe/c-q4")) |
| 29 | +atoms = [C, C] |
| 30 | +positions = [[0.0, 0.0, 0.0], [1//3, 2//3, 0.0]] |
| 31 | +model = model_PBE(lattice, atoms, positions) |
| 32 | +basis = PlaneWaveBasis(model; Ecut=15, kgrid=[5, 5, 1]) |
| 33 | +nbandsalg = AdaptiveBands(basis.model; n_bands_converge=15) |
| 34 | +scfres = self_consistent_field(basis; nbandsalg, tol=1e-5); |
| 35 | + |
| 36 | +# Plot bandstructure of the system |
| 37 | + |
| 38 | +bands = compute_bands(scfres; kline_density=10) |
| 39 | +plot_bandstructure(bands) |
| 40 | + |
| 41 | +# ## Wannierization with Wannier.jl |
| 42 | +# |
| 43 | +# Now we use the `wannier_model` routine to generate a Wannier.jl model |
| 44 | +# that can be used to perform the wannierization procedure. |
| 45 | +# For now, this model generation produces file in the Wannier90 convention, |
| 46 | +# where all files are named with the same prefix and only differ by |
| 47 | +# their extensions. By default all generated input and output files are stored |
| 48 | +# in the subfolder "wannierjl" under the prefix "wannier" (i.e. "wannierjl/wannier.win", |
| 49 | +# "wannierjl/wannier.wout", etc.). A different file prefix can be given with the |
| 50 | +# keyword argument `fileprefix` as shown below. |
| 51 | +# |
| 52 | +# We now produce a simple Wannier model for 5 MLFWs. |
| 53 | +# |
| 54 | +# For a good MLWF, we need to provide initial projections that resemble the expected shape |
| 55 | +# of the Wannier functions. |
| 56 | +# Here we will use: |
| 57 | +# - 3 bond-centered 2s hydrogenic orbitals for the expected σ bonds |
| 58 | +# - 2 atom-centered 2pz hydrogenic orbitals for the expected π bands |
| 59 | + |
| 60 | +using Wannier # Needed to make Wannier.Model available |
| 61 | + |
| 62 | +# From chemical intuition, we know that the bonds with the lowest energy are: |
| 63 | +# - the 3 σ bonds, |
| 64 | +# - the π and π* bonds. |
| 65 | +# We provide relevant initial projections to help Wannierization |
| 66 | +# converge to functions with a similar shape. |
| 67 | +s_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 0, 0, C.Z) |
| 68 | +pz_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 1, 0, C.Z) |
| 69 | +projections = [ |
| 70 | + ## Note: fractional coordinates for the centers! |
| 71 | + ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds |
| 72 | + s_guess((positions[1] + positions[2]) / 2), |
| 73 | + s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), |
| 74 | + s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), |
| 75 | + ## 2 atom-centered 2pz hydrogenic orbitals |
| 76 | + pz_guess(positions[1]), |
| 77 | + pz_guess(positions[2]), |
| 78 | +] |
| 79 | + |
| 80 | +# Wannierize: |
| 81 | +wannier_model = Wannier.Model(scfres; |
| 82 | + fileprefix="wannier/graphene", |
| 83 | + n_bands=scfres.n_bands_converge, |
| 84 | + n_wannier=5, |
| 85 | + projections, |
| 86 | + dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1) # maximum frozen window, for example 1 eV above Fermi level |
| 87 | + |
| 88 | +# Once we have the `wannier_model`, we can use the functions in the Wannier.jl package: |
| 89 | +# |
| 90 | +# Compute MLWF: |
| 91 | +U = disentangle(wannier_model, max_iter=200); |
| 92 | + |
| 93 | +# Inspect localization before and after Wannierization: |
| 94 | +omega(wannier_model) |
| 95 | +omega(wannier_model, U) |
| 96 | + |
| 97 | +# Build a Wannier interpolation model: |
| 98 | +kpath = irrfbz_path(model) |
| 99 | +interp_model = Wannier.InterpModel(wannier_model; kpath=kpath) |
| 100 | + |
| 101 | +# And so on... |
| 102 | +# Refer to the Wannier.jl documentation for further examples. |
| 103 | + |
| 104 | +# (Delete temporary files when done.) |
| 105 | +rm("wannier", recursive=true) |
| 106 | + |
| 107 | +# ### Custom initial guesses |
| 108 | +# |
| 109 | +# We can also provide custom initial guesses for Wannierization, |
| 110 | +# by passing a callable function in the `projections` array. |
| 111 | +# The function receives the basis and a list of points (fractional coordinates in reciprocal space), |
| 112 | +# and returns the Fourier transform of the initial guess function evaluated at each point. |
| 113 | +# |
| 114 | +# For example, we could use Gaussians for the σ and pz guesses with the following code: |
| 115 | +s_guess(center) = DFTK.GaussianWannierProjection(center) |
| 116 | +function pz_guess(center) |
| 117 | + ## Approximate with two Gaussians offset by 0.5 Å from the center of the atom |
| 118 | + offset = model.inv_lattice * [0, 0, austrip(0.5u"Å")] |
| 119 | + center1 = center + offset |
| 120 | + center2 = center - offset |
| 121 | + ## Build the custom projector |
| 122 | + (basis, qs) -> DFTK.GaussianWannierProjection(center1)(basis, qs) - DFTK.GaussianWannierProjection(center2)(basis, qs) |
| 123 | +end |
| 124 | +## Feed to Wannier via the `projections` as before... |
| 125 | + |
| 126 | +# This example assumes that Wannier.jl version 0.3.2 is used, |
| 127 | +# and may need to be updated to accomodate for changes in Wannier.jl. |
| 128 | +# |
| 129 | +# Note: Some parameters supported by Wannier90 may have to be passed to Wannier.jl differently, |
| 130 | +# for example the max number of iterations is passed to `disentangle` in Wannier.jl, |
| 131 | +# but as `num_iter` to `run_wannier90`. |
| 132 | + |
| 133 | +# ## Wannierization with Wannier90 |
| 134 | +# |
| 135 | +# We can use the `run_wannier90` routine to generate all required files and perform the wannierization procedure: |
| 136 | + |
| 137 | +using wannier90_jll # Needed to make run_wannier90 available |
| 138 | +run_wannier90(scfres; |
| 139 | + fileprefix="wannier/graphene", |
| 140 | + n_wannier=5, |
| 141 | + projections, |
| 142 | + num_print_cycles=25, |
| 143 | + num_iter=200, |
| 144 | + ## |
| 145 | + dis_win_max=19.0, |
| 146 | + dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1, # 1 eV above Fermi level |
| 147 | + dis_num_iter=300, |
| 148 | + dis_mix_ratio=1.0, |
| 149 | + ## |
| 150 | + wannier_plot=true, |
| 151 | + wannier_plot_format="cube", |
| 152 | + wannier_plot_supercell=5, |
| 153 | + write_xyz=true, |
| 154 | + translate_home_cell=true, |
| 155 | + ); |
| 156 | + |
| 157 | +# As can be observed standard optional arguments for the disentanglement |
| 158 | +# can be passed directly to `run_wannier90` as keyword arguments. |
| 159 | + |
| 160 | +# (Delete temporary files.) |
| 161 | +rm("wannier", recursive=true) |
0 commit comments