This structure provides an interface to the .wfc files.

The format is documented in: format of wfc files.

The relevant part is copied here for convenience:

• INTEGER :: ik k-point index (1 to number of k-points)
• REAL(8) :: xk(3) k-point coordinates
• INTEGER :: ispin spin index for LSDA case: ispin=1 for spin-up, ispin=2 for spin-down for unpolarized or non-colinear cases, ispin=1 always
• LOGICAL :: gamma_only if .true. write or read only half of the plane waves
• REAL(8) :: scalef scale factor applied to wavefunctions
• INTEGER :: ngw number of plane waves (PW)
• INTEGER :: igwx max number of PW (may be larger than ngw, not sure why)
• INTEGER :: npol number of spin states for PWs: 2 for non-colinear case, 1 otherwise
• INTEGER :: nbnd number of wavefunctions
• REAL(8) :: b1(3), b2(3), b3(3) primitive reciprocal lattice vectors
• INTEGER :: mill(3,igwx) miller indices: h=mill(1,i), k=mill(2,i), l=mill(3,i) the i-th PW has wave vector (k+G)(:)=xk(:)+h*b1(:)+k*b2(:)+ l*b3(:)
• COMPLEX(8) :: evc(npol*igwx,nbnd) wave functions in the PW basis set The first index runs on PW components, the second index runs on band states. For non-colinear case, each PW has a spin component first igwx components have PW with up spin, second igwx components have PW with down spin

Some variables have been renamed in the package to exorcise the acronyms.

• ik -> i_kpoint
• xk -> k_coordinates
• ispin -> i_spin
• scalef -> scale_factor
• ngw -> n_planewaves
• igwx -> max_n_planewaves
• npol -> n_polerizations
• nbnd -> n_band
• b1, b2, b3 -> b_1, b_2, b_3

WFC(wave_function_filename)

Load the metadata of a wave function from a wfc.dat (terrible filenames) file. This does not load the wave functions.

This is not guaranteed to work when Quantum Espresso is compiled with a Fortran compiler that is not gfort. Fortran does not enforce a standard on the size of it's primitive types in binary files.

Example:

Assume that you are in docs directory (change the path otherwise).

julia> path_to_si = "../test/test_data/test_5";
julia> wave_function = WFC(joinpath(path_to_si, "si.save/wfc1.dat"));
julia> wave_function.n_band, wave_function.gamma_only, wave_function.n_planewaves
(4, false, 537)

wave_functions_from_directory(save_dir)

Read all the wave function files in an directory for metadata. This does not load the wave functions.

The wfcX.dat files only gives you the content of a single wave function, but what you really want is probably reading the entire .save directory and convert them to orbitals on grids.

julia> wave_functions_list = wave_functions_from_directory(joinpath(path_to_si, "si.save"));
julia> length(wave_functions_list)
64

orbital_set_from_save(
    wave_functions_list,
    [domain_scaling_factor=2,]
)

Load the wavefunctions and parse them into orbitals.

The domain_scaling_factor is possibly either from Nyquist-Shannon sampling theorem. or from the energy cutoff for the density. Settting this to 2 gives the same FFT as QE.

julia> u = orbital_set_from_save(wave_functions_list);
julia> length(elements(u))
64

singleorbitalfromwavefunctions(w, l, k, n)

Use this if you want to load a single orbital. Mostly likely you will want orbital_set_from_save to load them all at once.

Create an orbital on the nth band from a WFC object w. The orbital will be defined on a reciprocal lattice l.

k is the kpoint, which can be out of the first Brillouin zone. If k is outside the first Brillouin zone, we will fold the kpoint into the first Brillouin zone and translate the wave function in on the reciprocal lattice (phase shift in real lattice).

The orbital is represented as values on the reciprocal lattice.

load_evc!(wfc, [bands=1:wfc.n_band])

Load the evc (frequency space wave-functions) into wave_functions. You probably shouldn't have to touch this.

i_kpoint_map(wave_function_list)

Construct a mapping between the mystical i_kpoints and kpoints from a list of WFC objects. Also figure out the Brillouin zone. For technical reasons, these two things are difficult to do separately.

