Source code for triqs_dftkit.wien2k.converter
##########################################################################
#
# 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
def _parse_kpath_label_line(line):
"""
Parse one line of the high-symmetry point block that dmftproj appends to
case.outband. The format is the fixed Fortran (2i6,a): columns 1-6 hold a
running counter, columns 7-12 the 1-based position of the point along the
band path, and columns 13 onwards the label.
Returns (counter, pos, label), or None if the line does not have that form.
"""
line = line.rstrip('\n')
if len(line) < 13:
return None
try:
counter = int(line[0:6])
pos = int(line[6:12])
except ValueError:
return None
label = line[12:].strip()
if not label:
return None
return counter, pos, label
def _read_kpath_labels(band_file, n_k):
"""
Read the high-symmetry k-path labels appended to the end of case.outband
and map them onto the flattened band k-point index.
dmftproj writes one line per high-symmetry point after the projector data,
in fixed Fortran format (2i6,a), e.g.
1 1GAMMA
2 122X
Returns (labels, idx) where labels is a list of label strings and idx is a
0-based numpy int array giving, for each label, the position of that
high-symmetry point in the n_k band path. Returns (None, None) if no label
block is present or if the block is incomplete.
"""
# The block is a handful of lines at the very end of a file that holds all
# the projectors, so read the tail rather than the whole file.
n_bytes = 8192
with open(band_file, 'rb') as R:
R.seek(0, os.SEEK_END)
size = R.tell()
R.seek(max(0, size - n_bytes))
lines = R.read().decode('utf-8', 'replace').splitlines()
if size > n_bytes:
# the first line of the chunk is in general cut in the middle
lines = lines[1:]
while lines and not lines[-1].strip():
lines.pop()
# Walk backwards from the end of the file: the label block is the trailing
# run of (2i6,a) lines, and the counter of its first line is 1.
labels = []
idx = []
first_counter = None
for line in reversed(lines):
parsed = _parse_kpath_label_line(line)
if parsed is None:
break
first_counter, pos, label = parsed
labels.append(label)
idx.append(pos - 1)
if first_counter == 1:
break
labels.reverse()
idx.reverse()
if not labels:
return None, None
if first_counter != 1:
mpi.report("convert_bands_input : WARNING : the high-symmetry point block in %s does not start at counter 1, so it is truncated or corrupted; skipping k-path labels." % band_file)
return None, None
idx = numpy.array(idx, dtype=int)
if idx[0] < 0 or idx[-1] >= n_k or numpy.any(numpy.diff(idx) < 0):
mpi.report("convert_bands_input : WARNING : inconsistent high-symmetry point indices in %s; skipping k-path labels." % band_file)
return None, None
return labels, idx
[docs]
class Converter(ConverterTools):
"""
Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkDFT class.
"""
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!
kpts_labels, kpts_labels_idx = _read_kpath_labels(self.band_file, n_k)
# 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']
if kpts_labels is not None:
things_to_save += ['kpts_labels', 'kpts_labels_idx']
mpi.report(" Stored %i high-symmetry k-path labels: %s" % (len(kpts_labels), ', '.join(kpts_labels)))
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_dftkit.wien2k.Converter.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]