Source code for triqs_dft_tools.converters.wien2k


##########################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
##########################################################################
"""
Wien2k converter
"""

from types import *
import numpy
from h5 import *
from .converter_tools import *
import os.path


[docs] class Wien2kConverter(ConverterTools): """ Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkDFT class. """
[docs] def __init__(self, filename, hdf_filename=None, dft_subgrp='dft_input', symmcorr_subgrp='dft_symmcorr_input', parproj_subgrp='dft_parproj_input', symmpar_subgrp='dft_symmpar_input', bands_subgrp='dft_bands_input', misc_subgrp='dft_misc_input', transp_subgrp='dft_transp_input', repacking=False): """ Initialise the class. Parameters ---------- filename : string Base name of DFT files. hdf_filename : string, optional Name of hdf5 archive to be created. dft_subgrp : string, optional Name of subgroup storing necessary DFT data. symmcorr_subgrp : string, optional Name of subgroup storing correlated-shell symmetry data. parproj_subgrp : string, optional Name of subgroup storing partial projector data. symmpar_subgrp : string, optional Name of subgroup storing partial-projector symmetry data. bands_subgrp : string, optional Name of subgroup storing band data. misc_subgrp : string, optional Name of subgroup storing miscellaneous DFT data. transp_subgrp : string, optional Name of subgroup storing transport data. repacking : boolean, optional Does the hdf5 archive need to be repacked to save space? """ assert isinstance(filename, str), "Wien2kConverter: Please provide the DFT files' base name as a string." if hdf_filename is None: hdf_filename = filename + '.h5' self.hdf_file = hdf_filename self.dft_file = filename + '.ctqmcout' self.symmcorr_file = filename + '.symqmc' self.parproj_file = filename + '.parproj' self.symmpar_file = filename + '.sympar' self.band_file = filename + '.outband' self.bandwin_file = filename + '.oubwin' self.struct_file = filename + '.struct' self.outputs_file = filename + '.outputs' self.pmat_file = filename + '.pmat' self.dft_subgrp = dft_subgrp self.symmcorr_subgrp = symmcorr_subgrp self.parproj_subgrp = parproj_subgrp self.symmpar_subgrp = symmpar_subgrp self.bands_subgrp = bands_subgrp self.misc_subgrp = misc_subgrp self.transp_subgrp = transp_subgrp self.fortran_to_replace = {'D': 'E'} # Checks if h5 file is there and repacks it if wanted: if (os.path.exists(self.hdf_file) and repacking): ConverterTools.repack(self)
[docs] def convert_dft_input(self): """ Reads the appropriate files and stores the data for the - dft_subgrp - symmcorr_subgrp - misc_subgrp 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.dft_file) # R is a generator : each R.Next() will return the next number in the # file R = ConverterTools.read_fortran_file( self, self.dft_file, self.fortran_to_replace) try: energy_unit = next(R) # read the energy convertion factor # read the number of k points n_k = int(next(R)) k_dep_projection = 1 # flag for spin-polarised calculation SP = int(next(R)) # flag for spin-orbit calculation SO = int(next(R)) charge_below = next(R) # total charge below energy window # total density required, for setting the chemical potential density_required = next(R) symm_op = 1 # Use symmetry groups for the k-sum # the information on the non-correlated shells is not important # here, maybe skip: # number of shells (e.g. Fe d, As p, O p) in the unit cell, n_shells = int(next(R)) # corresponds to index R in formulas # now read the information about the shells (atom, sort, l, dim): shell_entries = ['atom', 'sort', 'l', 'dim'] shells = [{name: int(val) for name, val in zip( shell_entries, R)} for ish in range(n_shells)] # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, n_corr_shells = int(next(R)) # corresponds to index R in formulas # now read the information about the shells (atom, sort, l, dim, SO # flag, irep): corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep'] corr_shells = [{name: int(val) for name, val in zip( corr_shell_entries, R)} for icrsh in range(n_corr_shells)] # determine the number of inequivalent correlated shells and maps, # needed for further reading n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence( self, corr_shells) use_rotations = 1 rot_mat = [numpy.identity( corr_shells[icrsh]['dim'], complex) for icrsh in range(n_corr_shells)] # read the matrices rot_mat_time_inv = [0 for i in range(n_corr_shells)] for icrsh in range(n_corr_shells): for i in range(corr_shells[icrsh]['dim']): # read real part: for j in range(corr_shells[icrsh]['dim']): rot_mat[icrsh][i, j] = next(R) # read imaginary part: for i in range(corr_shells[icrsh]['dim']): for j in range(corr_shells[icrsh]['dim']): rot_mat[icrsh][i, j] += 1j * next(R) if (SP == 1): # read time inversion flag: rot_mat_time_inv[icrsh] = int(next(R)) # Read here the info for the transformation of the basis: n_reps = [1 for i in range(n_inequiv_shells)] dim_reps = [0 for i in range(n_inequiv_shells)] T = [] for ish in range(n_inequiv_shells): # number of representatives ("subsets"), e.g. t2g and eg n_reps[ish] = int(next(R)) dim_reps[ish] = [int(next(R)) for i in range( n_reps[ish])] # dimensions of the subsets # The transformation matrix: # is of dimension 2l+1 without SO, and 2*(2l+1) with SO! 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], complex)) # now read it from file: for i in range(lmax): for j in range(lmax): T[ish][i, j] = next(R) for i in range(lmax): for j in range(lmax): T[ish][i, j] += 1j * next(R) # Spin blocks to be read: n_spin_blocs = SP + 1 - SO # read the list of n_orbitals for all k points n_orbitals = numpy.zeros([n_k, n_spin_blocs], int) for isp in range(n_spin_blocs): for ik in range(n_k): n_orbitals[ik, isp] = int(next(R)) # Initialise the projectors: proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max( [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex) # Read the projectors from the file: for ik in range(n_k): for icrsh in range(n_corr_shells): n_orb = corr_shells[icrsh]['dim'] # first Real part for BOTH spins, due to conventions in # dmftproj: for isp in range(n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik][isp]): proj_mat[ik, isp, icrsh, i, j] = next(R) # now Imag part: for isp in range(n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik][isp]): proj_mat[ik, isp, icrsh, i, j] += 1j * next(R) # now define the arrays for weights and hopping ... # w(k_index), default normalisation bz_weights = numpy.ones([n_k], float) / float(n_k) hopping = numpy.zeros([n_k, n_spin_blocs, numpy.max( n_orbitals), numpy.max(n_orbitals)], complex) # weights in the file for ik in range(n_k): bz_weights[ik] = next(R) # if the sum over spins is in the weights, take it out again!! sm = sum(bz_weights) bz_weights[:] /= sm # Grab the H # we use now the convention of a DIAGONAL Hamiltonian -- convention # for Wien2K. for isp in range(n_spin_blocs): for ik in range(n_k): n_orb = n_orbitals[ik, isp] for i in range(n_orb): hopping[ik, isp, i, i] = next(R) * energy_unit # keep some things that we need for reading parproj: things_to_set = ['n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'n_spin_blocs', 'n_orbitals', 'n_k', 'SO', 'SP', 'energy_unit'] for it in things_to_set: setattr(self, it, locals()[it]) except StopIteration: # a more explicit error if the file is corrupted. raise IOError("wien2k : reading file %s failed!" % self.dft_file) R.close() # Reading done! #new variable: dft_code - this determines which DFT code the inputs come from. #used for certain routines within dft_tools if treating the inputs differently is required. dft_code = 'wien2k' # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp) # The subgroup containing the data. If it does not exist, it is # created. If it exists, the data is overwritten! 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', 'dft_code'] for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it] # Symmetries are used, so now convert symmetry information for # *correlated* orbitals: self.convert_symmetry_input(orbits=self.corr_shells, symm_file=self.symmcorr_file, symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP) self.convert_misc_input()
[docs] def convert_parproj_input(self): """ Reads the appropriate files and stores the data for the - parproj_subgrp - symmpar_subgrp in the hdf5 archive. """ if not (mpi.is_master_node()): return # get needed data from hdf file with HDFArchive(self.hdf_file, 'a') as ar: things_to_read = ['SP', 'SO', 'n_shells', 'n_k', 'n_orbitals', 'shells'] for it in things_to_read: if not hasattr(self, it): setattr(self, it, ar[self.dft_subgrp][it]) self.n_spin_blocs = self.SP + 1 - self.SO mpi.report("Reading input from %s..." % self.parproj_file) dens_mat_below = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], complex) for ish in range(self.n_shells)] for isp in range(self.n_spin_blocs)] R = ConverterTools.read_fortran_file( self, self.parproj_file, self.fortran_to_replace) n_parproj = [int(next(R)) for i in range(self.n_shells)] n_parproj = numpy.array(n_parproj) # Initialise P, here a double list of matrices: proj_mat_all = numpy.zeros([self.n_k, self.n_spin_blocs, self.n_shells, max( n_parproj), max([sh['dim'] for sh in self.shells]), numpy.max(self.n_orbitals)], complex) rot_mat_all = [numpy.identity( self.shells[ish]['dim'], complex) for ish in range(self.n_shells)] rot_mat_all_time_inv = [0 for i in range(self.n_shells)] for ish in range(self.n_shells): # read first the projectors for this orbital: for ik in range(self.n_k): for ir in range(n_parproj[ish]): for isp in range(self.n_spin_blocs): # read real part: for i in range(self.shells[ish]['dim']): for j in range(self.n_orbitals[ik][isp]): proj_mat_all[ik, isp, ish, ir, i, j] = next(R) for isp in range(self.n_spin_blocs): # read imaginary part: for i in range(self.shells[ish]['dim']): for j in range(self.n_orbitals[ik][isp]): proj_mat_all[ik, isp, ish, ir, i, j] += 1j * next(R) # now read the Density Matrix for this orbital below the energy # window: for isp in range(self.n_spin_blocs): for i in range(self.shells[ish]['dim']): # read real part: for j in range(self.shells[ish]['dim']): dens_mat_below[isp][ish][i, j] = next(R) for isp in range(self.n_spin_blocs): # read imaginary part: for i in range(self.shells[ish]['dim']): for j in range(self.shells[ish]['dim']): dens_mat_below[isp][ish][i, j] += 1j * next(R) if (self.SP == 0): dens_mat_below[isp][ish] /= 2.0 # Global -> local rotation matrix for this shell: for i in range(self.shells[ish]['dim']): # read real part: for j in range(self.shells[ish]['dim']): rot_mat_all[ish][i, j] = next(R) for i in range(self.shells[ish]['dim']): # read imaginary part: for j in range(self.shells[ish]['dim']): rot_mat_all[ish][i, j] += 1j * next(R) if (self.SP): rot_mat_all_time_inv[ish] = int(next(R)) R.close() # Reading done! # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (self.parproj_subgrp in ar): ar.create_group(self.parproj_subgrp) # The subgroup containing the data. If it does not exist, it is # created. If it exists, the data is overwritten! things_to_save = ['dens_mat_below', 'n_parproj', 'proj_mat_all', 'rot_mat_all', 'rot_mat_all_time_inv'] for it in things_to_save: ar[self.parproj_subgrp][it] = locals()[it] # Symmetries are used, so now convert symmetry information for *all* # orbitals: self.convert_symmetry_input(orbits=self.shells, symm_file=self.symmpar_file, symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP)
[docs] def convert_bands_input(self): """ Reads the appropriate files and stores the data for the bands_subgrp in the hdf5 archive. """ if not (mpi.is_master_node()): return try: # get needed data from hdf file with HDFArchive(self.hdf_file, 'a') as ar: things_to_read = ['SP', 'SO', 'n_corr_shells', 'n_shells', 'corr_shells', 'shells', 'energy_unit'] for it in things_to_read: if not hasattr(self, it): setattr(self, it, ar[self.dft_subgrp][it]) self.n_spin_blocs = self.SP + 1 - self.SO mpi.report("Reading input from %s..." % self.band_file) R = ConverterTools.read_fortran_file( self, self.band_file, self.fortran_to_replace) n_k = int(next(R)) # read the list of n_orbitals for all k points n_orbitals = numpy.zeros([n_k, self.n_spin_blocs], int) for isp in range(self.n_spin_blocs): for ik in range(n_k): n_orbitals[ik, isp] = int(next(R)) # Initialise the projectors: proj_mat = numpy.zeros([n_k, self.n_spin_blocs, self.n_corr_shells, max( [crsh['dim'] for crsh in self.corr_shells]), numpy.max(n_orbitals)], complex) # Read the projectors from the file: for ik in range(n_k): for icrsh in range(self.n_corr_shells): n_orb = self.corr_shells[icrsh]['dim'] # first Real part for BOTH spins, due to conventions in # dmftproj: for isp in range(self.n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik, isp]): proj_mat[ik, isp, icrsh, i, j] = next(R) # now Imag part: for isp in range(self.n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik, isp]): proj_mat[ik, isp, icrsh, i, j] += 1j * next(R) hopping = numpy.zeros([n_k, self.n_spin_blocs, numpy.max( n_orbitals), numpy.max(n_orbitals)], complex) # Grab the H # we use now the convention of a DIAGONAL Hamiltonian!!!! for isp in range(self.n_spin_blocs): for ik in range(n_k): n_orb = n_orbitals[ik, isp] for i in range(n_orb): hopping[ik, isp, i, i] = next(R) * self.energy_unit # now read the partial projectors: n_parproj = [int(next(R)) for i in range(self.n_shells)] n_parproj = numpy.array(n_parproj) # Initialise P, here a double list of matrices: proj_mat_all = numpy.zeros([n_k, self.n_spin_blocs, self.n_shells, max(n_parproj), max( [sh['dim'] for sh in self.shells]), numpy.max(n_orbitals)], complex) for ish in range(self.n_shells): for ik in range(n_k): for ir in range(n_parproj[ish]): for isp in range(self.n_spin_blocs): # read real part: for i in range(self.shells[ish]['dim']): for j in range(n_orbitals[ik, isp]): proj_mat_all[ik, isp, ish, ir, i, j] = next(R) # read imaginary part: for i in range(self.shells[ish]['dim']): for j in range(n_orbitals[ik, isp]): proj_mat_all[ik, isp, ish, ir, i, j] += 1j * next(R) R.close() except KeyError: raise IOError("convert_bands_input : Needed data not found in hdf file. Consider calling convert_dft_input first!") except StopIteration: # a more explicit error if the file is corrupted. raise IOError("wien2k : reading file %s failed!" % self.band_file) # Reading done! # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (self.bands_subgrp in ar): ar.create_group(self.bands_subgrp) # The subgroup containing the data. If it does not exist, it is # created. If it exists, the data is overwritten! things_to_save = ['n_k', 'n_orbitals', 'proj_mat', 'hopping', 'n_parproj', 'proj_mat_all'] for it in things_to_save: ar[self.bands_subgrp][it] = locals()[it]
[docs] def convert_misc_input(self): """ Reads additional information on: - the band window from :file:`case.oubwin`, - lattice parameters from :file:`case.struct`, - symmetries from :file:`case.outputs`, if those Wien2k files are present and stores the data in the hdf5 archive. This function is automatically called by :meth:`convert_dft_input <triqs_dft_tools.converters.wien2k.Wien2kConverter.convert_dft_input>`. """ if not (mpi.is_master_node()): return # Check if SP, SO and n_k are already in h5 with HDFArchive(self.hdf_file, 'r') as ar: if not (self.dft_subgrp in ar): raise IOError("convert_misc_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp) SP = ar[self.dft_subgrp]['SP'] SO = ar[self.dft_subgrp]['SO'] n_k = ar[self.dft_subgrp]['n_k'] things_to_save = [] # Read relevant data from .oubwin/up/dn files ############################################# # band_window: Contains the index of the lowest and highest band within the # projected subspace (used by dmftproj) for each k-point. if (SP == 0 and SO == 0): # read .oubwin file files = [self.bandwin_file] elif (SP == 1 and SO == 0): # read .oubwinup and .oubwindn files = [self.bandwin_file + 'up', self.bandwin_file + 'dn'] elif (SP == 1 and SO == 1): # read either .oubwinup or .oubwindn if os.path.exists(self.bandwin_file + 'up'): files = [self.bandwin_file + 'up'] elif os.path.exists(self.bandwin_file + 'dn'): files = [self.bandwin_file + 'dn'] else: assert 0, "convert_misc_input: If SO and SP are 1 provide either .oubwinup or .oubwindn file" else: assert 0, "convert_misc_input: Reading oubwin error! Check SP and SO, if SO=1 SP must be 1." band_window = [None for isp in range(SP + 1 - SO)] for isp, f in enumerate(files): if os.path.exists(f): mpi.report("Reading input from %s..." % f) R = ConverterTools.read_fortran_file( self, f, self.fortran_to_replace) n_k_oubwin = int(next(R)) if (n_k_oubwin != n_k): mpi.report( "convert_misc_input : WARNING : n_k in case.oubwin is different from n_k in case.klist") assert int( next(R)) == SO, "convert_misc_input: SO is inconsistent in oubwin file!" band_window[isp] = numpy.zeros((n_k_oubwin, 2), dtype=int) for ik in range(n_k_oubwin): next(R) band_window[isp][ik, 0] = next(R) # lowest band band_window[isp][ik, 1] = next(R) # highest band next(R) things_to_save.append('band_window') R.close() # Reading done! # Read relevant data from .struct file ###################################### # lattice_type: bravais lattice type as defined by Wien2k # lattice_constants: unit cell parameters in a. u. # lattice_angles: unit cell angles in rad if (os.path.exists(self.struct_file)): mpi.report("Reading input from %s..." % self.struct_file) with open(self.struct_file) as R: try: R.readline() lattice_type = R.readline().split()[0] R.readline() temp = R.readline() lattice_constants = numpy.array( [float(temp[0 + 10 * i:10 + 10 * i].strip()) for i in range(3)]) lattice_angles = numpy.array( [float(temp[30 + 10 * i:40 + 10 * i].strip()) for i in range(3)]) * numpy.pi / 180.0 things_to_save.extend( ['lattice_type', 'lattice_constants', 'lattice_angles']) except IOError: raise IOError("convert_misc_input: reading file %s failed" % self.struct_file) # Read relevant data from .outputs file ####################################### # rot_symmetries: matrix representation of all (space group) symmetry # operations if (os.path.exists(self.outputs_file)): mpi.report("Reading input from %s..." % self.outputs_file) rot_symmetries = [] with open(self.outputs_file) as R: try: while 1: temp = R.readline().strip(' ').split() if (temp[0] == 'PGBSYM:'): n_symmetries = int(temp[-1]) break for i in range(n_symmetries): while 1: if (R.readline().strip().split()[0] == 'Symmetry'): break sym_i = numpy.zeros((3, 3), dtype=float) for ir in range(3): temp = R.readline().strip().split() for ic in range(3): sym_i[ir, ic] = float(temp[ic]) R.readline() rot_symmetries.append(sym_i) things_to_save.extend(['n_symmetries', 'rot_symmetries']) things_to_save.append('rot_symmetries') except IOError: raise IOError("convert_misc_input: reading file %s failed" % self.outputs_file) # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (self.misc_subgrp in ar): ar.create_group(self.misc_subgrp) for it in things_to_save: ar[self.misc_subgrp][it] = locals()[it]
[docs] def convert_transport_input(self): """ Reads the necessary information for transport calculations on: - the optical band window and the velocity matrix elements from :file:`case.pmat` and stores the data in the hdf5 archive. """ if not (mpi.is_master_node()): return # Check if SP, SO and n_k are already in h5 with HDFArchive(self.hdf_file, 'r') as ar: if not (self.dft_subgrp in ar): raise IOError("convert_transport_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp) SP = ar[self.dft_subgrp]['SP'] SO = ar[self.dft_subgrp]['SO'] n_k = ar[self.dft_subgrp]['n_k'] # Read relevant data from .pmat/up/dn files ########################################### # band_window_optics: Contains the index of the lowest and highest band within the # band window (used by optics) for each k-point. # velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics # and each k-point. if (SP == 0 and SO == 0): # read .pmat file files = [self.pmat_file] elif (SP == 1 and SO == 0): # read .pmatup and pmatdn files = [self.pmat_file + 'up', self.pmat_file + 'dn'] elif (SP == 1 and SO == 1): # read either .pmatup or .pmatdn if os.path.exists(self.pmat_file + 'up'): files = [self.pmat_file + 'up'] elif os.path.exists(self.pmat_file + 'dn'): files = [self.pmat_file + 'dn'] else: assert 0, "convert_transport_input: If SO and SP are 1 provide either .pmatup or .pmatdn file" else: assert 0, "convert_transport_input: Reading velocity file error! Check SP and SO, if SO=1 SP must be 1." velocities_k = [[] for f in files] band_window_optics = [] for isp, f in enumerate(files): if not os.path.exists(f): raise IOError("convert_transport_input: File %s does not exist" % f) mpi.report("Reading input from %s..." % f) R = ConverterTools.read_fortran_file( self, f, {'D': 'E', '(': '', ')': '', ',': ' '}) band_window_optics_isp = [] for ik in range(n_k): next(R) nu1 = int(next(R)) nu2 = int(next(R)) band_window_optics_isp.append((nu1, nu2)) n_bands = nu2 - nu1 + 1 for _ in range(4): next(R) if n_bands <= 0: velocity_xyz = numpy.zeros((1, 1, 3), dtype=complex) else: velocity_xyz = numpy.zeros( (n_bands, n_bands, 3), dtype=complex) for nu_i in range(n_bands): for nu_j in range(nu_i, n_bands): for i in range(3): velocity_xyz[nu_i][nu_j][ i] = next(R) + next(R) * 1j if (nu_i != nu_j): velocity_xyz[nu_j][nu_i][i] = velocity_xyz[ nu_i][nu_j][i].conjugate() velocities_k[isp].append(velocity_xyz) band_window_optics.append(numpy.array(band_window_optics_isp)) # check the number of kpoints is correct last_kpt = None try: last_kpt = next(R) except: pass if last_kpt is not None: raise ValueError("The number of kpoints in case.pmat is larger than number of kpoints of Green's function. Please rerun dmftproj") R.close() # Reading done! # Put data to HDF5 file with HDFArchive(self.hdf_file, 'a') as ar: if not (self.transp_subgrp in ar): ar.create_group(self.transp_subgrp) # The subgroup containing the data. If it does not exist, it is # created. If it exists, the data is overwritten!!! things_to_save = ['band_window_optics', 'velocities_k'] for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it]
[docs] def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): """ Reads and stores symmetrisation data from symm_file, which can be is case.sympar or case.symqmc. Parameters ---------- orbits : list of dicts This is either shells or corr_shells depending on whether the symmetry information is for correlated shells or partial projectors. symm_file : string Name of the file containing symmetry data. This is case.symqmc for correlated shells and case.sympar for partial projectors. symm_subgrp : string, optional Name of subgroup storing symmetry data. SO : integer Is spin-orbit coupling considered? SP : integer Is the system spin-polarised? """ if not (mpi.is_master_node()): return mpi.report("Reading input from %s..." % symm_file) n_orbits = len(orbits) R = ConverterTools.read_fortran_file( self, symm_file, self.fortran_to_replace) try: n_symm = int(next(R)) # Number of symmetry operations n_atoms = int(next(R)) # number of atoms involved perm = [[int(next(R)) for i in range(n_atoms)] for j in range(n_symm)] # list of permutations of the atoms if SP: # time inversion for SO coupling time_inv = [int(next(R)) for j in range(n_symm)] else: time_inv = [0 for j in range(n_symm)] # Now read matrices: mat = [] for i_symm in range(n_symm): mat.append([numpy.zeros([orbits[orb]['dim'], orbits[orb][ 'dim']], complex) for orb in range(n_orbits)]) for orb in range(n_orbits): for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): # real part mat[i_symm][orb][i, j] = next(R) for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): mat[i_symm][orb][i, j] += 1j * \ next(R) # imaginary part mat_tinv = [numpy.identity(orbits[orb]['dim'], complex) for orb in range(n_orbits)] if ((SO == 0) and (SP == 0)): # here we need an additional time inversion operation, so read # it: for orb in range(n_orbits): for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): # real part mat_tinv[orb][i, j] = next(R) for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): mat_tinv[orb][i, j] += 1j * \ next(R) # imaginary part except StopIteration: # a more explicit error if the file is corrupted. raise IOError("wien2k : reading file %s failed!" %symm_file) R.close() # Reading done! # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (symm_subgrp in ar): ar.create_group(symm_subgrp) things_to_save = ['n_symm', 'n_atoms', 'perm', 'orbits', 'SO', 'SP', 'time_inv', 'mat', 'mat_tinv'] for it in things_to_save: ar[symm_subgrp][it] = locals()[it]