This mapping may be different from what is used for .mmn and .amn files.

julia> k_map, brillouin_zone = i_kpoint_map(wave_functions_list);
julia> k_map[1]
GridVector{BrillouinZone3D}:
    coefficients: [0, 0, 0]

The UNKXXXXX.1 files are produced by pw2wannier90.x to show orbitals in real space.

This structure provides an interface to the UNK files.

• nx, ny, nz: Number of grid points in x, y, z direction.
• k: The kpoint at which the unk orbital is defined.
• n_band: The number of bands stored.
• psi_r: Each row of this matrix corresponds to a band. Each column corresponds to a grid point (not sure which one yet).
• filename: The name of the file from which the UNK object is constructed.

UNK(unk_filename)

Load a unk file named unk_filename as a UNK object.

The wave functions are read such that each row is an orbital to avoid transposing in memory for SCDM.

This takes 10s to load a 1.4G unk file using 12 threads and an NVMe drive. Since the unk files are too clumsy to cinlude in VCS, we will not provide the files.

Example:

wave_function = UNK(joinpath(path_to_si, "unk/UNK00001.1"))

wannier_from_unk_dir(
    unk_dir,
    wave_functions_list,
    [domain_scaling_factor=2,]
)

Load the wavefunctions and parse them into orbitals. Since the UNKXXXXX.1 does not carry much information, we need the wave_functions_list from the .save directory as well.

Example:

u_real = wannier_from_unk_dir(joinpath(path_to_si, "unk"), wave_functions_list)

single_orbital_from_unk(unk, h, k, n)

Use this if you want to load a single orbital. Mostly likely you will want wannier_from_unk_dir to load them all at once.

Create an orbital on the nth band from a UNK structure unk. The orbital will be defined on a homecell h.

k is the kpoint, which can be out of the first Brillouin zone.

The kpoint is not folded and the wave function is not phase shifted. Currently, this is done at a later stage after the fft. TODO: Phase shift the orbital.

The orbital is represented as values on the homecell.

This structure provides an interface to the MMN files.

• The first line should tell a joke. The parser won't parse unless tickled.

The format of the MMN files can be found in the user guide of wannier90.x.

Γ-point:

wannier90.x uses MMN files also for the X, Y, Z matrices of gamma point calculations. This is nowhere documented.

The $g$-vectors (taken from wannier90 user guide):

The last three integers specify the $G$ vector, in reciprocal lattice units, that brings the k-point specified by the second integer, and that thus lives inside the first Brillouin zone, to the actual $k + b$ we need.

MMN(mmn_filename, [ad_hoc_gamma = false])

Load a MMN file. Set ad_hoc_gamma = true if it is a gamma-only calculation so that the gamma-trick shenanigans are dealt with.

julia> mmn = MMN(joinpath(path_to_si, "output/pw2wan/si.mmn"))
julia> mmn.n_kpoint, mmn.n_neighbor
(64, 8)

The .amn files are used by Wannier90 for storing the gauge transform.

AMN(amn_filename)

Create an AMN object from a file.

julia> amn = AMN(joinpath(path_to_si, "output/pw2wan/si.amn"));
julia> amn.n_kpoint
64

Gauge(grid, amn, k_map, [orthonormalization = true])

Create a gauge from an AMN object and a k_map.

An AMN object also just stores the raw data, and we need to convert this to a Gauge to use it. Like .mmn files, we need a mapping of k-points as well as the Brillouin zone to construct the gauge, which can then be indexed with k-points.

julia> U = Gauge(brillouin_zone, amn, k_map);
julia> U[brillouin_zone[0, 0, 0]][:, :]
4×4 Matrix{ComplexF64}:
  -0.424021+0.0687693im  -0.530967+0.0861146im  -0.469845+0.0762009im  -0.540275+0.0876234im
 0.00775601-0.282049im     -0.3875-0.232102im    0.714604+0.092279im   -0.246713+0.369214im
   0.837758+0.116721im   -0.437746+0.134205im   -0.142503-0.20453im    -0.103362-0.045631im
  0.0694477+0.124822im    0.516605+0.173543im    0.152318-0.411003im   -0.694672+0.0889086im