mirror of
https://github.com/triqs/dft_tools
synced 2024-12-23 04:43:42 +01:00
Revision of Wannier90Converter:
1) write the Hamiltonian in Bloch space for charge self-consistent calculations, and 2) spin-orbit coupling if bloch_basis = True: * if "seedname_u.dat" (and in case of disentanglement "seedname_u_dis.dat") present, write hopping in Bloch basis * "proj_mat" transforming from Bloch to orbital space * diagonal hoppings are directly read from "seedname.eig" * fermi weights and band_window of Wannier Hamiltonian are read from DFT output and "seedname.nnkp", written into new subgroup "dft_misc_input" * automatic calculation of "density_required" * implemented for Quantum Espresso (read from "seedname.nscf.out" if verbosity = 'high') and VASP (read from "OUTCAR"/"LOCPROJ") * spin-orbit coupling SO = 1 implemented * substitute k_mesh and bz_weights with kpts and kpt_weights, respectively (previous names kept for compatibility) * updated tests
This commit is contained in:
commit
92e24a9ada
@ -24,14 +24,15 @@
|
||||
# Wannier90 to HDF5 converter for the SumkDFT class of dfttools/TRIQS;
|
||||
#
|
||||
# written by Gabriele Sclauzero (Materials Theory, ETH Zurich), Dec 2015 -- Jan 2016,
|
||||
# and updated by Maximilian Merkel (Materials Theory, ETH Zurich), Aug 2020,
|
||||
# updated by Maximilian Merkel (Materials Theory, ETH Zurich), Aug 2020 -- Jan 2021,
|
||||
# and by Sophie Beck (Materials Theory, ETH Zurich), Sep 2020 -- Jan 2021,
|
||||
# under the supervision of Claude Ederer (Materials Theory).
|
||||
# Partially based on previous work by K. Dymkovski and the DFT_tools/TRIQS team.
|
||||
#
|
||||
# Limitations of the current implementation:
|
||||
# - the case with SO=1 is not considered at the moment
|
||||
# - the T rotation matrices are not used in this implementation
|
||||
# - projectors for uncorrelated shells (proj_mat_all) cannot be set
|
||||
# - in bloch_basis mode, the number of wannier functions per shell has to be equal
|
||||
# to the dim of the shell
|
||||
#
|
||||
# Things to be improved/checked:
|
||||
# - the case with SP=1 might work, but was never tested (do we need to define
|
||||
@ -46,11 +47,11 @@
|
||||
|
||||
|
||||
import numpy
|
||||
import math
|
||||
from h5 import *
|
||||
from .converter_tools import *
|
||||
from itertools import product
|
||||
import os.path
|
||||
from itertools import product
|
||||
|
||||
from h5 import HDFArchive
|
||||
from .converter_tools import ConverterTools
|
||||
import triqs.utility.mpi as mpi
|
||||
|
||||
class Wannier90Converter(ConverterTools):
|
||||
@ -59,8 +60,8 @@ class Wannier90Converter(ConverterTools):
|
||||
"""
|
||||
|
||||
def __init__(self, seedname, hdf_filename=None, dft_subgrp='dft_input',
|
||||
symmcorr_subgrp='dft_symmcorr_input', repacking=False,
|
||||
rot_mat_type='hloc_diag'):
|
||||
symmcorr_subgrp='dft_symmcorr_input', misc_subgrp='dft_misc_input',
|
||||
repacking=False, rot_mat_type='hloc_diag', bloch_basis=False):
|
||||
"""
|
||||
Initialise the class.
|
||||
|
||||
@ -74,11 +75,15 @@ class Wannier90Converter(ConverterTools):
|
||||
Name of subgroup storing necessary DFT data
|
||||
symmcorr_subgrp : string, optional
|
||||
Name of subgroup storing correlated-shell symmetry data
|
||||
misc_subgrp : string, optional
|
||||
Name of subgroup storing miscellaneous DFT data.
|
||||
repacking : boolean, optional
|
||||
Does the hdf5 archive need to be repacked to save space?
|
||||
rot_mat_type : string, optional
|
||||
Type of rot_mat used
|
||||
Can be 'hloc_diag', 'wannier', 'none'
|
||||
bloch_basis : boolean, optional
|
||||
Should the Hamiltonian be written in Bloch rather than Wannier basis?
|
||||
"""
|
||||
|
||||
self._name = "Wannier90Converter"
|
||||
@ -93,11 +98,13 @@ class Wannier90Converter(ConverterTools):
|
||||
self.w90_seed = seedname
|
||||
self.dft_subgrp = dft_subgrp
|
||||
self.symmcorr_subgrp = symmcorr_subgrp
|
||||
self.misc_subgrp = misc_subgrp
|
||||
self.fortran_to_replace = {'D': 'E'}
|
||||
# threshold below which matrix elements from wannier90 should be
|
||||
# considered equal
|
||||
self._w90zero = 2.e-6
|
||||
self.rot_mat_type = rot_mat_type
|
||||
self.bloch_basis = bloch_basis
|
||||
if self.rot_mat_type not in ('hloc_diag', 'wannier', 'none'):
|
||||
raise ValueError('Parameter rot_mat_type invalid, should be one of'
|
||||
+ '"hloc_diag", "wannier", "none"')
|
||||
@ -113,14 +120,14 @@ class Wannier90Converter(ConverterTools):
|
||||
- dft_subgrp
|
||||
- symmcorr_subgrp
|
||||
|
||||
in the hdf5 archive.
|
||||
in the hdf5 archive.
|
||||
|
||||
"""
|
||||
|
||||
# Read and write only on the master node
|
||||
if not (mpi.is_master_node()):
|
||||
return
|
||||
mpi.report("Reading input from %s..." % self.inp_file)
|
||||
mpi.report("\nReading input from %s..." % self.inp_file)
|
||||
|
||||
# R is a generator : each R.Next() will return the next number in the
|
||||
# file
|
||||
@ -139,8 +146,10 @@ class Wannier90Converter(ConverterTools):
|
||||
else:
|
||||
# some default grid, if everything else fails...
|
||||
nki = [8, 8, 8]
|
||||
# read the total number of electrons per cell
|
||||
# read the total number of electrons per cell if not in bloch basis
|
||||
# in bloch basis, this is later calculated from the partial occupations
|
||||
density_required = float(next(R))
|
||||
|
||||
# we do not read shells, because we have no additional shells beyond correlated ones,
|
||||
# and the data will be copied from corr_shells into shells (see below)
|
||||
# number of corr. shells (e.g. Fe d, Ce f) in the unit cell,
|
||||
@ -149,6 +158,7 @@ class Wannier90Converter(ConverterTools):
|
||||
# l, dim, SO flag, irep):
|
||||
corr_shells = [{name: int(val) for name, val in zip(
|
||||
corr_shell_entries, R)} for icrsh in range(n_corr_shells)]
|
||||
|
||||
try:
|
||||
self.fermi_energy = float(next(R))
|
||||
except:
|
||||
@ -169,17 +179,40 @@ class Wannier90Converter(ConverterTools):
|
||||
for ish in range(n_shells):
|
||||
shells.append({key: corr_shells[ish].get(
|
||||
key, None) for key in shell_entries})
|
||||
###
|
||||
SP = 0 # NO spin-polarised calculations for now
|
||||
SO = 0 # NO spin-orbit calculation for now
|
||||
|
||||
# Determine if any shell requires SO
|
||||
if any([corr_shell['SO']==1 for corr_shell in corr_shells]):
|
||||
SO = 1
|
||||
SP = 1
|
||||
print('Spin-orbit interaction turned on')
|
||||
else:
|
||||
SO = 0
|
||||
SP = 0
|
||||
# Only one block supported - either non-spin-polarized or spin-orbit coupled
|
||||
assert SP == SO, 'Spin-polarized calculations not implemented'
|
||||
|
||||
charge_below = 0 # total charge below energy window NOT used for now
|
||||
energy_unit = 1.0 # should be understood as eV units
|
||||
###
|
||||
|
||||
# this is more general
|
||||
n_spin = SP + 1 - SO
|
||||
n_spin_blocs = SP + 1 - SO
|
||||
assert n_spin_blocs > 0, 'Input error, if SO=1, SP must be 1.'
|
||||
|
||||
if SO == 1:
|
||||
for shell_list in [shells, corr_shells]:
|
||||
for entry in shell_list:
|
||||
entry['dim'] *= 2
|
||||
|
||||
dim_corr_shells = sum([sh['dim'] for sh in corr_shells])
|
||||
mpi.report(
|
||||
"Total number of WFs expected in the correlated shells: %d" % dim_corr_shells)
|
||||
mpi.report('Total number of WFs expected in the correlated shells: {0:d}'.format(dim_corr_shells))
|
||||
|
||||
# build the k-point mesh, if its size was given on input (kmesh_mode >= 0),
|
||||
# otherwise it is built according to the data in the hr file (see below)
|
||||
# If output is in bloch_basis, we use k mesh from seedname_u.mat for consistency
|
||||
if kmesh_mode >= 0 and not self.bloch_basis:
|
||||
n_k, kpts, kpt_weights = self.kmesh_build(nki, kmesh_mode)
|
||||
self.n_k = n_k
|
||||
self.kpts = kpts
|
||||
|
||||
# determine the number of inequivalent correlated shells and maps,
|
||||
# needed for further processing
|
||||
@ -190,15 +223,6 @@ class Wannier90Converter(ConverterTools):
|
||||
shells_map = [inequiv_to_corr[corr_to_inequiv[ish]]
|
||||
for ish in range(n_corr_shells)]
|
||||
mpi.report("Mapping: " + format(shells_map))
|
||||
mpi.report("Subtracting %f eV from the Fermi level." % self.fermi_energy)
|
||||
|
||||
# build the k-point mesh, if its size was given on input (kmesh_mode >= 0),
|
||||
# otherwise it is built according to the data in the hr file (see
|
||||
# below)
|
||||
if kmesh_mode >= 0:
|
||||
n_k, k_mesh, bz_weights = self.kmesh_build(nki, kmesh_mode)
|
||||
self.n_k = n_k
|
||||
self.k_mesh = k_mesh
|
||||
|
||||
# not used in this version: reset to dummy values?
|
||||
n_reps = [1 for i in range(n_inequiv_shells)]
|
||||
@ -207,33 +231,37 @@ class Wannier90Converter(ConverterTools):
|
||||
for ish in range(n_inequiv_shells):
|
||||
ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1
|
||||
lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
|
||||
T.append(numpy.zeros([lmax, lmax], numpy.complex_))
|
||||
T.append(numpy.zeros([lmax, lmax], dtype=complex))
|
||||
|
||||
spin_w90name = ['_up', '_down']
|
||||
hamr_full = []
|
||||
umat_full = []
|
||||
udismat_full = []
|
||||
bandmat_full = []
|
||||
|
||||
# TODO: generalise to SP=1 (only partially done)
|
||||
rot_mat_time_inv = [0 for i in range(n_corr_shells)]
|
||||
|
||||
# Second, let's read the file containing the Hamiltonian in WF basis
|
||||
# produced by Wannier90
|
||||
for isp in range(n_spin):
|
||||
for isp in range(n_spin_blocs):
|
||||
# begin loop on isp
|
||||
|
||||
# build filename according to wannier90 conventions
|
||||
if SP == 1:
|
||||
if SP == 1 and SO == 0:
|
||||
mpi.report(
|
||||
"Reading information for spin component n. %d" % isp)
|
||||
hr_file = self.w90_seed + spin_w90name[isp] + '_hr.dat'
|
||||
file_seed = self.w90_seed + spin_w90name[isp]
|
||||
else:
|
||||
hr_file = self.w90_seed + '_hr.dat'
|
||||
file_seed = self.w90_seed
|
||||
# now grab the data from the H(R) file
|
||||
mpi.report(
|
||||
"The Hamiltonian in MLWF basis is extracted from %s ..." % hr_file)
|
||||
nr, rvec, rdeg, nw, hamr = self.read_wannier90hr(hr_file)
|
||||
"\nThe Hamiltonian in MLWF basis is extracted from %s files..." % file_seed)
|
||||
(nr, rvec, rdeg, nw, hamr, u_mat, udis_mat,
|
||||
band_mat, k_mesh_from_umat) = self.read_wannier90data(file_seed, kmesh_mode)
|
||||
# number of R vectors, their indices, their degeneracy, number of
|
||||
# WFs, H(R)
|
||||
mpi.report("... done: %d R vectors, %d WFs found" % (nr, nw))
|
||||
mpi.report("\n... done: %d R vectors, %d WFs found" % (nr, nw))
|
||||
|
||||
if isp == 0:
|
||||
# set or check some quantities that must be the same for both
|
||||
@ -241,7 +269,14 @@ class Wannier90Converter(ConverterTools):
|
||||
self.nrpt = nr
|
||||
|
||||
# k-point grid: (if not defined before)
|
||||
if kmesh_mode == -1:
|
||||
if self.bloch_basis:
|
||||
mpi.report('Reading k mesh from seedname_u.mat file')
|
||||
kpts = k_mesh_from_umat
|
||||
n_k = len(kpts)
|
||||
kpt_weights = numpy.full(n_k, 1/n_k)
|
||||
self.n_k = n_k
|
||||
self.kpts = kpts
|
||||
elif kmesh_mode == -1:
|
||||
# the size of the k-point mesh is determined from the
|
||||
# largest R vector
|
||||
nki = [2 * rvec[:, idir].max() + 1 for idir in range(3)]
|
||||
@ -249,9 +284,10 @@ class Wannier90Converter(ConverterTools):
|
||||
# wannier90 convention: if we have nki k-points along the i-th direction,
|
||||
# then we should get 2*(nki/2)+nki%2 R points along that
|
||||
# direction
|
||||
n_k, k_mesh, bz_weights = self.kmesh_build(nki)
|
||||
self.n_k = n_k
|
||||
self.k_mesh = k_mesh
|
||||
n_k, kpts, kpt_weights = self.kmesh_build(nki)
|
||||
self.n_k = n_k
|
||||
self.kpts = kpts
|
||||
|
||||
|
||||
# set the R vectors and their degeneracy
|
||||
self.rvec = rvec
|
||||
@ -273,9 +309,10 @@ class Wannier90Converter(ConverterTools):
|
||||
|
||||
# we assume spin up and spin down always have same total number
|
||||
# of WFs
|
||||
n_orbitals = numpy.ones(
|
||||
[self.n_k, n_spin], numpy.int) * self.nwfs
|
||||
|
||||
# get second dimension of udis_mat which corresponds to number of bands in window
|
||||
# n_bands_max corresponds to numpy.max(n_orbitals)
|
||||
n_bands_max = udis_mat.shape[1]
|
||||
n_orbitals = numpy.full([self.n_k, n_spin_blocs], n_bands_max)
|
||||
else:
|
||||
# consistency check between the _up and _down file contents
|
||||
if nr != self.nrpt:
|
||||
@ -286,8 +323,9 @@ class Wannier90Converter(ConverterTools):
|
||||
"Different number of WFs for spin-up/spin-down!")
|
||||
|
||||
hamr_full.append(hamr)
|
||||
# FIXME: when do we actually need deepcopy()?
|
||||
# hamr_full.append(deepcopy(hamr))
|
||||
umat_full.append(u_mat)
|
||||
udismat_full.append(udis_mat)
|
||||
bandmat_full.append(band_mat)
|
||||
|
||||
for ir in range(nr):
|
||||
# checks if the Hamiltonian is real (it should, if
|
||||
@ -322,39 +360,92 @@ class Wannier90Converter(ConverterTools):
|
||||
"Rotations for spin component n. %d do not match!" % isp)
|
||||
# end loop on isp
|
||||
|
||||
mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki))
|
||||
# Reads misc input needed for CSC calculations
|
||||
if self.bloch_basis:
|
||||
if os.path.isfile(self.w90_seed + '.nscf.out'):
|
||||
fermi_weight_file = self.w90_seed + '.nscf.out'
|
||||
print('Reading DFT band occupations from Quantum Espresso output {}'.format(fermi_weight_file))
|
||||
elif os.path.isfile('LOCPROJ'):
|
||||
fermi_weight_file = 'LOCPROJ'
|
||||
print('Reading DFT band occupations from Vasp output {}'.format(fermi_weight_file))
|
||||
else:
|
||||
raise IOError('seedname.nscf.out or LOCPROJ required in bloch_basis mode')
|
||||
|
||||
f_weights, band_window, self.fermi_energy = self.convert_misc_input(fermi_weight_file,
|
||||
self.w90_seed + '.nnkp', n_spin_blocs)
|
||||
# Get density from k-point weighted average and sum over all spins and bands
|
||||
density_required = numpy.sum(f_weights.T * kpt_weights) * (2 - SP)
|
||||
print('Overwriting required density with DFT result {:.5f}'.format(density_required))
|
||||
print('and using the DFT Fermi energy {:.5f} eV\n'.format(self.fermi_energy))
|
||||
|
||||
if not self.bloch_basis:
|
||||
mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki))
|
||||
|
||||
# if calculations are spin-polarized, then renormalize k-point weights
|
||||
if SP == 1:
|
||||
bz_weights = 0.5 * bz_weights
|
||||
if SP == 1 and SO == 0:
|
||||
kpt_weights *= 0.5
|
||||
|
||||
# Third, compute the hoppings in reciprocal space
|
||||
hopping = numpy.zeros([self.n_k, n_spin, numpy.max(
|
||||
n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
|
||||
for isp in range(n_spin):
|
||||
# make Fourier transform H(R) -> H(k) : it can be done one spin at
|
||||
# a time
|
||||
hamk = self.fourier_ham(self.nwfs, hamr_full[isp])
|
||||
# copy the H(k) in the right place of hoppings... is there a better
|
||||
# way to do this??
|
||||
for ik in range(self.n_k):
|
||||
#hopping[ik,isp,:,:] = deepcopy(hamk[ik][:,:])*energy_unit
|
||||
hopping[ik, isp, :, :] = hamk[ik][:, :] * energy_unit
|
||||
|
||||
# Then, initialise the projectors
|
||||
k_dep_projection = 0 # we always have the same number of WFs at each k-point
|
||||
proj_mat = numpy.zeros([self.n_k, n_spin, n_corr_shells, max(
|
||||
[crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
|
||||
# Third, initialise the projectors
|
||||
k_dep_projection = 0 # at the moment not really used, but might get important
|
||||
proj_mat = numpy.zeros([self.n_k, n_spin_blocs, n_corr_shells, max(
|
||||
[crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], dtype=complex)
|
||||
iorb = 0
|
||||
# Projectors simply consist in identity matrix blocks selecting those MLWFs that
|
||||
# correspond to the specific correlated shell indexed by icrsh.
|
||||
# Projectors are either identity matrix blocks to use with Wannier basis
|
||||
# OR correspond to the overlap between Kohn-Sham and Wannier orbitals as
|
||||
# P_{nu,alpha](k) = <w_{alpha,k}|psi_{nu,k}>
|
||||
# NOTE: we assume that the correlated orbitals appear at the beginning of the H(R)
|
||||
# file and that the ordering of MLWFs matches the corr_shell info from
|
||||
# the input.
|
||||
for icrsh in range(n_corr_shells):
|
||||
norb = corr_shells[icrsh]['dim']
|
||||
proj_mat[:, :, icrsh, 0:norb, iorb:iorb +
|
||||
norb] = numpy.identity(norb, numpy.complex_)
|
||||
iorb += norb
|
||||
for isp in range(n_spin_blocs):
|
||||
# now combine udismat and umat
|
||||
u_total = numpy.einsum('abc,acd->abd',udismat_full[isp],umat_full[isp])
|
||||
# transpose and write into proj_mat
|
||||
u_temp = numpy.transpose(u_total.conj(),(0,2,1))
|
||||
for icrsh in range(n_corr_shells):
|
||||
dim = corr_shells[icrsh]['dim']
|
||||
proj_mat[:, isp, icrsh, 0:dim, :] = u_temp[:,iorb:iorb+dim,:]
|
||||
iorb += dim
|
||||
|
||||
# Then, compute the hoppings in reciprocal space
|
||||
hopping = numpy.zeros([self.n_k, n_spin_blocs, numpy.max(n_orbitals), numpy.max(n_orbitals)], dtype=complex)
|
||||
for isp in range(n_spin_blocs):
|
||||
# if bloch_basis is True, use Kohn-Sham eigenvalues as hamk
|
||||
# this ensures that the calculation of the band-correlation energy
|
||||
# is consistent with SumkDFT's calc_density_correction
|
||||
if self.bloch_basis:
|
||||
# diagonal Kohn-Sham bands
|
||||
# TODO: test for system with self.nwfs > dim_corr_shells
|
||||
hamk = [numpy.diag(bandmat_full[isp][ik]) for ik in range(self.n_k)]
|
||||
|
||||
# Sanity check if the local Hamiltonian with the projector method
|
||||
# corresponds to W90 result
|
||||
wannier_ham = self.fourier_ham(hamr_full[isp])
|
||||
for ik in range(self.n_k):
|
||||
proj_mat_flattened = proj_mat[ik, isp].reshape(self.nwfs, numpy.max(n_orbitals))
|
||||
downfolded_ham = proj_mat_flattened.dot(hamk[ik].dot(proj_mat_flattened.conj().T))
|
||||
|
||||
if not numpy.allclose(downfolded_ham, wannier_ham[ik], atol=self._w90zero, rtol=0):
|
||||
mpi.report('WARNING: mismatch between downfolded Hamiltonian and '
|
||||
+ f'Fourier transformed H(R). First occurred at kpt {ik}:')
|
||||
|
||||
with numpy.printoptions(formatter={'complexfloat': '{:+.3f}'.format}):
|
||||
mpi.report('Downfolded Hamiltonian, P H_eig P')
|
||||
mpi.report(downfolded_ham)
|
||||
mpi.report('\nWannier Hamiltonian, Fourier(H(r))')
|
||||
mpi.report(wannier_ham[ik])
|
||||
break
|
||||
# else for an isolated set of bands use fourier transform of H(R)
|
||||
else:
|
||||
# make Fourier transform H(R) -> H(k) : it can be done one spin at a time
|
||||
hamk = self.fourier_ham(hamr_full[isp])
|
||||
# finally write hamk into hoppings
|
||||
for ik in range(self.n_k):
|
||||
hopping[ik, isp] = hamk[ik] - numpy.identity(numpy.max(n_orbitals)) * self.fermi_energy
|
||||
hopping *= energy_unit
|
||||
mpi.report("Subtracting {:.5f} eV from the Fermi level.".format(self.fermi_energy))
|
||||
|
||||
# bz_weights required by triqs h5 standard but soon to be replaced by kpt_weights
|
||||
bz_weights = kpt_weights
|
||||
|
||||
# Finally, save all required data into the HDF archive:
|
||||
# use_rotations is supposed to be an int = 0, 1, no bool
|
||||
@ -367,18 +458,25 @@ class Wannier90Converter(ConverterTools):
|
||||
things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required',
|
||||
'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat',
|
||||
'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping',
|
||||
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
|
||||
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr', 'kpt_weights', 'kpts']
|
||||
for it in things_to_save:
|
||||
ar[self.dft_subgrp][it] = locals()[it]
|
||||
|
||||
def read_wannier90hr(self, hr_filename="wannier_hr.dat"):
|
||||
if self.bloch_basis:
|
||||
# Store Fermi weights to 'dft_misc_input'
|
||||
if not (self.misc_subgrp in ar):
|
||||
ar.create_group(self.misc_subgrp)
|
||||
ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights
|
||||
ar[self.misc_subgrp]['band_window'] = band_window+1 # Change to 1-based index
|
||||
|
||||
def read_wannier90data(self, wannier_seed="wannier", kmesh_mode=0):
|
||||
"""
|
||||
Method for reading the seedname_hr.dat file produced by Wannier90 (http://wannier.org)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
hr_filename : string
|
||||
full name of the H(R) file produced by Wannier90 (usually seedname_hr.dat)
|
||||
wannier_seed : string
|
||||
seedname to read H(R) file produced by Wannier90 (usually seedname_hr.dat)
|
||||
|
||||
Returns
|
||||
-------
|
||||
@ -392,21 +490,30 @@ class Wannier90Converter(ConverterTools):
|
||||
number of Wannier functions found
|
||||
h_of_r : list of numpy.array
|
||||
<w_i|H(R)|w_j> = Hamilonian matrix elements in the Wannier basis
|
||||
|
||||
u_mat : numpy.array
|
||||
U_mn^k = unitary matrix elements which mix the Kohn-Sham states
|
||||
udis_mat : numpy.array
|
||||
U^dis(k) = rectangular matrix for entangled bands
|
||||
band_mat : numpy.array
|
||||
\epsilon_nk = Kohn-Sham eigenvalues (in eV) needed for entangled bands
|
||||
If not self.bloch_basis unused and therefore None
|
||||
k_mesh : numpy.array
|
||||
The k mesh read from the seedname_u.mat file to ensure consistency
|
||||
If not self.bloch_basis unused and therefore None
|
||||
"""
|
||||
|
||||
# Read only from the master node
|
||||
if not (mpi.is_master_node()):
|
||||
return
|
||||
|
||||
hr_filename = wannier_seed + '_hr.dat'
|
||||
try:
|
||||
with open(hr_filename, "r") as hr_filedesc:
|
||||
with open(hr_filename, 'r') as hr_filedesc:
|
||||
hr_data = hr_filedesc.readlines()
|
||||
hr_filedesc.close()
|
||||
except IOError:
|
||||
mpi.report("The file %s could not be read!" % hr_filename)
|
||||
|
||||
mpi.report("Reading %s..." % hr_filename + hr_data[0])
|
||||
mpi.report('reading {:20}...{}'.format(hr_filename,hr_data[0].strip('\n')))
|
||||
|
||||
try:
|
||||
# reads number of Wannier functions per spin
|
||||
@ -415,11 +522,79 @@ class Wannier90Converter(ConverterTools):
|
||||
except ValueError:
|
||||
mpi.report("Could not read number of WFs or R vectors")
|
||||
|
||||
k_mesh = None
|
||||
|
||||
if not self.bloch_basis:
|
||||
# For kmesh_mode == -1, size of automatic k mesh known when rvec_idx
|
||||
# have been read. For kmesh_mode >= 0, kpts have been determined already
|
||||
if kmesh_mode >= 0:
|
||||
n_k = self.n_k
|
||||
else:
|
||||
# first, read u matrices from 'seedname_u.mat'
|
||||
u_filename = wannier_seed + '_u.mat'
|
||||
with open(u_filename,'r') as u_file:
|
||||
u_data = u_file.readlines()
|
||||
# reads number of kpoints and number of wannier functions
|
||||
n_k, num_wf_u, _ = map(int, u_data[1].split())
|
||||
assert num_wf_u == num_wf, '#WFs must be identical for *_u.mat and *_hr.dat'
|
||||
mpi.report('reading {:20}...{}'.format(u_filename,u_data[0].strip('\n')))
|
||||
|
||||
# Reads k mesh from all lines with 3 floats
|
||||
k_mesh = numpy.loadtxt((line for line in u_data if line.count('.') == 3))
|
||||
assert k_mesh.shape == (n_k, 3)
|
||||
# Reads u matrices from all lines with 2 floats
|
||||
u_mat = numpy.loadtxt((line for line in u_data if line.count('.') == 2))
|
||||
assert u_mat.shape == (n_k*num_wf*num_wf, 2)
|
||||
|
||||
mpi.report('Writing h5 archive in projector formalism: H(k) defined in KS Bloch basis')
|
||||
|
||||
try:
|
||||
# read 'seedname_u_dis.mat'
|
||||
udis_filename = wannier_seed + '_u_dis.mat'
|
||||
# if it exists the Kohn-Sham eigenvalues and the window are needed
|
||||
band_filename = wannier_seed + '.eig'
|
||||
wout_filename = wannier_seed + '.wout'
|
||||
|
||||
with open(udis_filename,'r') as udis_file:
|
||||
udis_data = udis_file.readlines()
|
||||
disentangle = True
|
||||
except IOError:
|
||||
disentangle = False
|
||||
mpi.report('WARNING: File {} missing.'.format(udis_filename))
|
||||
mpi.report('Assuming an isolated set of bands. Check if this is what you want!')
|
||||
|
||||
# read Kohn-Sham eigenvalues from 'seedname.eig'
|
||||
mpi.report('Reading {}'.format(band_filename))
|
||||
band_data = numpy.loadtxt(band_filename, usecols=2)
|
||||
|
||||
if disentangle:
|
||||
# reads number of kpoints, number of wannier functions and bands
|
||||
num_k_udis, num_wf_udis, num_ks_bands = map(int, udis_data[1].split())
|
||||
|
||||
assert num_k_udis == n_k, '#k points must be identical for *.inp and *_u_dis.mat'
|
||||
assert num_wf_udis == num_wf, '#WFs must be identical for *_u.mat and *_hr.dat'
|
||||
|
||||
mpi.report('Found {:22}...{}, '.format(udis_filename,udis_data[0].strip('\n')))
|
||||
udis_data = numpy.loadtxt(udis_data, usecols=(0, 1), skiprows=2)
|
||||
|
||||
# read disentanglement window from 'seedname.wout'
|
||||
with open(wout_filename) as wout_file:
|
||||
for line in wout_file:
|
||||
if 'Outer:' in line:
|
||||
content = line.split()
|
||||
index = content.index('Outer:') + 1
|
||||
dis_window_min = float(content[index])
|
||||
dis_window_max = float(content[index+2])
|
||||
break
|
||||
mpi.report('and {} for disentanglement energy window.'.format(wout_filename))
|
||||
else:
|
||||
num_ks_bands = num_wf
|
||||
|
||||
# allocate arrays to save the R vector indexes and degeneracies and the
|
||||
# Hamiltonian
|
||||
rvec_idx = numpy.zeros((nrpt, 3), dtype=int)
|
||||
rvec_deg = numpy.zeros(nrpt, dtype=int)
|
||||
h_of_r = [numpy.zeros((num_wf, num_wf), dtype=numpy.complex_)
|
||||
h_of_r = [numpy.zeros((num_wf, num_wf), dtype=complex)
|
||||
for n in range(nrpt)]
|
||||
|
||||
# variable currpos points to the current line in the file
|
||||
@ -437,7 +612,7 @@ class Wannier90Converter(ConverterTools):
|
||||
ir += 1
|
||||
# for each direct lattice vector R read the block of the
|
||||
# Hamiltonian H(R)
|
||||
for ir, jj, ii in product(list(range(nrpt)), list(range(num_wf)), list(range(num_wf))):
|
||||
for ir, jj, ii in product(range(nrpt), range(num_wf), range(num_wf)):
|
||||
# advance one line, split the line into tokens
|
||||
currpos += 1
|
||||
cline = hr_data[currpos].split()
|
||||
@ -445,8 +620,7 @@ class Wannier90Converter(ConverterTools):
|
||||
if int(cline[3]) != ii + 1 or int(cline[4]) != jj + 1:
|
||||
mpi.report(
|
||||
"Inconsistent indices at %s%s of R n. %s" % (ii, jj, ir))
|
||||
rcurr = numpy.array(
|
||||
[int(cline[0]), int(cline[1]), int(cline[2])])
|
||||
rcurr = numpy.array([int(cline[0]), int(cline[1]), int(cline[2])])
|
||||
if ii == 0 and jj == 0:
|
||||
rvec_idx[ir] = rcurr
|
||||
rprec = rcurr
|
||||
@ -457,16 +631,57 @@ class Wannier90Converter(ConverterTools):
|
||||
"Inconsistent indices for R vector n. %s" % ir)
|
||||
|
||||
# fill h_of_r with the matrix elements of the Hamiltonian
|
||||
if not numpy.any(rcurr) and ii == jj:
|
||||
h_of_r[ir][ii, jj] = complex(float(cline[5]) - self.fermi_energy, float(cline[6]))
|
||||
else:
|
||||
h_of_r[ir][ii, jj] = complex(float(cline[5]), float(cline[6]))
|
||||
h_of_r[ir][ii, jj] = complex(float(cline[5]), float(cline[6]))
|
||||
|
||||
except ValueError:
|
||||
mpi.report("Wrong data or structure in file %s" % hr_filename)
|
||||
|
||||
# first, get the input for u_mat
|
||||
if self.bloch_basis:
|
||||
u_mat = u_mat[:, 0] + 1j * u_mat[:, 1]
|
||||
u_mat = u_mat.reshape((n_k, num_wf, num_wf)).transpose((0, 2, 1))
|
||||
else:
|
||||
if kmesh_mode == -1:
|
||||
n_k = numpy.prod([2 * rvec_idx[:, idir].max() + 1 for idir in range(3)])
|
||||
# Wannier basis; fill u_mat with identity
|
||||
u_mat = numpy.zeros([n_k, num_wf, num_wf], dtype=complex)
|
||||
for ik in range(n_k):
|
||||
u_mat[ik,:,:] = numpy.identity(num_wf,dtype=complex)
|
||||
|
||||
# now, check what is needed in the case of disentanglement:
|
||||
# The file seedname_u_dis.mat contains only the bands inside the window
|
||||
# and then fills the rest up with zeros. Therefore, we need to put the
|
||||
# entries from udis_data in the correct position in udis_mat, i.e.
|
||||
# shifting by the number of bands below dis_window_min
|
||||
if self.bloch_basis:
|
||||
# reshape band_data
|
||||
band_mat = band_data.reshape(n_k, num_ks_bands)
|
||||
else:
|
||||
# Not used in wannier basis
|
||||
band_mat = None
|
||||
|
||||
if self.bloch_basis and disentangle:
|
||||
# Determine which bands are inside the band window
|
||||
inside_window = numpy.logical_and(band_mat >= dis_window_min,
|
||||
band_mat <= dis_window_max)
|
||||
n_inside_per_k = numpy.sum(inside_window, axis=1)
|
||||
|
||||
# Reformats udis_data as complex, without header
|
||||
udis_data = udis_data[:, 0] + 1j * udis_data[:, 1]
|
||||
udis_data = udis_data.reshape((n_k, num_wf*num_ks_bands+1))[:, 1:]
|
||||
udis_data = udis_data.reshape((n_k, num_wf, num_ks_bands))
|
||||
|
||||
#initiate U disentanglement matrices and fill from file "seedname_u_dis.mat"
|
||||
udis_mat = numpy.zeros([n_k, num_ks_bands, num_wf], dtype=complex)
|
||||
for ik in range(n_k):
|
||||
udis_mat[ik, inside_window[ik]] = udis_data[ik, :, :n_inside_per_k[ik]].T
|
||||
assert numpy.allclose(udis_data[ik, :, n_inside_per_k[ik]:], 0)
|
||||
else:
|
||||
# no disentanglement; fill udis_mat with identity
|
||||
udis_mat = numpy.array([numpy.identity(num_wf,dtype=complex)] * n_k)
|
||||
|
||||
# return the data into variables
|
||||
return nrpt, rvec_idx, rvec_deg, num_wf, h_of_r
|
||||
return nrpt, rvec_idx, rvec_deg, num_wf, h_of_r, u_mat, udis_mat, band_mat, k_mesh
|
||||
|
||||
def find_rot_mat(self, n_sh, sh_lst, sh_map, ham0):
|
||||
"""
|
||||
@ -588,7 +803,7 @@ class Wannier90Converter(ConverterTools):
|
||||
|
||||
return succeeded, rot_mat
|
||||
|
||||
def kmesh_build(self, msize=None, mmode=0):
|
||||
def kmesh_build(self, msize, mmode=0):
|
||||
"""
|
||||
Method for the generation of the k-point mesh.
|
||||
Right now it only supports the option for generating a full grid containing k=0,0,0.
|
||||
@ -604,7 +819,7 @@ class Wannier90Converter(ConverterTools):
|
||||
-------
|
||||
nkpt : integer
|
||||
total number of k-points in the mesh
|
||||
k_mesh : numpy.array[nkpt,3] of floats
|
||||
kpts : numpy.array[nkpt,3] of floats
|
||||
the coordinates of all k-points
|
||||
wk : numpy.array[nkpt] of floats
|
||||
the weight of each k-point
|
||||
@ -614,30 +829,26 @@ class Wannier90Converter(ConverterTools):
|
||||
if mmode != 0:
|
||||
raise ValueError("Mesh generation mode not supported: %s" % mmode)
|
||||
|
||||
assert len(msize) == 3
|
||||
|
||||
# a regular mesh including Gamma point
|
||||
# total number of k-points
|
||||
nkpt = msize[0] * msize[1] * msize[2]
|
||||
kmesh = numpy.zeros((nkpt, 3), dtype=float)
|
||||
ii = 0
|
||||
for ix, iy, iz in product(list(range(msize[0])), list(range(msize[1])), list(range(msize[2]))):
|
||||
kmesh[ii, :] = [float(ix) / msize[0], float(iy) /
|
||||
msize[1], float(iz) / msize[2]]
|
||||
ii += 1
|
||||
nkpt = numpy.prod(msize)
|
||||
kpts = numpy.array(list(product(range(msize[0]), range(msize[1]), range(msize[2]))))
|
||||
kpts = kpts / numpy.array(msize)
|
||||
# weight is equal for all k-points because wannier90 uses uniform grid on whole BZ
|
||||
# (normalization is always 1 and takes into account spin degeneracy)
|
||||
wk = numpy.ones([nkpt], dtype=float) / float(nkpt)
|
||||
wk = numpy.full(nkpt, 1/nkpt)
|
||||
|
||||
return nkpt, kmesh, wk
|
||||
return nkpt, kpts, wk
|
||||
|
||||
def fourier_ham(self, norb, h_of_r):
|
||||
def fourier_ham(self, h_of_r):
|
||||
"""
|
||||
Method for obtaining H(k) from H(R) via Fourier transform
|
||||
The R vectors and k-point mesh are read from global module variables
|
||||
|
||||
Parameters
|
||||
----------
|
||||
norb : integer
|
||||
number of orbitals
|
||||
h_of_r : list of numpy.array[norb,norb]
|
||||
Hamiltonian H(R) in Wannier basis
|
||||
|
||||
@ -648,14 +859,122 @@ class Wannier90Converter(ConverterTools):
|
||||
|
||||
"""
|
||||
|
||||
twopi = 2 * numpy.pi
|
||||
h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_)
|
||||
h_of_k = [numpy.zeros((self.nwfs, self.nwfs), dtype=complex)
|
||||
for ik in range(self.n_k)]
|
||||
ridx = numpy.array(list(range(self.nrpt)))
|
||||
for ik, ir in product(list(range(self.n_k)), ridx):
|
||||
rdotk = twopi * numpy.dot(self.k_mesh[ik], self.rvec[ir])
|
||||
factor = (math.cos(rdotk) + 1j * math.sin(rdotk)) / \
|
||||
float(self.rdeg[ir])
|
||||
rdotk = numpy.dot(self.kpts[ik], self.rvec[ir])
|
||||
factor = numpy.exp(2j * numpy.pi * rdotk) / self.rdeg[ir]
|
||||
h_of_k[ik][:, :] += factor * h_of_r[ir][:, :]
|
||||
|
||||
return h_of_k
|
||||
|
||||
def convert_misc_input(self, out_filename, nnkp_filename, n_spin_blocs):
|
||||
"""
|
||||
Reads input from DFT code calculations to get occupations
|
||||
|
||||
Parameters
|
||||
----------
|
||||
out_filename : string
|
||||
filename of DFT output file containing occupation data
|
||||
nnkp_file : string
|
||||
filename of Wannier postproc_setup run
|
||||
n_spin_blocs : int
|
||||
SP + 1 - SO
|
||||
|
||||
Returns
|
||||
-------
|
||||
fermi_weights : numpy.array[self.n_k, n_spin_blocs ,n_orbitals]
|
||||
occupations from DFT calculation
|
||||
band_window : numpy.array[self.n_k, n_spin_blocs ,n_orbitals]
|
||||
band indices of correlated subspace
|
||||
fermi_energy : float
|
||||
the DFT Kohn-Sham Fermi energy
|
||||
"""
|
||||
|
||||
# Read only from the master node
|
||||
if not (mpi.is_master_node()):
|
||||
return
|
||||
|
||||
assert n_spin_blocs == 1, 'spin-polarized not implemented'
|
||||
assert 'nscf.out' in out_filename or out_filename == 'LOCPROJ'
|
||||
|
||||
occupations = []
|
||||
if 'nscf.out' in out_filename:
|
||||
occupations = []
|
||||
with open(out_filename,'r') as out_file:
|
||||
out_data = out_file.readlines()
|
||||
# Reads number of Kohn-Sham states and total number of electrons
|
||||
for line in out_data:
|
||||
if 'number of electrons' in line:
|
||||
n_total_electrons = float(line.split()[-1])
|
||||
|
||||
if 'number of Kohn-Sham states' in line:
|
||||
n_ks = int(line.split()[-1])
|
||||
|
||||
if 'Fermi energy' in line:
|
||||
fermi_energy = float(line.split()[-2])
|
||||
|
||||
# get occupations
|
||||
for ct, line in enumerate(out_data):
|
||||
if line.strip() == 'End of band structure calculation':
|
||||
break
|
||||
|
||||
assert 'k = ' in out_data[ct + 2], 'Cannot read occupations. Set verbosity = "high" in {}'.format(out_filename)
|
||||
out_data = out_data[ct+2:]
|
||||
|
||||
# block size of eigenvalues + occupations per k-point
|
||||
n_block = int(2*numpy.ceil(n_ks/8)+5)
|
||||
|
||||
for ik in range(self.n_k):
|
||||
# get data
|
||||
k_block = [line.split() for line in out_data[ik*n_block+2:ik*n_block+n_block-1]]
|
||||
# second half corresponds to occupations
|
||||
occs = k_block[int(len(k_block)/2)+1:]
|
||||
flattened_occs = [float(item) for sublist in occs for item in sublist]
|
||||
occupations.append(flattened_occs)
|
||||
else:
|
||||
# Reads LOCPROJ
|
||||
with open(out_filename, 'r') as file:
|
||||
header = file.readline()
|
||||
n_ks = int(header.split()[2])
|
||||
fermi_energy = float(header.split()[4])
|
||||
|
||||
occupations = numpy.loadtxt((line for line in file if 'orbital' in line), usecols=5)
|
||||
occupations = occupations.reshape((self.n_k, n_ks))
|
||||
|
||||
# assume that all bands contribute, then remove from exclude_bands; python indexing
|
||||
corr_bands = list(range(n_ks))
|
||||
# read exclude_bands from "seedname.nnkp" file
|
||||
with open(nnkp_filename, 'r') as nnkp_file:
|
||||
read = False
|
||||
skip = False
|
||||
for line in nnkp_file:
|
||||
if line.strip() == 'begin exclude_bands':
|
||||
read = True
|
||||
# skip one more line that contains total number of excluded bands
|
||||
skip = True
|
||||
continue
|
||||
elif line.strip() == 'end exclude_bands':
|
||||
read = False
|
||||
continue
|
||||
elif skip:
|
||||
skip = False
|
||||
continue
|
||||
elif read:
|
||||
# wannier index -1
|
||||
corr_bands.remove(int(line)-1)
|
||||
|
||||
# For now, it is only supported to exclude the lowest and highest bands
|
||||
# If bands in the middle are supposed to be excluded, this doesn't work with the band_window
|
||||
# We'd need to manually add rows of zeros to the projectors
|
||||
if numpy.any(numpy.diff(corr_bands) != 1):
|
||||
raise NotImplementedError('Can only exclude the lowest or highest bands')
|
||||
band_window = numpy.array([[(min(corr_bands), max(corr_bands))]*self.n_k]*n_spin_blocs)
|
||||
|
||||
included_occupations = numpy.array(occupations)[:, corr_bands]
|
||||
# Adds spin dimension without changing the array
|
||||
f_weights = included_occupations.reshape(included_occupations.shape[0], 1,
|
||||
included_occupations.shape[1])
|
||||
|
||||
return f_weights, band_window, fermi_energy
|
||||
|
Binary file not shown.
Binary file not shown.
Loading…
Reference in New Issue
Block a user