From 5839ce5bf19930d0bf097614c73df320e0ed23f4 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 30 Oct 2023 11:06:44 +0100 Subject: [PATCH 01/29] Add basic Wannier.jl integration --- Project.toml | 1 + docs/make.jl | 2 +- examples/wannier.jl | 148 ++++++++++++++ examples/wannier90.jl | 72 ------- src/DFTK.jl | 5 + src/external/stubs.jl | 20 +- src/external/wannier.jl | 36 ++++ src/external/wannier90.jl | 290 ++------------------------- src/external/wannier_shared.jl | 351 +++++++++++++++++++++++++++++++++ 9 files changed, 572 insertions(+), 353 deletions(-) create mode 100644 examples/wannier.jl delete mode 100644 examples/wannier90.jl create mode 100644 src/external/wannier.jl create mode 100644 src/external/wannier_shared.jl diff --git a/Project.toml b/Project.toml index ebbcf111fb..ce0ecea82c 100644 --- a/Project.toml +++ b/Project.toml @@ -117,6 +117,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" +Wannier = "2b19380a-1f7e-4d7d-b1b8-8aa60b3321c9" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" wannier90_jll = "c5400fa0-8d08-52c2-913f-1e3f656c1ce9" diff --git a/docs/make.jl b/docs/make.jl index 476f816d2d..7742be137d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -47,7 +47,7 @@ PAGES = [ # options we have "examples/atomsbase.jl", "examples/input_output.jl", - "examples/wannier90.jl", + "examples/wannier.jl", ], "Tipps and tricks" => [ # Resolving convergence issues, what solver to use, improving performance or diff --git a/examples/wannier.jl b/examples/wannier.jl new file mode 100644 index 0000000000..2b0e7235a0 --- /dev/null +++ b/examples/wannier.jl @@ -0,0 +1,148 @@ +# # Wannierization using Wannier.jl or Wannier90 +# +# DFTK features an interface with the programs +# [Wannier.jl](https://wannierjl.org) and [Wannier90](http://www.wannier.org/), +# in order to compute maximally-localized Wannier functions (MLWFs) +# from an initial self consistent field calculation. +# All processes are handled by calling the routine `get_wannier_model` (for Wannier.jl) or `run_wannier90` (for Wannier90). +# +# !!! warning "No guarantees on Wannier interface" +# This code is at an early stage and has so far not been fully tested. +# Bugs are likely and we welcome issues in case you find any! +# +# This example shows how to obtain the MLWFs corresponding +# to the first five bands of graphene. Since the bands 2 to 11 are entangled, +# 15 bands are first computed to obtain 5 MLWFs by a disantanglement procedure. + +using DFTK +using Plots +using Unitful +using UnitfulAtomic + +d = 10u"Å" +a = 2.641u"Å" # Graphene Lattice constant +lattice = [a -a/2 0; + 0 √3*a/2 0; + 0 0 d] + +C = ElementPsp(:C, psp=load_psp("hgh/pbe/c-q4")) +atoms = [C, C] +positions = [[0.0, 0.0, 0.0], [1//3, 2//3, 0.0]] +model = model_PBE(lattice, atoms, positions) +basis = PlaneWaveBasis(model; Ecut=15, kgrid=[5, 5, 1]) +nbandsalg = AdaptiveBands(basis.model; n_bands_converge=15) +scfres = self_consistent_field(basis; nbandsalg, tol=1e-5); + +# Plot bandstructure of the system + +plot_bandstructure(scfres; kline_density=10) + +# ## Wannierization with Wannier.jl +# +# Now we use the `get_wannier_model` routine to generate a Wannier.jl model +# that can be used to perform the wannierization procedure. +# For now, this model generation produces file in the Wannier90 convention, +# where all files are named with the same prefix and only differ by +# their extensions. By default all generated input and output files are stored +# in the subfolder "wannierjl" under the prefix "wannier" (i.e. "wannierjl/wannier.win", +# "wannierjl/wannier.wout", etc.). A different file prefix can be given with the +# keyword argument `fileprefix` as shown below. +# +# We now produce a simple Wannier model for 5 MLFWs. +# +# For a good MLWF, we need to provide initial projections that resemble the expected shape +# of the Wannier functions. +# Here we will use: +# - 3 bond-centered gaussians for the expected σ bonds +# - 2 atom-centered opposite gaussians (for each atom: a gaussian with the + sign and a gaussian with the - sign) +# for the expected pz orbitals + +using Wannier # Needed to make get_wannier_model available + +"""Helper function to produce a starting guess for pz orbitals""" +make_double_gaussian(atomindex) = begin + offset = model.inv_lattice * [0, 0, austrip(0.5u"Å")] + center1 = positions[atomindex] + offset + center2 = positions[atomindex] - offset + DFTK.opposite_gaussians_projection(center1, center2) +end + +wannier_model = get_wannier_model(scfres; + fileprefix="wannier/graphene", + n_bands=scfres.n_bands_converge, + n_wannier=5, + projections=[ + ## Note: fractional coordinates for the centers! + ## 3 bond-centered gaussians to imitate σ orbitals + DFTK.GaussianWannierProjection((positions[1] + positions[2]) / 2), + DFTK.GaussianWannierProjection((positions[1] + positions[2] + [0, -1, 0]) / 2), + DFTK.GaussianWannierProjection((positions[1] + positions[2] + [-1, -1, 0]) / 2), + ## 2 double gaussians to imitate pz orbitals + make_double_gaussian(1), + make_double_gaussian(2), + ], + dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1) # maximum frozen window, for example 1 eV above Fermi level + +# Once we have the `wannier_model`, we can use the functions in the Wannier.jl package: +# +# Compute MLWF: +U = disentangle(wannier_model, max_iter=200); + +# Inspect localization before and after Wannierization: +omega(wannier_model) +omega(wannier_model, U) + +# Build a Wannier interpolation model: (needs a kpath) +interp_model = Wannier.InterpModel(wannier_model; kpath=kpath) + +# And so on... +# Refer to the Wannier.jl documentation for further examples. + +# (Delete temporary files when done.) +rm("wannier", recursive=true) + +# This example assumes that Wannier.jl version 0.3.1 is used, +# and may need to be updated to accomodate for changes in Wannier.jl. +# +# Note: Some parameters supported by Wannier90 may have to be passed to Wannier.jl differently, +# for example the max number of iterations is passed to `disentangle` in Wannier.jl, +# but as `num_iter` to `run_wannier90`. + +# ## Wannierization with Wannier90 +# +# We can use the `run_wannier90` routine to generate all required files and perform the wannierization procedure: + +using wannier90_jll # Needed to make run_wannier90 available +run_wannier90(scfres; + fileprefix="wannier/graphene", + n_wannier=5, + projections=[ + ## Note: fractional coordinates for the centers! + ## 3 bond-centered gaussians to imitate σ orbitals + DFTK.GaussianWannierProjection((positions[1] + positions[2]) / 2), + DFTK.GaussianWannierProjection((positions[1] + positions[2] + [0, -1, 0]) / 2), + DFTK.GaussianWannierProjection((positions[1] + positions[2] + [-1, -1, 0]) / 2), + ## 2 double gaussians to imitate pz orbitals + make_double_gaussian(1), + make_double_gaussian(2), + ], + num_print_cycles=25, + num_iter=200, + ## + dis_win_max=19.0, + dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1, # 1 eV above Fermi level + dis_num_iter=300, + dis_mix_ratio=1.0, + ## + wannier_plot=true, + wannier_plot_format="cube", + wannier_plot_supercell=5, + write_xyz=true, + translate_home_cell=true, + ); + +# As can be observed standard optional arguments for the disentanglement +# can be passed directly to `run_wannier90` as keyword arguments. + +# (Delete temporary files.) +rm("wannier", recursive=true) diff --git a/examples/wannier90.jl b/examples/wannier90.jl deleted file mode 100644 index aa4e0e4d92..0000000000 --- a/examples/wannier90.jl +++ /dev/null @@ -1,72 +0,0 @@ -# # Wannierization using Wannier90 -# -# DFTK features an interface with the program -# [Wannier90](http://www.wannier.org/), -# in order to compute maximally-localized Wannier functions (MLWFs) -# from an initial self consistent field calculation. -# All processes are handled by calling the routine `run_wannier90`. -# -# !!! warning "No guarantees on Wannier90 interface" -# This code is at an early stage and has so far not been fully tested. -# Bugs are likely and we welcome issues in case you find any! -# -# This example shows how to obtain the MLWFs corresponding -# to the first five bands of graphene. Since the bands 2 to 11 are entangled, -# 15 bands are first computed to obtain 5 MLWFs by a disantanglement procedure. - -using DFTK -using Unitful -using UnitfulAtomic - -d = 10u"Å" -a = 2.641u"Å" # Graphene Lattice constant -lattice = [a -a/2 0; - 0 √3*a/2 0; - 0 0 d] - -C = ElementPsp(:C, psp=load_psp("hgh/pbe/c-q4")) -atoms = [C, C] -positions = [[0.0, 0.0, 0.0], [1//3, 2//3, 0.0]] -model = model_PBE(lattice, atoms, positions) -basis = PlaneWaveBasis(model; Ecut=15, kgrid=[5, 5, 1]) -nbandsalg = AdaptiveBands(basis.model; n_bands_converge=15) -scfres = self_consistent_field(basis; nbandsalg, tol=1e-5); - -# Plot bandstructure of the system - -plot_bandstructure(scfres; kline_density=10) - -# Now we use the `run_wannier90` routine to generate all files needed by -# wannier90 and to perform the wannierization procedure. -# In Wannier90's convention, all files are named with the same prefix and only differ by -# their extensions. By default all generated input and output files are stored -# in the subfolder "wannier90" under the prefix "wannier" (i.e. "wannier90/wannier.win", -# "wannier90/wannier.wout", etc.). A different file prefix can be given with the -# keyword argument `fileprefix` as shown below. -# -# We now solve for 5 MLWF using wannier90: - -using wannier90_jll # Needed to make run_wannier90 available -run_wannier90(scfres; - fileprefix="wannier/graphene", - n_wannier=5, - num_print_cycles=25, - num_iter=200, - ## - dis_win_max=19.0, - dis_froz_max=0.1, - dis_num_iter=300, - dis_mix_ratio=1.0, - ## - wannier_plot=true, - wannier_plot_format="cube", - wannier_plot_supercell=5, - write_xyz=true, - translate_home_cell=true, - ); - -# As can be observed standard optional arguments for the disentanglement -# can be passed directly to `run_wannier90` as keyword arguments. - -# (Delete temporary files.) -rm("wannier", recursive=true) diff --git a/src/DFTK.jl b/src/DFTK.jl index e4603e5b1c..67dd2927aa 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -183,10 +183,12 @@ include("pseudo/attach_psp.jl") export DFTKPotential export atomic_system, periodic_system # Reexport from AtomsBase +export get_wannier_model export run_wannier90 include("external/atomsbase.jl") include("external/interatomicpotentials.jl") include("external/stubs.jl") # Function stubs for conditionally defined methods +include("external/wannier_shared.jl") export compute_bands export plot_bandstructure @@ -238,6 +240,9 @@ function __init__() @require wannier90_jll="c5400fa0-8d08-52c2-913f-1e3f656c1ce9" begin include("external/wannier90.jl") end + @require Wannier="2b19380a-1f7e-4d7d-b1b8-8aa60b3321c9" begin + include("external/wannier.jl") + end @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin include("workarounds/cuda_arrays.jl") end diff --git a/src/external/stubs.jl b/src/external/stubs.jl index 321587c63d..00c099381d 100644 --- a/src/external/stubs.jl +++ b/src/external/stubs.jl @@ -1,13 +1,27 @@ # Stubs for conditionally defined functions """ -Wannerize the obtained bands using wannier90. By default all converged +Build a Wannier.jl model for the obtained bands. By default all converged bands from the `scfres` are employed (change with `n_bands` kwargs) -and `n_wannier = n_bands` wannier functions are computed using -random Gaussians as guesses. All keyword arguments supported by +and `n_wannier = n_bands` wannier functions are computed. +Random Gaussians are used as guesses by default, can be changed using the `projections` kwarg. +All keyword arguments supported by Wannier90 for the disentanglement may be added as keyword arguments. The function returns the `fileprefix`. +!!! warning "Experimental feature" + Currently this is an experimental feature, which has not yet been tested + to full depth. The interface is considered unstable and may change + incompatibly in the future. Use at your own risk and please report bugs + in case you encounter any. +""" +function get_wannier_model end + +""" +Wannerize the obtained bands using wannier90. +Same arguments as `get_wannier_model`. +The function returns the `fileprefix`. + !!! warning "Experimental feature" Currently this is an experimental feature, which has not yet been tested to full depth. The interface is considered unstable and may change diff --git a/src/external/wannier.jl b/src/external/wannier.jl new file mode 100644 index 0000000000..e4b0a33a46 --- /dev/null +++ b/src/external/wannier.jl @@ -0,0 +1,36 @@ +""" +Use Wannier.jl to read a .win file and produce the nnkp data (i.e. the b-vectors used in the Mmn matrix). +""" +function get_nnkpt_from_wannier(fileprefix) + model = Wannier.read_w90(fileprefix; amn=false, mmn=false, eig=false) + bvectors = model.bvectors + + nnkpts = reduce(vcat, + map(1:model.n_kpts) do ik + zp = zip(bvectors.kpb_k[:, ik], eachcol(bvectors.kpb_b[:, :, ik])) + map(zp) do (ikpb, G_shift) + (ik, ikpb, G_shift=copy(G_shift)) + end + end) + + (nntot=model.n_bvecs, nnkpts) +end + +""" +Build a Wannier.jl model that can be used for Wannierization. +""" +@timing function get_wannier_model(scfres; + n_bands=scfres.n_bands_converge, + n_wannier=n_bands, + projections::AbstractVector{<:WannierProjection}=default_wannier_centres(n_wannier), + fileprefix=joinpath("wannierjl", "wannier"), + wannier_plot=false, + kwargs...) + # Write the files + write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do + get_nnkpt_from_wannier(fileprefix) + end + + # Read Wannier.jl model + Wannier.read_w90(fileprefix) +end \ No newline at end of file diff --git a/src/external/wannier90.jl b/src/external/wannier90.jl index 40defa3020..0fdc1a79cb 100644 --- a/src/external/wannier90.jl +++ b/src/external/wannier90.jl @@ -1,288 +1,24 @@ -using Dates """ -Write a win file at the indicated prefix. -Parameters to Wannier90 can be added as kwargs: e.g. `num_iter=500`. +Run the Wannierization procedure with Wannier90. """ -function write_w90_win(fileprefix::String, basis::PlaneWaveBasis; - bands_plot=false, wannier_plot=false, kwargs...) - @assert :num_bands in keys(kwargs) # Required parameter - @assert :num_wann in keys(kwargs) # Required parameter - isnothing(basis.kgrid) && error("The basis must be constructed from a MP grid.") - - open(fileprefix *".win", "w") do fp - println(fp, "! Generated by DFTK.jl at $(now())\n") - - # Drop parameters - for (key, value) in pairs(kwargs) - @printf fp "%-20s = %-30s\n" key value - end - if wannier_plot - println(fp, "wvfn_formatted = True") - println(fp, "wannier_plot = True") - end - - # Delimiter for system - println(fp, "\n", "!"^20 * " System \n\n") - - # Lattice vectors are in rows in Wannier90 - println(fp, "begin unit_cell_cart\nbohr") - for vec in eachcol(basis.model.lattice) - @printf fp "%10.6f %10.6f %10.6f \n" vec... - end - println(fp, "end unit_cell_cart \n") - - println(fp, "begin atoms_frac") - for (element, position) in zip(basis.model.atoms, basis.model.positions) - @printf fp "%-2s %10.6f %10.6f %10.6f \n" element.symbol position... - end - println(fp, "end atoms_frac\n") - - # Delimiter for k-points block - println(fp, "!"^20 * " k_points\n") - - if bands_plot - kpath = irrfbz_path(basis.model) - length(kpath.paths) > 1 || @warn( - "Only first kpath branch considered in write_w90_win") - path = kpath.paths[1] - - println(fp, "begin kpoint_path") - for i in 1:length(path)-1 - A, B = path[i:i+1] # write segment A -> B - @printf(fp, "%s %10.6f %10.6f %10.6f ", A, round.(kpath.points[A], digits=5)...) - @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpath.points[B], digits=5)...) - end - println(fp, "end kpoint_path") - println(fp, "bands_plot = true\n") - end - - println(fp, "mp_grid : $(basis.kgrid[1]) $(basis.kgrid[2]) $(basis.kgrid[3])\n") - println(fp, "begin kpoints") - for kpt in basis.kpoints - @printf fp "%10.6f %10.6f %10.6f\n" kpt.coordinate... - end - println(fp, "end kpoints") - end - nothing -end - - -""" -Read the .nnkp file provided by the preprocessing routine of Wannier90 -(i.e. "wannier90.x -pp prefix") -Returns: -1) the array 'nnkpts' of k points, their respective nearest neighbors and associated - shifing vectors (non zero if the neighbor is located in another cell). -2) the number 'nntot' of neighbors per k point. - -TODO: add the possibility to exclude bands -""" -function read_w90_nnkp(fileprefix::String) - fn = fileprefix * ".nnkp" - !isfile(fn) && error("Expected file $fileprefix.nnkp not found.") - - # Extract nnkp block - lines = readlines(fn) - ibegin_nnkp = findfirst(lines .== "begin nnkpts") - iend_nnkp = findfirst(lines .== "end nnkpts") - @assert !isnothing(ibegin_nnkp) && !isnothing(iend_nnkp) - - nntot = parse(Int, lines[ibegin_nnkp + 1]) - nnkpts = map(ibegin_nnkp+2:iend_nnkp-1) do iline - splitted = parse.(Int, split(lines[iline], ' ', keepempty=false)) - @assert length(splitted) == 5 - - # Meaning of the entries: - # 1st: Index of k-point - # 2nd: Index of periodic image of k+b k-point - # 3rd: Shift vector to get k-point of ikpb to the actual k+b point required - (ik=splitted[1], ikpb=splitted[2], G_shift=splitted[3:5]) - end - (; nntot, nnkpts) -end - - -""" -Write the eigenvalues in a format readable by Wannier90. -""" -function write_w90_eig(fileprefix::String, eigenvalues; n_bands) - open("$fileprefix.eig", "w") do fp - for (k, εk) in enumerate(eigenvalues) - for (n, εnk) in enumerate(εk[1:n_bands]) - @printf fp "%3i %3i %25.18f \n" n k auconvert(Unitful.eV, εnk).val - end - end - end - nothing -end - - -@timing function write_w90_unk(fileprefix::String, basis, ψ; n_bands, spin=1) - fft_size = basis.fft_size - - for ik in krange_spin(basis, spin) - open(dirname(fileprefix) * (@sprintf "/UNK%05i.%i" ik spin), "w") do fp - println(fp, "$(fft_size[1]) $(fft_size[2]) $(fft_size[3]) $ik $n_bands") - for iband in 1:n_bands - ψnk_real = ifft(basis, basis.kpoints[ik], @view ψ[ik][:, iband]) - for iz in 1:fft_size[3], iy in 1:fft_size[2], ix in 1:fft_size[1] - println(fp, real(ψnk_real[ix, iy, iz]), " ", imag(ψnk_real[ix, iy, iz])) - end - end - end - end - nothing -end - - -@doc raw""" -Computes the matrix ``[M^{k,b}]_{m,n} = \langle u_{m,k} | u_{n,k+b} \rangle`` -for given `k`, `kpb` = ``k+b``. - -`G_shift` is the "shifting" vector, correction due to the periodicity conditions -imposed on ``k \to ψ_k``. -It is non zero if `kpb` is taken in another unit cell of the reciprocal lattice. -We use here that: -``u_{n(k + G_{\rm shift})}(r) = e^{-i*\langle G_{\rm shift},r \rangle} u_{nk}``. -""" -@views function overlap_Mmn_k_kpb(basis::PlaneWaveBasis, ψ, ik, ikpb, G_shift, n_bands) - # Search for common Fourier modes and their resp. indices in Bloch states k and kpb - # TODO Check if this can be improved using the G vector mapping in the kpoints - k = basis.kpoints[ik] - kpb = basis.kpoints[ikpb] - equivalent_G_vectors = [(iGk, index_G_vectors(basis, kpb, Gk + G_shift)) - for (iGk, Gk) in enumerate(G_vectors(basis, k))] - iGk = [eqG[1] for eqG in equivalent_G_vectors if !isnothing(eqG[2])] - iGkpb = [eqG[2] for eqG in equivalent_G_vectors if !isnothing(eqG[2])] - - # Compute overlaps - # TODO This should be improved ... - Mkb = zeros(ComplexF64, (n_bands, n_bands)) - for n in 1:n_bands - for m in 1:n_bands - # Select the coefficient in right order - Mkb[m, n] = dot(ψ[ik][iGk, m], ψ[ikpb][iGkpb, n]) - end - end - iszero(Mkb) && return Matrix(I, n_bands, n_bands) - Mkb -end - -@timing function write_w90_mmn(fileprefix::String, basis::PlaneWaveBasis, ψ, nnkp; n_bands) - open(fileprefix * ".mmn", "w") do fp - println(fp, "Generated by DFTK at $(now())") - println(fp, "$n_bands $(length(ψ)) $(nnkp.nntot)") - for (ik, ikpb, G_shift) in nnkp.nnkpts - @printf fp "%i %i %i %i %i \n" ik ikpb G_shift... - for ovlp in overlap_Mmn_k_kpb(basis, ψ, ik, ikpb, G_shift, n_bands) - @printf fp "%22.18f %22.18f \n" real(ovlp) imag(ovlp) - end - end - end - nothing -end - - -@doc raw""" -Compute the matrix ``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle`` - -``g^{per}_n`` are periodized gaussians whose respective centers are given as an - (num_bands,1) array [ [center 1], ... ]. - -Centers are to be given in lattice coordinates and G_vectors in reduced coordinates. -The dot product is computed in the Fourier space. - -Given an orbital ``g_n``, the periodized orbital is defined by : - ``g^{per}_n = \sum\limits_{R \in {\rm lattice}} g_n( \cdot - R)``. -The Fourier coefficient of ``g^{per}_n`` at any G -is given by the value of the Fourier transform of ``g_n`` in G. -""" -function compute_Ak_gaussian_guess(basis::PlaneWaveBasis, ψk, kpt, centers, n_bands) - n_wannier = length(centers) - # TODO This function should be improved in performance - - # associate a center with the fourier transform of the corresponding gaussian - fourier_gn(center, qs) = [exp( 2π*(-im*dot(q, center) - dot(q, q) / 4) ) for q in qs] - qs = vec(map(G -> G .+ kpt.coordinate, G_vectors(basis))) # all q = k+G in reduced coordinates - Ak = zeros(eltype(ψk), (n_bands, n_wannier)) - - # Compute Ak - for n in 1:n_wannier - # Functions are l^2 normalized in Fourier, in DFTK conventions. - norm_gn_per = norm(fourier_gn(centers[n], qs), 2) - # Fourier coeffs of gn_per for k+G in common with ψk - coeffs_gn_per = fourier_gn(centers[n], qs[kpt.mapping]) ./ norm_gn_per - # Compute overlap - for m in 1:n_bands - # TODO Check the ordering of m and n here! - Ak[m, n] = dot(ψk[:, m], coeffs_gn_per) - end - end - Ak -end - - -@timing function write_w90_amn(fileprefix::String, basis::PlaneWaveBasis, ψ; n_bands, centers) - open(fileprefix * ".amn", "w") do fp - println(fp, "Generated by DFTK at $(now())") - println(fp, "$n_bands $(length(basis.kpoints)) $(length(centers))") - - for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints)) - Ak = compute_Ak_gaussian_guess(basis, ψk, kpt, centers, n_bands) - for n in 1:size(Ak, 2) - for m in 1:size(Ak, 1) - @printf(fp, "%3i %3i %3i %22.18f %22.18f \n", - m, n, ik, real(Ak[m, n]), imag(Ak[m, n])) - end - end - end - end - nothing -end - -""" -Default random Gaussian guess for maximally-localised wannier functions -generated in reduced coordinates. -""" -default_wannier_centres(n_wannier) = [rand(1, 3) for _ in 1:n_wannier] - @timing function run_wannier90(scfres; - n_bands=scfres.n_bands_converge, - n_wannier=n_bands, - centers=default_wannier_centres(n_wannier), - fileprefix=joinpath("wannier90", "wannier"), - wannier_plot=false, kwargs...) - # TODO None of the routines consider spin at the moment - @assert scfres.basis.model.spin_polarization in (:none, :spinless) - @assert length(centers) == n_wannier - - # Undo symmetry operations to get full k-point list - scfres_unfold = unfold_bz(scfres) - basis = scfres_unfold.basis - ψ = scfres_unfold.ψ - - # Make wannier directory ... - dir, prefix = dirname(fileprefix), basename(fileprefix) - mkpath(dir) - - # Write input file and launch Wannier90 preprocessing - write_w90_win(fileprefix, basis; - num_wann=length(centers), num_bands=n_bands, wannier_plot, kwargs...) - wannier90_jll.wannier90(exe -> run(Cmd(`$exe -pp $prefix`; dir))) - nnkp = read_w90_nnkp(fileprefix) + n_bands=scfres.n_bands_converge, + n_wannier=n_bands, + projections::AbstractVector{<:WannierProjection}=default_wannier_centres(n_wannier), + fileprefix=joinpath("wannier90", "wannier"), + wannier_plot=false, + kwargs...) - # Files for main wannierization run - write_w90_eig(fileprefix, scfres_unfold.eigenvalues; n_bands) - write_w90_amn(fileprefix, basis, ψ; n_bands, centers) - write_w90_mmn(fileprefix, basis, ψ, nnkp; n_bands) + prefix = basename(fileprefix) - # Writing the unk files is expensive (requires FFTs), so only do if needed. - if wannier_plot - write_w90_unk(fileprefix, basis, ψ; n_bands, spin=1) + # Prepare files + write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do + wannier90_jll.wannier90(exe -> run(Cmd(`$exe -pp $prefix`; dir))) + read_w90_nnkp(fileprefix) end # Run Wannierisation procedure @timing "Wannierization" wannier90_jll.wannier90(exe -> run(Cmd(`$exe $prefix`; dir))) fileprefix -end +end \ No newline at end of file diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl new file mode 100644 index 0000000000..7ab8741caf --- /dev/null +++ b/src/external/wannier_shared.jl @@ -0,0 +1,351 @@ +# Shared functionality to operate with Wannier programs (Wannier.jl and Wannier90). +# The communication happens with the Wannier90 file formats which Wannier.jl also understands. + +using Dates + +@doc raw""" +Base type for Wannier projections. + +Wannierization searchs for a unitary matrix ``U_{m_n}``. +As a starting point for the search, we can provide an initial guess function ``g`` +for the shape of the Wannier functions, based on what we expect from knowledge of the problem or physical intuition. +This starting matrix is called ``[A_k]_{m,n}``, and is computed as follows: +``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle``. +The matrix will be orthonormalized by the chosen Wannier program, we don't need to do so ourselves. +""" +abstract type WannierProjection end + +# Dispatch function for WannierProjection implementations, evaluating g in Fourier space. +function get_fourier_projection_coefficients end + +""" +A Gaussian-shaped initial guess. Can be used as an approximation of an s- or σ-like orbital. +""" +struct GaussianWannierProjection <: WannierProjection + center::AbstractVector{Float64} +end + +function get_fourier_projection_coefficients(proj::GaussianWannierProjection, qs) + # associate a center with the fourier transform of the corresponding gaussian + # TODO: what is the normalization here? + [exp( 2π*(-im*dot(q, proj.center) - dot(q, q) / 4) ) for q in qs] +end + +function random_gaussian_projection() + GaussianWannierProjection(rand(3)) +end + +""" +A sum of two Gaussians. Can be used as an approximation of a p-like orbital. +""" +struct DoubleGaussianWannierProjection <: WannierProjection + factor1::Number + proj1::GaussianWannierProjection + factor2::Number + proj2::GaussianWannierProjection +end + +function get_fourier_projection_coefficients(proj::DoubleGaussianWannierProjection, qs) + get_fourier_projection_coefficients(proj.proj1, qs) .+ (proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, qs)) +end + +function DoubleGaussianWannierProjection(factor1::Number, center1::AbstractVector{Float64}, factor2::Number, center2::AbstractVector{Float64}) + DoubleGaussianWannierProjection( + factor1, + GaussianWannierProjection(center1), + factor2, + GaussianWannierProjection(center2), + ) +end + +function opposite_gaussians_projection(center1::AbstractVector{Float64}, center2::AbstractVector{Float64}) + DoubleGaussianWannierProjection(1, center1, -1, center2) +end + +""" +Write a win file at the indicated prefix. +Parameters to Wannier90 can be added as kwargs: e.g. `num_iter=500`. +""" +function write_w90_win(fileprefix::String, basis::PlaneWaveBasis; + bands_plot=false, wannier_plot=false, kwargs...) + @assert :num_bands in keys(kwargs) # Required parameter + @assert :num_wann in keys(kwargs) # Required parameter + isnothing(basis.kgrid) && error("The basis must be constructed from a MP grid.") + + open(fileprefix *".win", "w") do fp + println(fp, "! Generated by DFTK.jl at $(now())\n") + + # Drop parameters + for (key, value) in pairs(kwargs) + @printf fp "%-20s = %-30s\n" key value + end + if wannier_plot + println(fp, "wvfn_formatted = True") + println(fp, "wannier_plot = True") + end + + # Delimiter for system + println(fp, "\n", "!"^20 * " System \n\n") + + # Lattice vectors are in rows in Wannier90 + println(fp, "begin unit_cell_cart\nbohr") + for vec in eachcol(basis.model.lattice) + @printf fp "%10.6f %10.6f %10.6f \n" vec... + end + println(fp, "end unit_cell_cart \n") + + println(fp, "begin atoms_frac") + for (element, position) in zip(basis.model.atoms, basis.model.positions) + @printf fp "%-2s %10.6f %10.6f %10.6f \n" element.symbol position... + end + println(fp, "end atoms_frac\n") + + # Delimiter for k-points block + println(fp, "!"^20 * " k_points\n") + + if bands_plot + kpath = irrfbz_path(basis.model) + length(kpath.paths) > 1 || @warn( + "Only first kpath branch considered in write_w90_win") + path = kpath.paths[1] + + println(fp, "begin kpoint_path") + for i in 1:length(path)-1 + A, B = path[i:i+1] # write segment A -> B + @printf(fp, "%s %10.6f %10.6f %10.6f ", A, round.(kpath.points[A], digits=5)...) + @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpath.points[B], digits=5)...) + end + println(fp, "end kpoint_path") + println(fp, "bands_plot = true\n") + end + + println(fp, "mp_grid : $(basis.kgrid[1]) $(basis.kgrid[2]) $(basis.kgrid[3])\n") + println(fp, "begin kpoints") + for kpt in basis.kpoints + @printf fp "%10.10f %10.10f %10.10f\n" kpt.coordinate... + end + println(fp, "end kpoints") + end + nothing +end + + +""" +Read the .nnkp file provided by the preprocessing routine of Wannier90 +(i.e. "wannier90.x -pp prefix") +Returns: +1) the array 'nnkpts' of k points, their respective nearest neighbors and associated + shifing vectors (non zero if the neighbor is located in another cell). +2) the number 'nntot' of neighbors per k point. + +TODO: add the possibility to exclude bands +""" +function read_w90_nnkp(fileprefix::String) + fn = fileprefix * ".nnkp" + !isfile(fn) && error("Expected file $fileprefix.nnkp not found.") + + # Extract nnkp block + lines = readlines(fn) + ibegin_nnkp = findfirst(lines .== "begin nnkpts") + iend_nnkp = findfirst(lines .== "end nnkpts") + @assert !isnothing(ibegin_nnkp) && !isnothing(iend_nnkp) + + nntot = parse(Int, lines[ibegin_nnkp + 1]) + nnkpts = map(ibegin_nnkp+2:iend_nnkp-1) do iline + splitted = parse.(Int, split(lines[iline], ' ', keepempty=false)) + @assert length(splitted) == 5 + + # Meaning of the entries: + # 1st: Index of k-point + # 2nd: Index of periodic image of k+b k-point + # 3rd: Shift vector to get k-point of ikpb to the actual k+b point required + (ik=splitted[1], ikpb=splitted[2], G_shift=splitted[3:5]) + end + (; nntot, nnkpts) +end + + +""" +Write the eigenvalues in a format readable by Wannier90. +""" +function write_w90_eig(fileprefix::String, eigenvalues; n_bands) + open("$fileprefix.eig", "w") do fp + for (k, εk) in enumerate(eigenvalues) + for (n, εnk) in enumerate(εk[1:n_bands]) + @printf fp "%3i %3i %25.18f \n" n k auconvert(Unitful.eV, εnk).val + end + end + end + nothing +end + + +@timing function write_w90_unk(fileprefix::String, basis, ψ; n_bands, spin=1) + fft_size = basis.fft_size + + for ik in krange_spin(basis, spin) + open(dirname(fileprefix) * (@sprintf "/UNK%05i.%i" ik spin), "w") do fp + println(fp, "$(fft_size[1]) $(fft_size[2]) $(fft_size[3]) $ik $n_bands") + for iband in 1:n_bands + ψnk_real = ifft(basis, basis.kpoints[ik], @view ψ[ik][:, iband]) + for iz in 1:fft_size[3], iy in 1:fft_size[2], ix in 1:fft_size[1] + @printf fp "%25.18f %25.18f\n" real(ψnk_real[ix, iy, iz]) imag(ψnk_real[ix, iy, iz]) + end + end + end + end + nothing +end + + +@doc raw""" +Computes the matrix ``[M^{k,b}]_{m,n} = \langle u_{m,k} | u_{n,k+b} \rangle`` +for given `k`, `kpb` = ``k+b``. + +`G_shift` is the "shifting" vector, correction due to the periodicity conditions +imposed on ``k \to ψ_k``. +It is non zero if `kpb` is taken in another unit cell of the reciprocal lattice. +We use here that: +``u_{n(k + G_{\rm shift})}(r) = e^{-i*\langle G_{\rm shift},r \rangle} u_{nk}``. +""" +@views function overlap_Mmn_k_kpb(basis::PlaneWaveBasis, ψ, ik, ikpb, G_shift, n_bands) + # Search for common Fourier modes and their resp. indices in Bloch states k and kpb + # TODO Check if this can be improved using the G vector mapping in the kpoints + k = basis.kpoints[ik] + kpb = basis.kpoints[ikpb] + equivalent_G_vectors = [(iGk, index_G_vectors(basis, kpb, Gk + G_shift)) + for (iGk, Gk) in enumerate(G_vectors(basis, k))] + iGk = [eqG[1] for eqG in equivalent_G_vectors if !isnothing(eqG[2])] + iGkpb = [eqG[2] for eqG in equivalent_G_vectors if !isnothing(eqG[2])] + + # Compute overlaps + # TODO This should be improved ... + Mkb = zeros(ComplexF64, (n_bands, n_bands)) + for n in 1:n_bands + for m in 1:n_bands + # Select the coefficient in right order + Mkb[m, n] = dot(ψ[ik][iGk, m], ψ[ikpb][iGkpb, n]) + end + end + iszero(Mkb) && return Matrix(I, n_bands, n_bands) + Mkb +end + +@timing function write_w90_mmn(fileprefix::String, basis::PlaneWaveBasis, ψ, nnkp; n_bands) + open(fileprefix * ".mmn", "w") do fp + println(fp, "Generated by DFTK at $(now())") + println(fp, "$n_bands $(length(ψ)) $(nnkp.nntot)") + for (ik, ikpb, G_shift) in nnkp.nnkpts + @printf fp "%i %i %i %i %i \n" ik ikpb G_shift... + for ovlp in overlap_Mmn_k_kpb(basis, ψ, ik, ikpb, G_shift, n_bands) + @printf fp "%22.18f %22.18f \n" real(ovlp) imag(ovlp) + end + end + end + nothing +end + +@doc raw""" +Compute the matrix ``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle`` + +``g^{per}_n`` are periodized gaussians whose respective centers are given as an + (num_bands,1) array [ [center 1], ... ]. + +Centers are to be given in lattice coordinates and G_vectors in reduced coordinates. +The dot product is computed in the Fourier space. + +Given an orbital ``g_n``, the periodized orbital is defined by : + ``g^{per}_n = \sum\limits_{R \in {\rm lattice}} g_n( \cdot - R)``. +The Fourier coefficient of ``g^{per}_n`` at any G +is given by the value of the Fourier transform of ``g_n`` in G. +""" +function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::AbstractVector{<:WannierProjection}, n_bands) + n_wannier = length(projections) + # TODO This function should be improved in performance + + qs = vec(map(G -> G .+ kpt.coordinate, G_vectors(basis))) # all q = k+G in reduced coordinates + Ak = zeros(eltype(ψk), (n_bands, n_wannier)) + + # Compute Ak + for n in 1:n_wannier + proj = projections[n] + # Functions are l^2 normalized in Fourier, in DFTK conventions. + norm_gn_per = norm(get_fourier_projection_coefficients(proj, qs), 2) + # Fourier coeffs of gn_per for k+G in common with ψk + coeffs_gn_per = get_fourier_projection_coefficients(proj, qs[kpt.mapping]) ./ norm_gn_per + # Compute overlap + for m in 1:n_bands + # TODO Check the ordering of m and n here! + Ak[m, n] = dot(ψk[:, m], coeffs_gn_per) + end + end + Ak +end + + +@timing function write_w90_amn( + fileprefix::String, + basis::PlaneWaveBasis, + projections::AbstractVector{<:WannierProjection}, + ψ; + n_bands) + open(fileprefix * ".amn", "w") do fp + println(fp, "Generated by DFTK at $(now())") + println(fp, "$n_bands $(length(basis.kpoints)) $(length(projections))") + + for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints)) + Ak = compute_amn_kpoint(basis, ψk, kpt, projections, n_bands) + for n in 1:size(Ak, 2) + for m in 1:size(Ak, 1) + @printf(fp, "%3i %3i %3i %22.18f %22.18f \n", + m, n, ik, real(Ak[m, n]), imag(Ak[m, n])) + end + end + end + end + nothing +end + +""" +Default random Gaussian guess for maximally-localised wannier functions +generated in reduced coordinates. +""" +default_wannier_centres(n_wannier) = [random_gaussian_projection() for _ in 1:n_wannier] + +""" +Shared file writing code for Wannier.jl and Wannier90. +""" +@timing function write_wannier90_files(preprocess_call::Function, scfres; + n_bands, + n_wannier, + projections::AbstractVector{<:WannierProjection}, + fileprefix, + wannier_plot, + kwargs...) + # TODO None of the routines consider spin at the moment + @assert scfres.basis.model.spin_polarization in (:none, :spinless) + @assert length(projections) == n_wannier + + # Undo symmetry operations to get full k-point list + scfres_unfold = unfold_bz(scfres) + basis = scfres_unfold.basis + ψ = scfres_unfold.ψ + + # Make wannier directory ... + mkpath(dirname(fileprefix)) + + # Write input file and launch preprocessing + write_w90_win(fileprefix, basis; + num_wann=n_wannier, num_bands=n_bands, wannier_plot, kwargs...) + nnkp = preprocess_call() + + # Files for main wannierization run + write_w90_eig(fileprefix, scfres_unfold.eigenvalues; n_bands) + write_w90_amn(fileprefix, basis, projections, ψ; n_bands) + write_w90_mmn(fileprefix, basis, ψ, nnkp; n_bands) + + # Writing the unk files is expensive (requires FFTs), so only do if needed. + if wannier_plot + write_w90_unk(fileprefix, basis, ψ; n_bands, spin=1) + end +end \ No newline at end of file From c95ba0e99c4c31f355cd9b0c8353e5779b86bff9 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 30 Oct 2023 11:47:26 +0100 Subject: [PATCH 02/29] Add Wannier to docs deps --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index 2fd13363cd..e5d4516dfa 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -24,6 +24,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" +Wannier = "2b19380a-1f7e-4d7d-b1b8-8aa60b3321c9" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" wannier90_jll = "c5400fa0-8d08-52c2-913f-1e3f656c1ce9" From 1941d8ef4c33c3125ad32a5b0cf279c1619850f4 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 30 Oct 2023 13:19:26 +0100 Subject: [PATCH 03/29] Add kpath to example --- examples/wannier.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/wannier.jl b/examples/wannier.jl index 2b0e7235a0..4eec5c7251 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -92,7 +92,8 @@ U = disentangle(wannier_model, max_iter=200); omega(wannier_model) omega(wannier_model, U) -# Build a Wannier interpolation model: (needs a kpath) +# Build a Wannier interpolation model: +kpath = irrfbz_path(model) interp_model = Wannier.InterpModel(wannier_model; kpath=kpath) # And so on... From a9dc3ee842e5375303c39bbd4b40437764eb946e Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 30 Oct 2023 13:45:18 +0100 Subject: [PATCH 04/29] Fix wannier90 runs --- src/external/wannier90.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/external/wannier90.jl b/src/external/wannier90.jl index 0fdc1a79cb..9b5a8dcef8 100644 --- a/src/external/wannier90.jl +++ b/src/external/wannier90.jl @@ -10,7 +10,7 @@ Run the Wannierization procedure with Wannier90. wannier_plot=false, kwargs...) - prefix = basename(fileprefix) + prefix, dir = basename(fileprefix), dirname(fileprefix) # Prepare files write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do From ac9dbb3d37094fb32ae9e0ecf1a887294746af0e Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Wed, 8 Nov 2023 17:32:39 +0100 Subject: [PATCH 05/29] Add Wannier to test dependencies --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ce0ecea82c..3afbfd90c2 100644 --- a/Project.toml +++ b/Project.toml @@ -122,4 +122,4 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" wannier90_jll = "c5400fa0-8d08-52c2-913f-1e3f656c1ce9" [targets] -test = ["Test", "TestItemRunner", "ASEconvert", "Aqua", "AtomsIO", "AtomsIOPython", "CUDA", "CUDA_Runtime_jll", "ComponentArrays", "DoubleFloats", "FiniteDiff", "FiniteDifferences", "GenericLinearAlgebra", "IntervalArithmetic", "JLD2", "JSON3", "Logging", "Plots", "QuadGK", "Random", "KrylovKit", "WriteVTK", "wannier90_jll"] +test = ["Test", "TestItemRunner", "ASEconvert", "Aqua", "AtomsIO", "AtomsIOPython", "CUDA", "CUDA_Runtime_jll", "ComponentArrays", "DoubleFloats", "FiniteDiff", "FiniteDifferences", "GenericLinearAlgebra", "IntervalArithmetic", "JLD2", "JSON3", "Logging", "Plots", "QuadGK", "Random", "KrylovKit", "Wannier", "WriteVTK", "wannier90_jll"] From 1e880419377ce0c0eb476d479eed55e152cecaac Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 13 Nov 2023 18:47:40 +0100 Subject: [PATCH 06/29] Empty after Wannier 0.3.2 release From c9dc3f3b77597ac4cf27eb0c2a5f58dbc357c838 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 13 Nov 2023 18:56:37 +0100 Subject: [PATCH 07/29] Try updating CI to 1.8 to see if that fixes checks --- .github/workflows/ci.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 28ecc30ef2..04c4d18a64 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -37,7 +37,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: - version: '1.7' + version: '1.8' arch: x64 - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 From 502f32e7e22ea132eb9787647be205a86d22a98b Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Tue, 31 Oct 2023 16:31:18 +0100 Subject: [PATCH 08/29] Optimize Amn computation by using less G vectors --- src/external/wannier_shared.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 7ab8741caf..5b03d00b0e 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -46,7 +46,8 @@ struct DoubleGaussianWannierProjection <: WannierProjection end function get_fourier_projection_coefficients(proj::DoubleGaussianWannierProjection, qs) - get_fourier_projection_coefficients(proj.proj1, qs) .+ (proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, qs)) + proj.factor1 .* get_fourier_projection_coefficients(proj.proj1, qs) + .+ proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, qs) end function DoubleGaussianWannierProjection(factor1::Number, center1::AbstractVector{Float64}, factor2::Number, center2::AbstractVector{Float64}) @@ -269,10 +270,10 @@ function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::Abstra # Compute Ak for n in 1:n_wannier proj = projections[n] - # Functions are l^2 normalized in Fourier, in DFTK conventions. - norm_gn_per = norm(get_fourier_projection_coefficients(proj, qs), 2) + gn_per = get_fourier_projection_coefficients(proj, qs[kpt.mapping]) # Fourier coeffs of gn_per for k+G in common with ψk - coeffs_gn_per = get_fourier_projection_coefficients(proj, qs[kpt.mapping]) ./ norm_gn_per + # Functions are l^2 normalized in Fourier, in DFTK conventions. + coeffs_gn_per = gn_per ./ norm(gn_per) # Compute overlap for m in 1:n_bands # TODO Check the ordering of m and n here! From 9abd0f99ceb00554f8d0154c46c1cbae27e56ab9 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 13 Nov 2023 10:29:05 +0100 Subject: [PATCH 09/29] Add hydrogenic initial Wannier projection Co-authored-by: Junfeng Qiao --- src/DFTK.jl | 1 + src/common/hydrogenic.jl | 22 +++++++++++ src/external/wannier_shared.jl | 67 +++++++++++++++++++++++++++++++--- 3 files changed, 84 insertions(+), 6 deletions(-) create mode 100644 src/common/hydrogenic.jl diff --git a/src/DFTK.jl b/src/DFTK.jl index 67dd2927aa..707d32ecf6 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -37,6 +37,7 @@ include("common/versioninfo.jl") include("architecture.jl") include("common/zeros_like.jl") include("common/norm.jl") +include("common/hydrogenic.jl") export PspHgh export PspUpf diff --git a/src/common/hydrogenic.jl b/src/common/hydrogenic.jl new file mode 100644 index 0000000000..a215e087c0 --- /dev/null +++ b/src/common/hydrogenic.jl @@ -0,0 +1,22 @@ +@doc raw""" +Radial functions from solutions of Hydrogenic Schrödinger equation. +Same as Wannier90 user guide Table 3.3. +# Arguments +- `r`: radial grid +- `n`: principal quantum number +- `α`: diffusivity, ``\frac{Z}/{a}`` where ``Z`` is the atomic number and + ``a`` is the Bohr radius. +""" +function radial_hydrogenic(r::AbstractVector{T}, n::Integer, α::Real=1.0) where {T<:Real} + @assert n > 0 + if n == 1 + f = 2 * α^(3/2) * exp.(-α * r) + elseif n == 2 + f = 2^(-3/2) * α^(3/2) * (2 .- α * r) .* exp.(-α * r/2) + elseif n == 3 + f = sqrt(4/27) * α^(3/2) * (1 .- 2/3 * α * r .+ 2/27 * α^2 * r.^2) .* exp.(-α * r/3) + else + error("n = $n is not supported") + end + f +end \ No newline at end of file diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 5b03d00b0e..b7bab1feb4 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -25,10 +25,15 @@ struct GaussianWannierProjection <: WannierProjection center::AbstractVector{Float64} end -function get_fourier_projection_coefficients(proj::GaussianWannierProjection, qs) +function get_fourier_projection_coefficients(proj::GaussianWannierProjection, basis::PlaneWaveBasis, qs) # associate a center with the fourier transform of the corresponding gaussian # TODO: what is the normalization here? - [exp( 2π*(-im*dot(q, proj.center) - dot(q, q) / 4) ) for q in qs] + + model = basis.model + [begin + q_cart = recip_vector_red_to_cart(model, q) + exp( 2π*(-im*dot(q, proj.center) - dot(q_cart, q_cart) / 4) ) + end for q in qs] end function random_gaussian_projection() @@ -45,9 +50,9 @@ struct DoubleGaussianWannierProjection <: WannierProjection proj2::GaussianWannierProjection end -function get_fourier_projection_coefficients(proj::DoubleGaussianWannierProjection, qs) - proj.factor1 .* get_fourier_projection_coefficients(proj.proj1, qs) - .+ proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, qs) +function get_fourier_projection_coefficients(proj::DoubleGaussianWannierProjection, basis::PlaneWaveBasis, qs) + proj.factor1 .* get_fourier_projection_coefficients(proj.proj1, basis, qs) + .+ proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, basis, qs) end function DoubleGaussianWannierProjection(factor1::Number, center1::AbstractVector{Float64}, factor2::Number, center2::AbstractVector{Float64}) @@ -63,6 +68,56 @@ function opposite_gaussians_projection(center1::AbstractVector{Float64}, center2 DoubleGaussianWannierProjection(1, center1, -1, center2) end +@doc raw""" +A hydrogenic initial guess. + +`α` is the diffusivity, ``\frac{Z}/{a}`` where ``Z`` is the atomic number and + ``a`` is the Bohr radius. +""" +struct HydrogenicWannierProjection <: WannierProjection + center::AbstractVector{Float64} + n::Integer + l::Integer + m::Integer + α::Real +end + +function get_fourier_projection_coefficients(proj::HydrogenicWannierProjection, basis::PlaneWaveBasis, qs) + # TODO: Performance can probably be improved a lot here. + # TODO: Some of this logic could be used for the pswfc from the UPF PSPs. + + # this is same as QE pw2wannier90, r = exp(x) / α + xmin = -6.0 + dx = 0.025 + rmax = 10.0 + n_r = round(Int, (log(rmax) - xmin) / dx) + 1 + x = range(; start=xmin, length=n_r, step=dx) + # rgrid depends on α + r = exp.(x) ./ proj.α + dr = r .* dx + R = radial_hydrogenic(r, proj.n, proj.α) + r2_R_dr = r.^2 .* R .* dr + + model = basis.model + + [begin + q_cart = recip_vector_red_to_cart(model, q) + + qnorm = norm(q_cart) + radial_part = 0 + for (ir, rval) in enumerate(r) + @inbounds radial_part += r2_R_dr[ir] * sphericalbesselj_fast(proj.l, qnorm * rval) + end + + # both q and proj.center in reduced coordinates + center_offset = exp(-im*2π*dot(q, proj.center)) + # without the 4π normalization because the matrix will be orthonormalized anyway + angular_part = ylm_real(proj.l, proj.m, q_cart) * (-im)^proj.l + + center_offset * angular_part * radial_part + end for q in qs] +end + """ Write a win file at the indicated prefix. Parameters to Wannier90 can be added as kwargs: e.g. `num_iter=500`. @@ -270,7 +325,7 @@ function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::Abstra # Compute Ak for n in 1:n_wannier proj = projections[n] - gn_per = get_fourier_projection_coefficients(proj, qs[kpt.mapping]) + gn_per = get_fourier_projection_coefficients(proj, basis, qs[kpt.mapping]) # Fourier coeffs of gn_per for k+G in common with ψk # Functions are l^2 normalized in Fourier, in DFTK conventions. coeffs_gn_per = gn_per ./ norm(gn_per) From 4c2c6dd3b29a6d34604e8c366d14b05df28246df Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Tue, 14 Nov 2023 15:56:09 +0100 Subject: [PATCH 10/29] Update example to use hydrogenic projections --- examples/wannier.jl | 45 ++++++++++++++++++++------------------------- 1 file changed, 20 insertions(+), 25 deletions(-) diff --git a/examples/wannier.jl b/examples/wannier.jl index 4eec5c7251..115a91e308 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -53,19 +53,14 @@ plot_bandstructure(scfres; kline_density=10) # For a good MLWF, we need to provide initial projections that resemble the expected shape # of the Wannier functions. # Here we will use: -# - 3 bond-centered gaussians for the expected σ bonds -# - 2 atom-centered opposite gaussians (for each atom: a gaussian with the + sign and a gaussian with the - sign) -# for the expected pz orbitals +# - 3 bond-centered 2s hydrogenic orbitals for the expected σ bonds +# - 2 atom-centered 2pz hydrogenic orbitals for the expected π bands using Wannier # Needed to make get_wannier_model available -"""Helper function to produce a starting guess for pz orbitals""" -make_double_gaussian(atomindex) = begin - offset = model.inv_lattice * [0, 0, austrip(0.5u"Å")] - center1 = positions[atomindex] + offset - center2 = positions[atomindex] - offset - DFTK.opposite_gaussians_projection(center1, center2) -end +## Helper functions to produce hydrogenic starting guesses +s_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 0, 0, C.Z) +pz_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 1, 0, C.Z) wannier_model = get_wannier_model(scfres; fileprefix="wannier/graphene", @@ -73,13 +68,13 @@ wannier_model = get_wannier_model(scfres; n_wannier=5, projections=[ ## Note: fractional coordinates for the centers! - ## 3 bond-centered gaussians to imitate σ orbitals - DFTK.GaussianWannierProjection((positions[1] + positions[2]) / 2), - DFTK.GaussianWannierProjection((positions[1] + positions[2] + [0, -1, 0]) / 2), - DFTK.GaussianWannierProjection((positions[1] + positions[2] + [-1, -1, 0]) / 2), - ## 2 double gaussians to imitate pz orbitals - make_double_gaussian(1), - make_double_gaussian(2), + ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds + s_guess((positions[1] + positions[2]) / 2), + s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), + s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), + ## 2 atom-centered 2pz hydrogenic orbitals + pz_guess(positions[1]), + pz_guess(positions[2]), ], dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1) # maximum frozen window, for example 1 eV above Fermi level @@ -102,7 +97,7 @@ interp_model = Wannier.InterpModel(wannier_model; kpath=kpath) # (Delete temporary files when done.) rm("wannier", recursive=true) -# This example assumes that Wannier.jl version 0.3.1 is used, +# This example assumes that Wannier.jl version 0.3.2 is used, # and may need to be updated to accomodate for changes in Wannier.jl. # # Note: Some parameters supported by Wannier90 may have to be passed to Wannier.jl differently, @@ -119,13 +114,13 @@ run_wannier90(scfres; n_wannier=5, projections=[ ## Note: fractional coordinates for the centers! - ## 3 bond-centered gaussians to imitate σ orbitals - DFTK.GaussianWannierProjection((positions[1] + positions[2]) / 2), - DFTK.GaussianWannierProjection((positions[1] + positions[2] + [0, -1, 0]) / 2), - DFTK.GaussianWannierProjection((positions[1] + positions[2] + [-1, -1, 0]) / 2), - ## 2 double gaussians to imitate pz orbitals - make_double_gaussian(1), - make_double_gaussian(2), + ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds + s_guess((positions[1] + positions[2]) / 2), + s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), + s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), + ## 2 atom-centered 2pz hydrogenic orbitals + pz_guess(positions[1]), + pz_guess(positions[2]), ], num_print_cycles=25, num_iter=200, From 1588782faf14c6782a2937f626f0d2dad16decad Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 09:49:05 +0100 Subject: [PATCH 11/29] Ensure that output of radial_hydrogenic is always of type T --- src/common/hydrogenic.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/common/hydrogenic.jl b/src/common/hydrogenic.jl index a215e087c0..a5208ffd16 100644 --- a/src/common/hydrogenic.jl +++ b/src/common/hydrogenic.jl @@ -7,14 +7,15 @@ Same as Wannier90 user guide Table 3.3. - `α`: diffusivity, ``\frac{Z}/{a}`` where ``Z`` is the atomic number and ``a`` is the Bohr radius. """ -function radial_hydrogenic(r::AbstractVector{T}, n::Integer, α::Real=1.0) where {T<:Real} +function radial_hydrogenic(r::AbstractVector{T}, n::Integer, α::Real=one(T)) where {T<:Real} + # T(...) is used for some literals to ensure that the output is of type T @assert n > 0 if n == 1 - f = 2 * α^(3/2) * exp.(-α * r) + f = 2 * α^(3/T(2)) * exp.(-α * r) elseif n == 2 - f = 2^(-3/2) * α^(3/2) * (2 .- α * r) .* exp.(-α * r/2) + f = 2^(-3/T(2)) * α^(3/T(2)) * (2 .- α * r) .* exp.(-α * r/2) elseif n == 3 - f = sqrt(4/27) * α^(3/2) * (1 .- 2/3 * α * r .+ 2/27 * α^2 * r.^2) .* exp.(-α * r/3) + f = sqrt(4/T(27)) * α^(3/T(2)) * (1 .- 2/T(3) * α * r .+ 2/T(27) * α^2 * r.^2) .* exp.(-α * r/3) else error("n = $n is not supported") end From ace1d4e85eb118b209f0ac654cece2b257481895 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 10:00:31 +0100 Subject: [PATCH 12/29] `get_wannier_model` -> `wannier_model` --- examples/wannier.jl | 8 ++++---- src/DFTK.jl | 2 +- src/external/stubs.jl | 4 ++-- src/external/wannier.jl | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/wannier.jl b/examples/wannier.jl index b1be740b4b..927234afbe 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -4,7 +4,7 @@ # [Wannier.jl](https://wannierjl.org) and [Wannier90](http://www.wannier.org/), # in order to compute maximally-localized Wannier functions (MLWFs) # from an initial self consistent field calculation. -# All processes are handled by calling the routine `get_wannier_model` (for Wannier.jl) or `run_wannier90` (for Wannier90). +# All processes are handled by calling the routine `wannier_model` (for Wannier.jl) or `run_wannier90` (for Wannier90). # # !!! warning "No guarantees on Wannier interface" # This code is at an early stage and has so far not been fully tested. @@ -39,7 +39,7 @@ plot_bandstructure(scfres; kline_density=10) # ## Wannierization with Wannier.jl # -# Now we use the `get_wannier_model` routine to generate a Wannier.jl model +# Now we use the `wannier_model` routine to generate a Wannier.jl model # that can be used to perform the wannierization procedure. # For now, this model generation produces file in the Wannier90 convention, # where all files are named with the same prefix and only differ by @@ -56,13 +56,13 @@ plot_bandstructure(scfres; kline_density=10) # - 3 bond-centered 2s hydrogenic orbitals for the expected σ bonds # - 2 atom-centered 2pz hydrogenic orbitals for the expected π bands -using Wannier # Needed to make get_wannier_model available +using Wannier # Needed to make wannier_model available ## Helper functions to produce hydrogenic starting guesses s_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 0, 0, C.Z) pz_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 1, 0, C.Z) -wannier_model = get_wannier_model(scfres; +wannier_model = wannier_model(scfres; fileprefix="wannier/graphene", n_bands=scfres.n_bands_converge, n_wannier=5, diff --git a/src/DFTK.jl b/src/DFTK.jl index d9a8552ccd..32b034e0a8 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -185,7 +185,7 @@ include("pseudo/list_psp.jl") include("pseudo/attach_psp.jl") export atomic_system, periodic_system # Reexport from AtomsBase -export get_wannier_model +export wannier_model export run_wannier90 include("external/atomsbase.jl") include("external/stubs.jl") # Function stubs for conditionally defined methods diff --git a/src/external/stubs.jl b/src/external/stubs.jl index 00c099381d..e1160b894d 100644 --- a/src/external/stubs.jl +++ b/src/external/stubs.jl @@ -15,11 +15,11 @@ The function returns the `fileprefix`. incompatibly in the future. Use at your own risk and please report bugs in case you encounter any. """ -function get_wannier_model end +function wannier_model end """ Wannerize the obtained bands using wannier90. -Same arguments as `get_wannier_model`. +Same arguments as [`wannier_model`](@ref). The function returns the `fileprefix`. !!! warning "Experimental feature" diff --git a/src/external/wannier.jl b/src/external/wannier.jl index e4b0a33a46..d8d0d059d2 100644 --- a/src/external/wannier.jl +++ b/src/external/wannier.jl @@ -19,7 +19,7 @@ end """ Build a Wannier.jl model that can be used for Wannierization. """ -@timing function get_wannier_model(scfres; +@timing function wannier_model(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, projections::AbstractVector{<:WannierProjection}=default_wannier_centres(n_wannier), From bfdb38a2edb0784ac6c57fdbba5f46c898ccf299 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 10:03:37 +0100 Subject: [PATCH 13/29] Remove random_gaussian_projection and fix centres typo --- src/external/wannier.jl | 2 +- src/external/wannier90.jl | 2 +- src/external/wannier_shared.jl | 6 +----- 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/external/wannier.jl b/src/external/wannier.jl index d8d0d059d2..46ee91e3bb 100644 --- a/src/external/wannier.jl +++ b/src/external/wannier.jl @@ -22,7 +22,7 @@ Build a Wannier.jl model that can be used for Wannierization. @timing function wannier_model(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector{<:WannierProjection}=default_wannier_centres(n_wannier), + projections::AbstractVector{<:WannierProjection}=default_wannier_centers(n_wannier), fileprefix=joinpath("wannierjl", "wannier"), wannier_plot=false, kwargs...) diff --git a/src/external/wannier90.jl b/src/external/wannier90.jl index 9b5a8dcef8..e1b91cd597 100644 --- a/src/external/wannier90.jl +++ b/src/external/wannier90.jl @@ -5,7 +5,7 @@ Run the Wannierization procedure with Wannier90. @timing function run_wannier90(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector{<:WannierProjection}=default_wannier_centres(n_wannier), + projections::AbstractVector{<:WannierProjection}=default_wannier_centers(n_wannier), fileprefix=joinpath("wannier90", "wannier"), wannier_plot=false, kwargs...) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 78b8732ab1..926cde10c8 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -36,10 +36,6 @@ function get_fourier_projection_coefficients(proj::GaussianWannierProjection, ba end for q in qs] end -function random_gaussian_projection() - GaussianWannierProjection(rand(3)) -end - """ A sum of two Gaussians. Can be used as an approximation of a p-like orbital. """ @@ -366,7 +362,7 @@ end Default random Gaussian guess for maximally-localised wannier functions generated in reduced coordinates. """ -default_wannier_centres(n_wannier) = [random_gaussian_projection() for _ in 1:n_wannier] +default_wannier_centers(n_wannier) = [GaussianWannierProjection(rand(3)) for _ in 1:n_wannier] """ Shared file writing code for Wannier.jl and Wannier90. From 34603d71308d0b47bc99adfd45bc3a7bb02dc807 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 14:31:45 +0100 Subject: [PATCH 14/29] Move Wannier compat to extensions, remove WannierProjection abstract type --- Project.toml | 4 + examples/wannier.jl | 4 +- .../wannier90.jl => ext/DFTKWannier90Ext.jl | 14 +++- .../wannier.jl => ext/DFTKWannierExt.jl | 13 +++- src/DFTK.jl | 6 -- src/external/wannier_shared.jl | 75 ++++++------------- 6 files changed, 50 insertions(+), 66 deletions(-) rename src/external/wannier90.jl => ext/DFTKWannier90Ext.jl (61%) rename src/external/wannier.jl => ext/DFTKWannierExt.jl (76%) diff --git a/Project.toml b/Project.toml index 8b4e990fcb..23e07234be 100644 --- a/Project.toml +++ b/Project.toml @@ -50,10 +50,14 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [weakdeps] JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +Wannier = "2b19380a-1f7e-4d7d-b1b8-8aa60b3321c9" +wannier90_jll = "c5400fa0-8d08-52c2-913f-1e3f656c1ce9" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [extensions] DFTKJSON3Ext = "JSON3" +DFTKWannierExt = "Wannier" +DFTKWannier90Ext = "wannier90_jll" DFTKWriteVTK = "WriteVTK" [compat] diff --git a/examples/wannier.jl b/examples/wannier.jl index 927234afbe..9fd869dcc9 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -56,13 +56,13 @@ plot_bandstructure(scfres; kline_density=10) # - 3 bond-centered 2s hydrogenic orbitals for the expected σ bonds # - 2 atom-centered 2pz hydrogenic orbitals for the expected π bands -using Wannier # Needed to make wannier_model available +using Wannier # Needed to make Wannier.Model available ## Helper functions to produce hydrogenic starting guesses s_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 0, 0, C.Z) pz_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 1, 0, C.Z) -wannier_model = wannier_model(scfres; +wannier_model = Wannier.Model(scfres; fileprefix="wannier/graphene", n_bands=scfres.n_bands_converge, n_wannier=5, diff --git a/src/external/wannier90.jl b/ext/DFTKWannier90Ext.jl similarity index 61% rename from src/external/wannier90.jl rename to ext/DFTKWannier90Ext.jl index e1b91cd597..cd0fa006cc 100644 --- a/src/external/wannier90.jl +++ b/ext/DFTKWannier90Ext.jl @@ -1,11 +1,15 @@ +module DFTKWannier90Ext + +using DFTK +import wannier90_jll """ Run the Wannierization procedure with Wannier90. """ -@timing function run_wannier90(scfres; +@DFTK.timing function run_wannier90(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector{<:WannierProjection}=default_wannier_centers(n_wannier), + projections::AbstractVector=DFTK.default_wannier_centers(n_wannier), fileprefix=joinpath("wannier90", "wannier"), wannier_plot=false, kwargs...) @@ -13,12 +17,14 @@ Run the Wannierization procedure with Wannier90. prefix, dir = basename(fileprefix), dirname(fileprefix) # Prepare files - write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do + DFTK.write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do wannier90_jll.wannier90(exe -> run(Cmd(`$exe -pp $prefix`; dir))) - read_w90_nnkp(fileprefix) + DFTK.read_w90_nnkp(fileprefix) end # Run Wannierisation procedure @timing "Wannierization" wannier90_jll.wannier90(exe -> run(Cmd(`$exe $prefix`; dir))) fileprefix +end + end \ No newline at end of file diff --git a/src/external/wannier.jl b/ext/DFTKWannierExt.jl similarity index 76% rename from src/external/wannier.jl rename to ext/DFTKWannierExt.jl index 46ee91e3bb..c9b0e2ee7f 100644 --- a/src/external/wannier.jl +++ b/ext/DFTKWannierExt.jl @@ -1,3 +1,8 @@ +module DFTKWannierExt + +using DFTK +import Wannier + """ Use Wannier.jl to read a .win file and produce the nnkp data (i.e. the b-vectors used in the Mmn matrix). """ @@ -19,18 +24,20 @@ end """ Build a Wannier.jl model that can be used for Wannierization. """ -@timing function wannier_model(scfres; +@DFTK.timing function Wannier.Model(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector{<:WannierProjection}=default_wannier_centers(n_wannier), + projections::AbstractVector=DFTK.default_wannier_centers(n_wannier), fileprefix=joinpath("wannierjl", "wannier"), wannier_plot=false, kwargs...) # Write the files - write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do + DFTK.write_wannier90_files(scfres; n_bands, n_wannier, projections, fileprefix, wannier_plot, kwargs...) do get_nnkpt_from_wannier(fileprefix) end # Read Wannier.jl model Wannier.read_w90(fileprefix) +end + end \ No newline at end of file diff --git a/src/DFTK.jl b/src/DFTK.jl index 32b034e0a8..a4a6eb267f 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -240,12 +240,6 @@ function __init__() # TODO Keep these requires for now as there are open PRs changing these files. @require JLD2="033835bb-8acc-5ee8-8aae-3f567f8a3819" include("external/jld2io.jl") @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plotting.jl") - @require wannier90_jll="c5400fa0-8d08-52c2-913f-1e3f656c1ce9" begin - include("external/wannier90.jl") - end - @require Wannier="2b19380a-1f7e-4d7d-b1b8-8aa60b3321c9" begin - include("external/wannier.jl") - end # TODO Move out into extension module @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 926cde10c8..149e9dfd54 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -3,29 +3,14 @@ using Dates -@doc raw""" -Base type for Wannier projections. - -Wannierization searchs for a unitary matrix ``U_{m_n}``. -As a starting point for the search, we can provide an initial guess function ``g`` -for the shape of the Wannier functions, based on what we expect from knowledge of the problem or physical intuition. -This starting matrix is called ``[A_k]_{m,n}``, and is computed as follows: -``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle``. -The matrix will be orthonormalized by the chosen Wannier program, we don't need to do so ourselves. -""" -abstract type WannierProjection end - -# Dispatch function for WannierProjection implementations, evaluating g in Fourier space. -function get_fourier_projection_coefficients end - """ A Gaussian-shaped initial guess. Can be used as an approximation of an s- or σ-like orbital. """ -struct GaussianWannierProjection <: WannierProjection +struct GaussianWannierProjection center::AbstractVector{Float64} end -function get_fourier_projection_coefficients(proj::GaussianWannierProjection, basis::PlaneWaveBasis, qs) +function (proj::GaussianWannierProjection)(basis::PlaneWaveBasis, qs) # associate a center with the fourier transform of the corresponding gaussian # TODO: what is the normalization here? @@ -37,31 +22,12 @@ function get_fourier_projection_coefficients(proj::GaussianWannierProjection, ba end """ -A sum of two Gaussians. Can be used as an approximation of a p-like orbital. +A p-like initial guess, using the difference of two Gaussians with different centers. """ -struct DoubleGaussianWannierProjection <: WannierProjection - factor1::Number - proj1::GaussianWannierProjection - factor2::Number - proj2::GaussianWannierProjection -end - -function get_fourier_projection_coefficients(proj::DoubleGaussianWannierProjection, basis::PlaneWaveBasis, qs) - proj.factor1 .* get_fourier_projection_coefficients(proj.proj1, basis, qs) - .+ proj.factor2 .* get_fourier_projection_coefficients(proj.proj2, basis, qs) -end - -function DoubleGaussianWannierProjection(factor1::Number, center1::AbstractVector{Float64}, factor2::Number, center2::AbstractVector{Float64}) - DoubleGaussianWannierProjection( - factor1, - GaussianWannierProjection(center1), - factor2, - GaussianWannierProjection(center2), - ) -end - -function opposite_gaussians_projection(center1::AbstractVector{Float64}, center2::AbstractVector{Float64}) - DoubleGaussianWannierProjection(1, center1, -1, center2) +function opposite_gaussians_projection(center1::AbstractVector, center2::AbstractVector) + function inner(basis, qs) + GaussianWannierProjection(center1)(basis, qs) - GaussianWannierProjection(center2)(basis, qs) + end end @doc raw""" @@ -70,7 +36,7 @@ A hydrogenic initial guess. `α` is the diffusivity, ``\frac{Z}/{a}`` where ``Z`` is the atomic number and ``a`` is the Bohr radius. """ -struct HydrogenicWannierProjection <: WannierProjection +struct HydrogenicWannierProjection center::AbstractVector{Float64} n::Integer l::Integer @@ -78,7 +44,7 @@ struct HydrogenicWannierProjection <: WannierProjection α::Real end -function get_fourier_projection_coefficients(proj::HydrogenicWannierProjection, basis::PlaneWaveBasis, qs) +function (proj::HydrogenicWannierProjection)(basis::PlaneWaveBasis, qs) # TODO: Performance can probably be improved a lot here. # TODO: Some of this logic could be used for the pswfc from the UPF PSPs. @@ -298,20 +264,27 @@ end end @doc raw""" -Compute the matrix ``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle`` +Compute the starting matrix for Wannierization. -``g^{per}_n`` are periodized gaussians whose respective centers are given as an - (num_bands,1) array [ [center 1], ... ]. +Wannierization searches for a unitary matrix ``U_{m_n}``. +As a starting point for the search, we can provide an initial guess function ``g`` +for the shape of the Wannier functions, based on what we expect from knowledge of the problem or physical intuition. +This starting matrix is called ``[A_k]_{m,n}``, and is computed as follows: +``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle``. +The matrix will be orthonormalized by the chosen Wannier program, we don't need to do so ourselves. Centers are to be given in lattice coordinates and G_vectors in reduced coordinates. The dot product is computed in the Fourier space. Given an orbital ``g_n``, the periodized orbital is defined by : ``g^{per}_n = \sum\limits_{R \in {\rm lattice}} g_n( \cdot - R)``. -The Fourier coefficient of ``g^{per}_n`` at any G +The Fourier coefficient of ``g^{per}_n`` at any G is given by the value of the Fourier transform of ``g_n`` in G. + +Each projection is a callable object that accepts the basis and some qpoints as an argument, +and returns the Fourier transform of ``g_n`` at the qpoints. """ -function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::AbstractVector{<:WannierProjection}, n_bands) +function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::AbstractVector, n_bands) n_wannier = length(projections) # TODO This function should be improved in performance @@ -321,7 +294,7 @@ function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::Abstra # Compute Ak for n in 1:n_wannier proj = projections[n] - gn_per = get_fourier_projection_coefficients(proj, basis, qs[kpt.mapping]) + gn_per = proj(basis, qs[kpt.mapping]) # Fourier coeffs of gn_per for k+G in common with ψk # Functions are l^2 normalized in Fourier, in DFTK conventions. coeffs_gn_per = gn_per ./ norm(gn_per) @@ -338,7 +311,7 @@ end @timing function write_w90_amn( fileprefix::String, basis::PlaneWaveBasis, - projections::AbstractVector{<:WannierProjection}, + projections::AbstractVector, ψ; n_bands) open(fileprefix * ".amn", "w") do fp @@ -370,7 +343,7 @@ Shared file writing code for Wannier.jl and Wannier90. @timing function write_wannier90_files(preprocess_call::Function, scfres; n_bands, n_wannier, - projections::AbstractVector{<:WannierProjection}, + projections::AbstractVector, fileprefix, wannier_plot, kwargs...) From b5fecedf3961cb94ff8bea37e488761c50f1aa46 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 15:09:15 +0100 Subject: [PATCH 15/29] Fix wannier90 ext precompilation --- ext/DFTKWannier90Ext.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/DFTKWannier90Ext.jl b/ext/DFTKWannier90Ext.jl index cd0fa006cc..97bb860568 100644 --- a/ext/DFTKWannier90Ext.jl +++ b/ext/DFTKWannier90Ext.jl @@ -23,7 +23,7 @@ Run the Wannierization procedure with Wannier90. end # Run Wannierisation procedure - @timing "Wannierization" wannier90_jll.wannier90(exe -> run(Cmd(`$exe $prefix`; dir))) + @DFTK.timing "Wannierization" wannier90_jll.wannier90(exe -> run(Cmd(`$exe $prefix`; dir))) fileprefix end From 6abe0648e482649f43c08fd31dd7fac25e5c54f8 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 15:45:22 +0100 Subject: [PATCH 16/29] Add test for type-stability of radial_hydrogenic --- test/hydrogenic.jl | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 test/hydrogenic.jl diff --git a/test/hydrogenic.jl b/test/hydrogenic.jl new file mode 100644 index 0000000000..9ea629400f --- /dev/null +++ b/test/hydrogenic.jl @@ -0,0 +1,8 @@ +@testitem "Type-stability of radial_hydrogenic" begin + using DFTK: radial_hydrogenic + + for n = [1, 2, 3] + @inferred radial_hydrogenic([1.0, 1.0], n) + @inferred radial_hydrogenic([1f0, 1f0], n) + end +end \ No newline at end of file From 5cf8726d1d8dd9082e80d7436e6d8b3afe41cb6e Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 15:48:42 +0100 Subject: [PATCH 17/29] Don't use begin and end in list comprehensions --- src/external/wannier_shared.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 149e9dfd54..856997367e 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -15,10 +15,10 @@ function (proj::GaussianWannierProjection)(basis::PlaneWaveBasis, qs) # TODO: what is the normalization here? model = basis.model - [begin + map(qs) do q q_cart = recip_vector_red_to_cart(model, q) exp( 2π*(-im*dot(q, proj.center) - dot(q_cart, q_cart) / 4) ) - end for q in qs] + end end """ @@ -62,7 +62,7 @@ function (proj::HydrogenicWannierProjection)(basis::PlaneWaveBasis, qs) model = basis.model - [begin + map(qs) do q q_cart = recip_vector_red_to_cart(model, q) qnorm = norm(q_cart) @@ -77,7 +77,7 @@ function (proj::HydrogenicWannierProjection)(basis::PlaneWaveBasis, qs) angular_part = ylm_real(proj.l, proj.m, q_cart) * (-im)^proj.l center_offset * angular_part * radial_part - end for q in qs] + end end """ From 26e321ed1529474924441f1c5842ed2369745854 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 16:16:26 +0100 Subject: [PATCH 18/29] Add basic Wannier.jl test --- Project.toml | 2 ++ test/external/wannier.jl | 68 ++++++++++++++++++++++++++++++++++++++ test/external/wannier90.jl | 32 ------------------ 3 files changed, 70 insertions(+), 32 deletions(-) create mode 100644 test/external/wannier.jl delete mode 100644 test/external/wannier90.jl diff --git a/Project.toml b/Project.toml index 23e07234be..95e8e7309d 100644 --- a/Project.toml +++ b/Project.toml @@ -106,6 +106,8 @@ Statistics = "1" TimerOutputs = "0.5.12" Unitful = "1" UnitfulAtomic = "1" +Wannier = "0.3.2" +wannier90_jll = "3.1" WriteVTK = "1" julia = "1.9" diff --git a/test/external/wannier.jl b/test/external/wannier.jl new file mode 100644 index 0000000000..ce0f3acecd --- /dev/null +++ b/test/external/wannier.jl @@ -0,0 +1,68 @@ +@testitem "Test run_wannier90" tags=[:dont_test_mpi, :dont_test_windows] setup=[TestCases] begin + using DFTK + using wannier90_jll + silicon = TestCases.silicon + + model = model_LDA(silicon.lattice, silicon.atoms, silicon.positions) + basis = PlaneWaveBasis(model; Ecut=5, kgrid=[4, 4, 4], kshift=[1, 1, 1]/2) + nbandsalg = AdaptiveBands(model; n_bands_converge=12) + scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) + + fileprefix = "wannier90_outputs/Si" + run_wannier90(scfres; fileprefix, + n_wannier=8, bands_plot=true, + num_print_cycles=50, num_iter=500, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_num_iter=120, + dis_mix_ratio=1.0, + wannier_plot=true) + + @test isfile("wannier90_outputs/Si.amn") + @test isfile("wannier90_outputs/Si.chk") + @test isfile("wannier90_outputs/Si.eig") + @test isfile("wannier90_outputs/Si.mmn") + @test isfile("wannier90_outputs/Si.nnkp") + @test isfile("wannier90_outputs/Si.win") + @test isfile("wannier90_outputs/Si.wout") + @test !isfile("wannier90_outputs/Si.werr") + + # remove produced files + rm("wannier90_outputs", recursive=true) +end + +@testitem "Test calling Wannier.Model with scfres" tags=[:dont_test_mpi] setup=[TestCases] begin + using DFTK + using Wannier + silicon = TestCases.silicon + + model = model_LDA(silicon.lattice, silicon.atoms, silicon.positions) + basis = PlaneWaveBasis(model; Ecut=5, kgrid=[4, 4, 4], kshift=[1, 1, 1]/2) + nbandsalg = AdaptiveBands(model; n_bands_converge=12) + scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) + + fileprefix = "wannierjl_outputs/Si" + wannier_model = Wannier.Model(scfres; fileprefix, + n_wannier=8, bands_plot=true, + num_print_cycles=50, num_iter=500, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_mix_ratio=1.0, + wannier_plot=true) + + wannier_model.U .= disentangle(wannier_model; max_iter=500) + + # for now the Wannier.jl compat writes the amn, eig, mmn and win files + @test isfile("wannierjl_outputs/Si.amn") + @test !isfile("wannierjl_outputs/Si.chk") + @test isfile("wannierjl_outputs/Si.eig") + @test isfile("wannierjl_outputs/Si.mmn") + @test !isfile("wannierjl_outputs/Si.nnkp") + @test isfile("wannierjl_outputs/Si.win") + @test !isfile("wannierjl_outputs/Si.wout") # NO .wout file, the output is in the `wannier_model` + @test !isfile("wannierjl_outputs/Si.werr") + + # remove produced files + rm("wannierjl_outputs", recursive=true) +end + diff --git a/test/external/wannier90.jl b/test/external/wannier90.jl deleted file mode 100644 index b46094dd51..0000000000 --- a/test/external/wannier90.jl +++ /dev/null @@ -1,32 +0,0 @@ -@testitem "Test run_wannier90" tags=[:dont_test_mpi, :dont_test_windows] setup=[TestCases] begin - using DFTK - using wannier90_jll - silicon = TestCases.silicon - - model = model_LDA(silicon.lattice, silicon.atoms, silicon.positions) - basis = PlaneWaveBasis(model; Ecut=5, kgrid=[4, 4, 4], kshift=[1, 1, 1]/2) - nbandsalg = AdaptiveBands(model; n_bands_converge=12) - scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) - - fileprefix = "wannier90_outputs/Si" - run_wannier90(scfres; fileprefix, - n_wannier=8, bands_plot=true, - num_print_cycles=50, num_iter=500, - dis_win_max=17.185257, - dis_froz_max=6.8318033, - dis_num_iter=120, - dis_mix_ratio=1.0, - wannier_plot=true) - - @test isfile("wannier90_outputs/Si.amn") - @test isfile("wannier90_outputs/Si.chk") - @test isfile("wannier90_outputs/Si.eig") - @test isfile("wannier90_outputs/Si.mmn") - @test isfile("wannier90_outputs/Si.nnkp") - @test isfile("wannier90_outputs/Si.win") - @test isfile("wannier90_outputs/Si.wout") - @test !isfile("wannier90_outputs/Si.werr") - - # remove produced files - rm("wannier90_outputs", recursive=true) -end From 59477b65b76fa4a803054c347eff886998b3d0cb Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 16:21:11 +0100 Subject: [PATCH 19/29] Place kpt parameter right after basis --- src/external/wannier_shared.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 856997367e..6396b7349f 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -284,7 +284,7 @@ is given by the value of the Fourier transform of ``g_n`` in G. Each projection is a callable object that accepts the basis and some qpoints as an argument, and returns the Fourier transform of ``g_n`` at the qpoints. """ -function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::AbstractVector, n_bands) +function compute_amn_kpoint(basis::PlaneWaveBasis, kpt, ψk, projections::AbstractVector, n_bands) n_wannier = length(projections) # TODO This function should be improved in performance @@ -300,7 +300,6 @@ function compute_amn_kpoint(basis::PlaneWaveBasis, ψk, kpt, projections::Abstra coeffs_gn_per = gn_per ./ norm(gn_per) # Compute overlap for m in 1:n_bands - # TODO Check the ordering of m and n here! Ak[m, n] = dot(ψk[:, m], coeffs_gn_per) end end @@ -319,7 +318,7 @@ end println(fp, "$n_bands $(length(basis.kpoints)) $(length(projections))") for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints)) - Ak = compute_amn_kpoint(basis, ψk, kpt, projections, n_bands) + Ak = compute_amn_kpoint(basis, kpt, ψk, projections, n_bands) for n in 1:size(Ak, 2) for m in 1:size(Ak, 1) @printf(fp, "%3i %3i %3i %22.18f %22.18f \n", From 80e811b3edc1e76a4a0f2722277eaf4a38b25e2d Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 19:54:51 +0100 Subject: [PATCH 20/29] Remove unused wannier_model stub --- src/DFTK.jl | 1 - src/external/stubs.jl | 19 +++---------------- 2 files changed, 3 insertions(+), 17 deletions(-) diff --git a/src/DFTK.jl b/src/DFTK.jl index a4a6eb267f..8d135e2cc9 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -185,7 +185,6 @@ include("pseudo/list_psp.jl") include("pseudo/attach_psp.jl") export atomic_system, periodic_system # Reexport from AtomsBase -export wannier_model export run_wannier90 include("external/atomsbase.jl") include("external/stubs.jl") # Function stubs for conditionally defined methods diff --git a/src/external/stubs.jl b/src/external/stubs.jl index e1160b894d..578b0a6726 100644 --- a/src/external/stubs.jl +++ b/src/external/stubs.jl @@ -1,27 +1,14 @@ # Stubs for conditionally defined functions """ -Build a Wannier.jl model for the obtained bands. By default all converged +Wannerize the obtained bands using wannier90. By default all converged bands from the `scfres` are employed (change with `n_bands` kwargs) and `n_wannier = n_bands` wannier functions are computed. -Random Gaussians are used as guesses by default, can be changed using the `projections` kwarg. -All keyword arguments supported by +Random Gaussians are used as guesses by default, can be changed +using the `projections` kwarg. All keyword arguments supported by Wannier90 for the disentanglement may be added as keyword arguments. The function returns the `fileprefix`. -!!! warning "Experimental feature" - Currently this is an experimental feature, which has not yet been tested - to full depth. The interface is considered unstable and may change - incompatibly in the future. Use at your own risk and please report bugs - in case you encounter any. -""" -function wannier_model end - -""" -Wannerize the obtained bands using wannier90. -Same arguments as [`wannier_model`](@ref). -The function returns the `fileprefix`. - !!! warning "Experimental feature" Currently this is an experimental feature, which has not yet been tested to full depth. The interface is considered unstable and may change From 62a6069e1f541bb896c8aedd8491c98c3583fbf5 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 20:03:08 +0100 Subject: [PATCH 21/29] Update wannier example a bit --- examples/wannier.jl | 41 +++++++++++++++++++---------------------- 1 file changed, 19 insertions(+), 22 deletions(-) diff --git a/examples/wannier.jl b/examples/wannier.jl index 9fd869dcc9..64ea85b7b6 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -58,24 +58,30 @@ plot_bandstructure(scfres; kline_density=10) using Wannier # Needed to make Wannier.Model available -## Helper functions to produce hydrogenic starting guesses +# From chemical intuition, we know that the bonds with the lowest energy are: +# - the 3 σ bonds, +# - the π and π* bonds. +# We provide relevant initial projections to help Wannierization +# converge to functions with a similar shape. s_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 0, 0, C.Z) pz_guess(center) = DFTK.HydrogenicWannierProjection(center, 2, 1, 0, C.Z) - +projections = [ + ## Note: fractional coordinates for the centers! + ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds + s_guess((positions[1] + positions[2]) / 2), + s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), + s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), + ## 2 atom-centered 2pz hydrogenic orbitals + pz_guess(positions[1]), + pz_guess(positions[2]), +] + +# Wannierize: wannier_model = Wannier.Model(scfres; fileprefix="wannier/graphene", n_bands=scfres.n_bands_converge, n_wannier=5, - projections=[ - ## Note: fractional coordinates for the centers! - ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds - s_guess((positions[1] + positions[2]) / 2), - s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), - s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), - ## 2 atom-centered 2pz hydrogenic orbitals - pz_guess(positions[1]), - pz_guess(positions[2]), - ], + projections, dis_froz_max=ustrip(auconvert(u"eV", scfres.εF))+1) # maximum frozen window, for example 1 eV above Fermi level # Once we have the `wannier_model`, we can use the functions in the Wannier.jl package: @@ -112,16 +118,7 @@ using wannier90_jll # Needed to make run_wannier90 available run_wannier90(scfres; fileprefix="wannier/graphene", n_wannier=5, - projections=[ - ## Note: fractional coordinates for the centers! - ## 3 bond-centered 2s hydrogenic orbitals to imitate σ bonds - s_guess((positions[1] + positions[2]) / 2), - s_guess((positions[1] + positions[2] + [0, -1, 0]) / 2), - s_guess((positions[1] + positions[2] + [-1, -1, 0]) / 2), - ## 2 atom-centered 2pz hydrogenic orbitals - pz_guess(positions[1]), - pz_guess(positions[2]), - ], + projections, num_print_cycles=25, num_iter=200, ## From 43f52f4419c5cb3b6b44c0246cfdc39dffe1c700 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 20:29:29 +0100 Subject: [PATCH 22/29] Use tmpdir for wannierization tests --- test/external/wannier.jl | 79 +++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 41 deletions(-) diff --git a/test/external/wannier.jl b/test/external/wannier.jl index ce0f3acecd..921252fbbe 100644 --- a/test/external/wannier.jl +++ b/test/external/wannier.jl @@ -8,27 +8,26 @@ nbandsalg = AdaptiveBands(model; n_bands_converge=12) scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) - fileprefix = "wannier90_outputs/Si" - run_wannier90(scfres; fileprefix, - n_wannier=8, bands_plot=true, - num_print_cycles=50, num_iter=500, - dis_win_max=17.185257, - dis_froz_max=6.8318033, - dis_num_iter=120, - dis_mix_ratio=1.0, - wannier_plot=true) + mktempdir() do tmpdir + run_wannier90(scfres; + fileprefix="$tmpdir/Si", + n_wannier=8, bands_plot=true, + num_print_cycles=50, num_iter=500, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_num_iter=120, + dis_mix_ratio=1.0, + wannier_plot=true) - @test isfile("wannier90_outputs/Si.amn") - @test isfile("wannier90_outputs/Si.chk") - @test isfile("wannier90_outputs/Si.eig") - @test isfile("wannier90_outputs/Si.mmn") - @test isfile("wannier90_outputs/Si.nnkp") - @test isfile("wannier90_outputs/Si.win") - @test isfile("wannier90_outputs/Si.wout") - @test !isfile("wannier90_outputs/Si.werr") - - # remove produced files - rm("wannier90_outputs", recursive=true) + @test isfile("$tmpdir/Si.amn") + @test isfile("$tmpdir/Si.chk") + @test isfile("$tmpdir/Si.eig") + @test isfile("$tmpdir/Si.mmn") + @test isfile("$tmpdir/Si.nnkp") + @test isfile("$tmpdir/Si.win") + @test isfile("$tmpdir/Si.wout") + @test !isfile("$tmpdir/Si.werr") + end end @testitem "Test calling Wannier.Model with scfres" tags=[:dont_test_mpi] setup=[TestCases] begin @@ -41,28 +40,26 @@ end nbandsalg = AdaptiveBands(model; n_bands_converge=12) scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) - fileprefix = "wannierjl_outputs/Si" - wannier_model = Wannier.Model(scfres; fileprefix, - n_wannier=8, bands_plot=true, - num_print_cycles=50, num_iter=500, - dis_win_max=17.185257, - dis_froz_max=6.8318033, - dis_mix_ratio=1.0, - wannier_plot=true) - - wannier_model.U .= disentangle(wannier_model; max_iter=500) + mktempdir() do tmpdir + wannier_model = Wannier.Model(scfres; + fileprefix="$tmpdir/Si", + n_wannier=8, bands_plot=true, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_mix_ratio=1.0, + wannier_plot=true) - # for now the Wannier.jl compat writes the amn, eig, mmn and win files - @test isfile("wannierjl_outputs/Si.amn") - @test !isfile("wannierjl_outputs/Si.chk") - @test isfile("wannierjl_outputs/Si.eig") - @test isfile("wannierjl_outputs/Si.mmn") - @test !isfile("wannierjl_outputs/Si.nnkp") - @test isfile("wannierjl_outputs/Si.win") - @test !isfile("wannierjl_outputs/Si.wout") # NO .wout file, the output is in the `wannier_model` - @test !isfile("wannierjl_outputs/Si.werr") + wannier_model.U .= disentangle(wannier_model; max_iter=500) - # remove produced files - rm("wannierjl_outputs", recursive=true) + # for now the Wannier.jl compat writes the amn, eig, mmn and win files + @test isfile("$tmpdir/Si.amn") + @test !isfile("$tmpdir/Si.chk") + @test isfile("$tmpdir/Si.eig") + @test isfile("$tmpdir/Si.mmn") + @test !isfile("$tmpdir/Si.nnkp") + @test isfile("$tmpdir/Si.win") + @test !isfile("$tmpdir/Si.wout") # NO .wout file, the output is in the `wannier_model` + @test !isfile("$tmpdir/Si.werr") + end end From f9b80c2a2d6c6e546be3c085a7309c39654a6fac Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 20:36:03 +0100 Subject: [PATCH 23/29] Tweak typing annotations --- ext/DFTKWannier90Ext.jl | 2 +- ext/DFTKWannierExt.jl | 2 +- src/external/wannier_shared.jl | 17 ++++++++--------- 3 files changed, 10 insertions(+), 11 deletions(-) diff --git a/ext/DFTKWannier90Ext.jl b/ext/DFTKWannier90Ext.jl index 97bb860568..1181811a5d 100644 --- a/ext/DFTKWannier90Ext.jl +++ b/ext/DFTKWannier90Ext.jl @@ -9,7 +9,7 @@ Run the Wannierization procedure with Wannier90. @DFTK.timing function run_wannier90(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector=DFTK.default_wannier_centers(n_wannier), + projections=DFTK.default_wannier_centers(n_wannier), fileprefix=joinpath("wannier90", "wannier"), wannier_plot=false, kwargs...) diff --git a/ext/DFTKWannierExt.jl b/ext/DFTKWannierExt.jl index c9b0e2ee7f..305a0a6caa 100644 --- a/ext/DFTKWannierExt.jl +++ b/ext/DFTKWannierExt.jl @@ -27,7 +27,7 @@ Build a Wannier.jl model that can be used for Wannierization. @DFTK.timing function Wannier.Model(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, - projections::AbstractVector=DFTK.default_wannier_centers(n_wannier), + projections=DFTK.default_wannier_centers(n_wannier), fileprefix=joinpath("wannierjl", "wannier"), wannier_plot=false, kwargs...) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index 6396b7349f..af73a7f437 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -7,7 +7,7 @@ using Dates A Gaussian-shaped initial guess. Can be used as an approximation of an s- or σ-like orbital. """ struct GaussianWannierProjection - center::AbstractVector{Float64} + center::AbstractVector end function (proj::GaussianWannierProjection)(basis::PlaneWaveBasis, qs) @@ -37,7 +37,7 @@ A hydrogenic initial guess. ``a`` is the Bohr radius. """ struct HydrogenicWannierProjection - center::AbstractVector{Float64} + center::AbstractVector n::Integer l::Integer m::Integer @@ -61,12 +61,11 @@ function (proj::HydrogenicWannierProjection)(basis::PlaneWaveBasis, qs) r2_R_dr = r.^2 .* R .* dr model = basis.model - map(qs) do q q_cart = recip_vector_red_to_cart(model, q) - qnorm = norm(q_cart) - radial_part = 0 + + radial_part = 0.0 for (ir, rval) in enumerate(r) @inbounds radial_part += r2_R_dr[ir] * sphericalbesselj_fast(proj.l, qnorm * rval) end @@ -284,7 +283,7 @@ is given by the value of the Fourier transform of ``g_n`` in G. Each projection is a callable object that accepts the basis and some qpoints as an argument, and returns the Fourier transform of ``g_n`` at the qpoints. """ -function compute_amn_kpoint(basis::PlaneWaveBasis, kpt, ψk, projections::AbstractVector, n_bands) +function compute_amn_kpoint(basis::PlaneWaveBasis, kpt, ψk, projections, n_bands) n_wannier = length(projections) # TODO This function should be improved in performance @@ -310,7 +309,7 @@ end @timing function write_w90_amn( fileprefix::String, basis::PlaneWaveBasis, - projections::AbstractVector, + projections, ψ; n_bands) open(fileprefix * ".amn", "w") do fp @@ -339,10 +338,10 @@ default_wannier_centers(n_wannier) = [GaussianWannierProjection(rand(3)) for _ i """ Shared file writing code for Wannier.jl and Wannier90. """ -@timing function write_wannier90_files(preprocess_call::Function, scfres; +@timing function write_wannier90_files(preprocess_call, scfres; n_bands, n_wannier, - projections::AbstractVector, + projections, fileprefix, wannier_plot, kwargs...) From 45f6344c366017bd0473e349a6bbcbda162face4 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 20:41:53 +0100 Subject: [PATCH 24/29] for in -> for = --- src/external/wannier_shared.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index af73a7f437..ff57a71b6e 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -127,7 +127,7 @@ function write_w90_win(fileprefix::String, basis::PlaneWaveBasis; path = kpath.paths[1] println(fp, "begin kpoint_path") - for i in 1:length(path)-1 + for i = 1:length(path)-1 A, B = path[i:i+1] # write segment A -> B @printf(fp, "%s %10.6f %10.6f %10.6f ", A, round.(kpath.points[A], digits=5)...) @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpath.points[B], digits=5)...) @@ -203,9 +203,9 @@ end for ik in krange_spin(basis, spin) open(dirname(fileprefix) * (@sprintf "/UNK%05i.%i" ik spin), "w") do fp println(fp, "$(fft_size[1]) $(fft_size[2]) $(fft_size[3]) $ik $n_bands") - for iband in 1:n_bands + for iband = 1:n_bands ψnk_real = ifft(basis, basis.kpoints[ik], @view ψ[ik][:, iband]) - for iz in 1:fft_size[3], iy in 1:fft_size[2], ix in 1:fft_size[1] + for iz = 1:fft_size[3], iy = 1:fft_size[2], ix = 1:fft_size[1] @printf fp "%25.18f %25.18f\n" real(ψnk_real[ix, iy, iz]) imag(ψnk_real[ix, iy, iz]) end end @@ -238,8 +238,8 @@ We use here that: # Compute overlaps # TODO This should be improved ... Mkb = zeros(ComplexF64, (n_bands, n_bands)) - for n in 1:n_bands - for m in 1:n_bands + for n = 1:n_bands + for m = 1:n_bands # Select the coefficient in right order Mkb[m, n] = dot(ψ[ik][iGk, m], ψ[ikpb][iGkpb, n]) end @@ -291,14 +291,14 @@ function compute_amn_kpoint(basis::PlaneWaveBasis, kpt, ψk, projections, n_band Ak = zeros(eltype(ψk), (n_bands, n_wannier)) # Compute Ak - for n in 1:n_wannier + for n = 1:n_wannier proj = projections[n] gn_per = proj(basis, qs[kpt.mapping]) # Fourier coeffs of gn_per for k+G in common with ψk # Functions are l^2 normalized in Fourier, in DFTK conventions. coeffs_gn_per = gn_per ./ norm(gn_per) # Compute overlap - for m in 1:n_bands + for m = 1:n_bands Ak[m, n] = dot(ψk[:, m], coeffs_gn_per) end end @@ -318,8 +318,8 @@ end for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints)) Ak = compute_amn_kpoint(basis, kpt, ψk, projections, n_bands) - for n in 1:size(Ak, 2) - for m in 1:size(Ak, 1) + for n = 1:size(Ak, 2) + for m = 1:size(Ak, 1) @printf(fp, "%3i %3i %3i %22.18f %22.18f \n", m, n, ik, real(Ak[m, n]), imag(Ak[m, n])) end @@ -333,7 +333,7 @@ end Default random Gaussian guess for maximally-localised wannier functions generated in reduced coordinates. """ -default_wannier_centers(n_wannier) = [GaussianWannierProjection(rand(3)) for _ in 1:n_wannier] +default_wannier_centers(n_wannier) = [GaussianWannierProjection(rand(3)) for _ = 1:n_wannier] """ Shared file writing code for Wannier.jl and Wannier90. From 0a3b59b87e0e71fa731a956c103031a40c17b7c2 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 21:40:29 +0100 Subject: [PATCH 25/29] Use a single SCF for the two Wannier tests --- ext/DFTKWannier90Ext.jl | 2 +- test/external/wannier.jl | 103 +++++++++++++++++++-------------------- 2 files changed, 50 insertions(+), 55 deletions(-) diff --git a/ext/DFTKWannier90Ext.jl b/ext/DFTKWannier90Ext.jl index 1181811a5d..6f299f1fe5 100644 --- a/ext/DFTKWannier90Ext.jl +++ b/ext/DFTKWannier90Ext.jl @@ -6,7 +6,7 @@ import wannier90_jll """ Run the Wannierization procedure with Wannier90. """ -@DFTK.timing function run_wannier90(scfres; +@DFTK.timing function DFTK.run_wannier90(scfres; n_bands=scfres.n_bands_converge, n_wannier=n_bands, projections=DFTK.default_wannier_centers(n_wannier), diff --git a/test/external/wannier.jl b/test/external/wannier.jl index 921252fbbe..ae08f76c1c 100644 --- a/test/external/wannier.jl +++ b/test/external/wannier.jl @@ -1,38 +1,7 @@ -@testitem "Test run_wannier90" tags=[:dont_test_mpi, :dont_test_windows] setup=[TestCases] begin - using DFTK - using wannier90_jll - silicon = TestCases.silicon - - model = model_LDA(silicon.lattice, silicon.atoms, silicon.positions) - basis = PlaneWaveBasis(model; Ecut=5, kgrid=[4, 4, 4], kshift=[1, 1, 1]/2) - nbandsalg = AdaptiveBands(model; n_bands_converge=12) - scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) - - mktempdir() do tmpdir - run_wannier90(scfres; - fileprefix="$tmpdir/Si", - n_wannier=8, bands_plot=true, - num_print_cycles=50, num_iter=500, - dis_win_max=17.185257, - dis_froz_max=6.8318033, - dis_num_iter=120, - dis_mix_ratio=1.0, - wannier_plot=true) - - @test isfile("$tmpdir/Si.amn") - @test isfile("$tmpdir/Si.chk") - @test isfile("$tmpdir/Si.eig") - @test isfile("$tmpdir/Si.mmn") - @test isfile("$tmpdir/Si.nnkp") - @test isfile("$tmpdir/Si.win") - @test isfile("$tmpdir/Si.wout") - @test !isfile("$tmpdir/Si.werr") - end -end - -@testitem "Test calling Wannier.Model with scfres" tags=[:dont_test_mpi] setup=[TestCases] begin +@testitem "Test Wannierization" tags=[:dont_test_mpi, :dont_test_windows] setup=[TestCases] begin using DFTK using Wannier + using wannier90_jll silicon = TestCases.silicon model = model_LDA(silicon.lattice, silicon.atoms, silicon.positions) @@ -40,26 +9,52 @@ end nbandsalg = AdaptiveBands(model; n_bands_converge=12) scfres = self_consistent_field(basis; nbandsalg, tol=1e-12) - mktempdir() do tmpdir - wannier_model = Wannier.Model(scfres; - fileprefix="$tmpdir/Si", - n_wannier=8, bands_plot=true, - dis_win_max=17.185257, - dis_froz_max=6.8318033, - dis_mix_ratio=1.0, - wannier_plot=true) - - wannier_model.U .= disentangle(wannier_model; max_iter=500) - - # for now the Wannier.jl compat writes the amn, eig, mmn and win files - @test isfile("$tmpdir/Si.amn") - @test !isfile("$tmpdir/Si.chk") - @test isfile("$tmpdir/Si.eig") - @test isfile("$tmpdir/Si.mmn") - @test !isfile("$tmpdir/Si.nnkp") - @test isfile("$tmpdir/Si.win") - @test !isfile("$tmpdir/Si.wout") # NO .wout file, the output is in the `wannier_model` - @test !isfile("$tmpdir/Si.werr") + @testset begin + @testset "Using wannier90" begin + mktempdir() do tmpdir + run_wannier90(scfres; + fileprefix="$tmpdir/Si", + n_wannier=8, bands_plot=true, + num_print_cycles=50, num_iter=500, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_num_iter=120, + dis_mix_ratio=1.0, + wannier_plot=true) + + @test isfile("$tmpdir/Si.amn") + @test isfile("$tmpdir/Si.chk") + @test isfile("$tmpdir/Si.eig") + @test isfile("$tmpdir/Si.mmn") + @test isfile("$tmpdir/Si.nnkp") + @test isfile("$tmpdir/Si.win") + @test isfile("$tmpdir/Si.wout") + @test !isfile("$tmpdir/Si.werr") + end + end + + @testset "Using Wannier.jl" begin + mktempdir() do tmpdir + wannier_model = Wannier.Model(scfres; + fileprefix="$tmpdir/Si", + n_wannier=8, bands_plot=true, + dis_win_max=17.185257, + dis_froz_max=6.8318033, + dis_mix_ratio=1.0, + wannier_plot=true) + + wannier_model.U .= disentangle(wannier_model; max_iter=500) + + # for now the Wannier.jl compat writes the amn, eig, mmn and win files + @test isfile("$tmpdir/Si.amn") + @test !isfile("$tmpdir/Si.chk") + @test isfile("$tmpdir/Si.eig") + @test isfile("$tmpdir/Si.mmn") + @test !isfile("$tmpdir/Si.nnkp") + @test isfile("$tmpdir/Si.win") + @test !isfile("$tmpdir/Si.wout") # NO .wout file, the output is in the `wannier_model` + @test !isfile("$tmpdir/Si.werr") + end + end end end - From f44c55f5ea4ca98dd17a3b6e9474af6070df86ba Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Mon, 11 Dec 2023 21:43:50 +0100 Subject: [PATCH 26/29] Remove opposite_gaussians_projection --- src/external/wannier_shared.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/external/wannier_shared.jl b/src/external/wannier_shared.jl index ff57a71b6e..7cdef6f7f0 100644 --- a/src/external/wannier_shared.jl +++ b/src/external/wannier_shared.jl @@ -21,15 +21,6 @@ function (proj::GaussianWannierProjection)(basis::PlaneWaveBasis, qs) end end -""" -A p-like initial guess, using the difference of two Gaussians with different centers. -""" -function opposite_gaussians_projection(center1::AbstractVector, center2::AbstractVector) - function inner(basis, qs) - GaussianWannierProjection(center1)(basis, qs) - GaussianWannierProjection(center2)(basis, qs) - end -end - @doc raw""" A hydrogenic initial guess. From 50bbf4025426ad02ba2d290ec345438c6176a0df Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Tue, 12 Dec 2023 12:00:36 +0100 Subject: [PATCH 27/29] Add custom projection example --- examples/wannier.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/examples/wannier.jl b/examples/wannier.jl index 64ea85b7b6..09371aef1a 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -103,6 +103,25 @@ interp_model = Wannier.InterpModel(wannier_model; kpath=kpath) # (Delete temporary files when done.) rm("wannier", recursive=true) +# ### Custom initial guesses +# +# We can also provide custom initial guesses for Wannierization, +# by passing a callable function in the `projections` array. +# The function receives the basis and a list of points (fractional coordinates in reciprocal space), +# and returns the Fourier transform of the initial guess function evaluated at each point. +# +# For example, we could use Gaussians for the σ and pz guesses with the following code: +s_guess(center) = DFTK.GaussianWannierProjection(center) +pz_guess(center) = begin + ## Approximate with two Gaussians offset by 0.5 Å from the center of the atom + offset = model.inv_lattice * [0, 0, austrip(0.5u"Å")] + center1 = center + offset + center2 = center - offset + ## Build the custom projector + (basis, qs) -> DFTK.GaussianWannierProjection(center1)(basis, qs) - DFTK.GaussianWannierProjection(center2)(basis, qs) +end +## Feed to Wannier via the `projections` as before... + # This example assumes that Wannier.jl version 0.3.2 is used, # and may need to be updated to accomodate for changes in Wannier.jl. # From b8f5e35965bd27095758345b9bca161df111cab1 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Tue, 12 Dec 2023 12:09:19 +0100 Subject: [PATCH 28/29] Remove wannier90 example that got added back because of the merge --- examples/wannier90.jl | 73 ------------------------------------------- 1 file changed, 73 deletions(-) delete mode 100644 examples/wannier90.jl diff --git a/examples/wannier90.jl b/examples/wannier90.jl deleted file mode 100644 index 5afb65837c..0000000000 --- a/examples/wannier90.jl +++ /dev/null @@ -1,73 +0,0 @@ -# # Wannierization using Wannier90 -# -# DFTK features an interface with the program -# [Wannier90](http://www.wannier.org/), -# in order to compute maximally-localized Wannier functions (MLWFs) -# from an initial self consistent field calculation. -# All processes are handled by calling the routine `run_wannier90`. -# -# !!! warning "No guarantees on Wannier90 interface" -# This code is at an early stage and has so far not been fully tested. -# Bugs are likely and we welcome issues in case you find any! -# -# This example shows how to obtain the MLWFs corresponding -# to the first five bands of graphene. Since the bands 2 to 11 are entangled, -# 15 bands are first computed to obtain 5 MLWFs by a disantanglement procedure. - -using DFTK -using Unitful -using UnitfulAtomic - -d = 10u"Å" -a = 2.641u"Å" # Graphene Lattice constant -lattice = [a -a/2 0; - 0 √3*a/2 0; - 0 0 d] - -C = ElementPsp(:C; psp=load_psp("hgh/pbe/c-q4")) -atoms = [C, C] -positions = [[0.0, 0.0, 0.0], [1//3, 2//3, 0.0]] -model = model_PBE(lattice, atoms, positions) -basis = PlaneWaveBasis(model; Ecut=15, kgrid=[5, 5, 1]) -nbandsalg = AdaptiveBands(basis.model; n_bands_converge=15) -scfres = self_consistent_field(basis; nbandsalg, tol=1e-5); - -# Plot bandstructure of the system - -bands = compute_bands(scfres; kline_density=10) -plot_bandstructure(bands) - -# Now we use the `run_wannier90` routine to generate all files needed by -# wannier90 and to perform the wannierization procedure. -# In Wannier90's convention, all files are named with the same prefix and only differ by -# their extensions. By default all generated input and output files are stored -# in the subfolder "wannier90" under the prefix "wannier" (i.e. "wannier90/wannier.win", -# "wannier90/wannier.wout", etc.). A different file prefix can be given with the -# keyword argument `fileprefix` as shown below. -# -# We now solve for 5 MLWF using wannier90: - -using wannier90_jll # Needed to make run_wannier90 available -run_wannier90(scfres; - fileprefix="wannier/graphene", - n_wannier=5, - num_print_cycles=25, - num_iter=200, - ## - dis_win_max=19.0, - dis_froz_max=0.1, - dis_num_iter=300, - dis_mix_ratio=1.0, - ## - wannier_plot=true, - wannier_plot_format="cube", - wannier_plot_supercell=5, - write_xyz=true, - translate_home_cell=true, - ); - -# As can be observed standard optional arguments for the disentanglement -# can be passed directly to `run_wannier90` as keyword arguments. - -# (Delete temporary files.) -rm("wannier", recursive=true) From 7bf3a2969de90700c5f289891eb292665993b765 Mon Sep 17 00:00:00 2001 From: Technici4n <13494793+Technici4n@users.noreply.github.com> Date: Tue, 12 Dec 2023 12:39:15 +0100 Subject: [PATCH 29/29] don't use begin with inline function definition --- examples/wannier.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/wannier.jl b/examples/wannier.jl index 73a8790a0d..9bb0589cee 100644 --- a/examples/wannier.jl +++ b/examples/wannier.jl @@ -113,7 +113,7 @@ rm("wannier", recursive=true) # # For example, we could use Gaussians for the σ and pz guesses with the following code: s_guess(center) = DFTK.GaussianWannierProjection(center) -pz_guess(center) = begin +function pz_guess(center) ## Approximate with two Gaussians offset by 0.5 Å from the center of the atom offset = model.inv_lattice * [0, 0, austrip(0.5u"Å")] center1 = center + offset