diff --git a/python/triqs_dft_tools/converters/wannier90.py b/python/triqs_dft_tools/converters/wannier90.py index 7f4cf8ab..0cadcf91 100644 --- a/python/triqs_dft_tools/converters/wannier90.py +++ b/python/triqs_dft_tools/converters/wannier90.py @@ -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) = # 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 = 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 diff --git a/test/python/w90_convert_hloc_diag.ref.h5 b/test/python/w90_convert_hloc_diag.ref.h5 index 6beacdd6..02faef6c 100644 Binary files a/test/python/w90_convert_hloc_diag.ref.h5 and b/test/python/w90_convert_hloc_diag.ref.h5 differ diff --git a/test/python/w90_convert_wannier.ref.h5 b/test/python/w90_convert_wannier.ref.h5 index 3829a31c..6f1a2268 100644 Binary files a/test/python/w90_convert_wannier.ref.h5 and b/test/python/w90_convert_wannier.ref.h5 differ