LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-23 21:03:45 +01:00
Code Issues Releases Wiki Activity
f134fc8fa3
dft_tools/python/triqs_dft_tools/converters/elk.py
Nils Wentzell f134fc8fa3 Replace numpy.lib.pad with numpy.pad
2024-06-21 12:01:14 -04:00

854 lines
39 KiB
Python
Raw Blame History

##########################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2019 by A. D. N. James, A. Hampel and M. Aichhorn
#
# 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/>.
#
##########################################################################
"""
Elk converter
"""
from types import *
import numpy
from .converter_tools import *
import os.path
from h5 import *
from locale import atof
from triqs_dft_tools.converters.elktools import readElkfiles as read_Elk
from triqs_dft_tools.converters.elktools import ElkConverterTools as Elk_tools
class ElkConverter(ConverterTools,Elk_tools,read_Elk):
"""
Conversion from Elk 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',
bc_subgrp='dft_bandchar_input', symmpar_subgrp='dft_symmpar_input',
bands_subgrp='dft_bands_input', misc_subgrp='dft_misc_input',
transp_subgrp='dft_transp_input',cont_subgrp='dft_contours_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), "ElkConverter: 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 = 'PROJ.OUT'
self.band_file = 'BAND.OUT'
self.eval_file = 'EIGVAL.OUT'
self.efermi_file = 'EFERMI.OUT'
self.kp_file = 'KPOINTS.OUT'
self.geom_file='GEOMETRY.OUT'
self.dft_subgrp = dft_subgrp
self.symmcorr_subgrp = symmcorr_subgrp
self.bc_subgrp = bc_subgrp
self.symmpar_subgrp = symmpar_subgrp
self.bands_subgrp = bands_subgrp
self.misc_subgrp = misc_subgrp
self.transp_subgrp = transp_subgrp
self.cont_subgrp = cont_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)
def check_dens(self,n_k,nstsv,occ,bz_weights,n_spin_blocs,band_window,SO):
"""
Check the charge density below the correlated energy window and up to the Fermi level
"""
density_required=0.0E-7
charge_below=0.0E-7
#calculate the valence charge and charge below lower energy window bound
#Elk does not use the tetrahedron method when calculating these charges
for ik in range(n_k):
for ist in range(nstsv):
#calculate the charge over all the bands
density_required+=occ[ik][ist]*bz_weights[ik]
for isp in range(n_spin_blocs):
#Convert occ list from elk to two index format for spins
jst=int((isp)*nstsv*0.5)
#Take lowest index in band_window for SO system
if(SO==0):
nst=band_window[isp][ik, 0]-1
else:
band=[band_window[0][ik, 0],band_window[1][ik, 0]]
nst=min(band)-1
#calculate the charge below energy window
for ist in range(jst,nst):
charge_below+=occ[ik][ist]*bz_weights[ik]
#return charges
return(density_required,charge_below)
def rotsym(self,n_shells,shells,n_symm,ind,basis,T,mat):
"""
Rotates the symmetry matrices into basis defined by the T unitary matrix
the outputted projectors are rotated to the irreducible representation
and then reduced in size to the orbitals used to construct the projectors.
"""
for ish in range(n_shells):
#check that the T matrix is not the Identity (i.e. not using spherical
#harmonics).
if(basis[ish]!=0):
#put mat into temporary matrix
temp=mat
#index range of lm values used to create the Wannier projectors
min_ind=numpy.min(ind[ish][:])
max_ind=numpy.max(ind[ish][:])+1
#dimension of lm values used to construct the projectors
dim=shells[ish]['dim']
#loop over all symmetries
for isym in range(n_symm):
#rotate symmetry matrix into basis defined by T
mat[isym][ish]=numpy.matmul(T[ish],mat[isym][ish])
mat[isym][ish]=numpy.matmul(mat[isym][ish],T[ish].conjugate().transpose())
#put desired subset of transformed symmetry matrix into temp matrix for symmetry isym
for id in range(len(ind[ish])):
i=ind[ish][id]
for jd in range(len(ind[ish][:])):
j=ind[ish][jd]
temp[isym][ish][id,jd]=mat[isym][ish][i,j]
#put temp matrix into mat
mat=temp
#reduce size of lm arrays in mat lm dim
for isym in range(n_symm):
dim=shells[ish]['dim']
mat[isym][ish]=mat[isym][ish][:dim,:dim]
return mat
def update_so_quatities(self,n_shells,shells,n_corr_shells,corr_shells,n_inequiv_shells,dim_reps,n_k,n_symm,n_orbitals,proj_mat,T,su2,mat,sym=True):
"""
Changes the array sizes and elements for arrays used in spin-orbit coupled calculations.
"""
#change dim for each shell
for ish in range(n_shells):
shells[ish]['dim'] *= 2
for ish in range(n_corr_shells):
corr_shells[ish]['dim'] *= 2
for ish in range(n_inequiv_shells):
dim_reps[ish]=[2*dim_reps[ish][i] for i in range(len(dim_reps[ish]))]
#Make temporary array of original n_orbitals
n_orbitals_orig=n_orbitals
#Make SO n_orbitals array
#loop over k-points
for ik in range(n_k):
#new orbital array
n_orbitals[ik,0]=max(n_orbitals[ik,:])
#reduce array size
n_orbitals=n_orbitals[:,:1]
#Resize proj_mat, mat, T
#make temporary projector array
proj_mat_tmp = numpy.zeros([n_k, 1, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex)
for ish in range(n_corr_shells):
#update proj_mat
for ik in range(n_k):
#extra array elements in "dim" dimension
size=int(0.5*corr_shells[ish]['dim'])
#put each spinor into tmp array and ensure elements are assigned correctly in case of change of max(n_orbitals)
proj_mat_tmp[ik][0][ish][0:size][0:n_orbitals_orig[ik,0]]=proj_mat[ik][0][ish][0:size][0:n_orbitals_orig[ik,0]]
#put other spinor projectors into extra "dim" elements
proj_mat_tmp[ik][0][ish][size:2*size][0:n_orbitals_orig[ik,1]]=proj_mat[ik][1][ish][0:size][0:n_orbitals_orig[ik,1]]
#update T
#extra array elements in each dimension
size=2*corr_shells[ish]['l']+1
#extend the arrays
T[ish]=numpy.pad(T[ish],((0,size),(0,size)),'constant',constant_values=(0.0))
#make block diagonal
T[ish][size:2*size,size:2*size]=T[ish][0:size,0:size]
#update the symmetries arrays if needed
if(sym):
#update mat - This includes the spin SU(2) matrix for spin-coupled calculations
#size of each quadrant in the lm symmetry array.
size=int(0.5*corr_shells[ish]['dim'])
#temporary spin block array for SU(2) spin operations on mat
spinmat = numpy.zeros([size,2,size,2],complex)
for isym in range(n_symm):
#expand size of array
mat[isym][ish]=numpy.pad(mat[isym][ish],((0,size),(0,size)),'constant',constant_values=(0.0))
#make arraye block diagonal
mat[isym][ish][size:2*size,size:2*size]=mat[isym][ish][0:size,0:size]
#apply SU(2) spin matrices to lm symmetries
#put mat into array of spin blocks
for i1,i2 in numpy.ndindex(2,2):
spinmat[0:size,i1,0:size,i2] = mat[isym][ish][i1*size:(i1+1)*size,i2*size:(i2+1)*size]
#apply the SU(2) spin matrices
for ilm,jlm in numpy.ndindex(size,size):
spinmat[ilm,:,jlm,:] = numpy.dot(su2[isym][:,:],spinmat[ilm,:,jlm,:])
#put spinmat into back in mat format
for i1,i2 in numpy.ndindex(2,2):
mat[isym][ish][i1*size:(i1+1)*size,i2*size:(i2+1)*size] = spinmat[0:size,i1,0:size,i2]
#assign arrays and delete temporary arrays
del proj_mat
proj_mat = proj_mat_tmp
del proj_mat_tmp, n_orbitals_orig
return shells,corr_shells,dim_reps,n_orbitals,proj_mat,T,mat
def sort_dft_eigvalues(self,n_spin_blocs,SO,n_k,n_orbitals,band_window,en,energy_unit):
"""
Rearranges the energy eigenvalue arrays into TRIQS format
"""
hopping = numpy.zeros([n_k, n_spin_blocs, numpy.max(n_orbitals), numpy.max(n_orbitals)], complex)
#loop over spin
for isp in range(n_spin_blocs):
#loop over k-points
for ik in range(n_k):
#loop over bands
for ist in range(0,n_orbitals[ik, isp]):
#converter index for spin polarised Elk indices and take SO into consideration
if(SO==0):
jst=int(band_window[isp][ik, 0]-1+ist)
else:
band=[band_window[0][ik, 0],band_window[1][ik, 0]]
jst=int(min(band)-1+ist)
#correlated window energies
hopping[ik,isp,ist,ist]=en[ik][jst]*energy_unit
return hopping
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
filext='.OUT'
dft_file='PROJ'+filext
mpi.report("Reading %s" % dft_file)
#Energy conversion - Elk uses Hartrees
energy_unit = 27.2113850560 # Elk uses hartrees
#The projectors change size per k-point
k_dep_projection = 1
#Symmetries are used
symm_op = 1 # Use symmetry groups for the k-sum
shells=[]
#read information about projectors calculated in the Elk calculation
[gen_info,n_corr_shells,n_inequiv_shells,corr_to_inequiv,inequiv_to_corr,corr_shells,n_reps,dim_reps,ind,basis,T] = read_Elk.read_proj(self,dft_file)
#get info for HDF5 file from gen_info
n_k=gen_info['n_k']
SP=gen_info['spinpol']-1
#Elk uses spinor wavefunctions. Therefore these two spinor wavefunctions have spin-orbit coupling incorporated in them. Here we read in the spinors
n_spin_blocs = SP + 1
SO=gen_info['SO']
n_atoms=gen_info['natm']
#Elk only calculates Wannier projectors (no theta projectors generated):
n_shells=n_corr_shells
for ish in range(n_shells):
shells.append(corr_shells[ish].copy())
#remove last 2 entries from corr_shlls
del shells[ish]['SO']
del shells[ish]['irep']
shells[ish]['dim'] = 2*shells[ish]['l']+1
#read eigenvalues calculated in the Elk calculation
mpi.report("Reading %s and EFERMI.OUT" % self.eval_file)
[en,occ,nstsv]=read_Elk.read_eig(self)
#read projectors calculated in the Elk calculation
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), nstsv], complex)
mpi.report("Reading projector(s)")
for ish in range(n_corr_shells):
[n_orbitals,band_window,rep,proj_mat]=read_Elk.read_projector(self,corr_shells,n_spin_blocs,ish,proj_mat,ind,T,basis,filext)
#read kpoints calculated in the Elk calculation
mpi.report("Reading %s" % self.kp_file)
[bz_weights,vkl]=read_Elk.read_kpoints(self)
#symmetry matrix
mpi.report("Reading GEOMETRY.OUT")
#read in atom positions, the symmetry operators (in lattice coordinates) and lattice vectors
[ns, na, atpos]=read_Elk.read_geometry(self)
#Read symmetry files
mpi.report("Reading SYMCRYS.OUT")
[n_symm,spinmat,symmat,tr] = read_Elk.readsym(self)
mpi.report("Reading LATTICE.OUT")
[amat,amatinv,bmat,bmatinv,cell_vol] = read_Elk.readlat(self)
#calculating atom permutations
perm = Elk_tools.gen_perm(self,n_symm,ns,na,n_atoms,symmat,tr,atpos)
#determine the cartesian lattice symmetries and the spin axis rotations
#required for the spinors (for SO for now)
su2 = []
symmatc=[]
for isym in range(n_symm):
#convert the lattice symmetry matrices into cartesian coordinates
tmp = numpy.matmul(amat,symmat[isym])
symmatc.append(numpy.matmul(tmp,amatinv))
#convert the spin symmetry matrices into cartesian coordinates
spinmatc = numpy.matmul(amat,spinmat[isym])
spinmatc = numpy.matmul(spinmatc,amatinv)
#calculate the rotation angle and spin axis vector in cartesian coordinates
[v,th] = self.rotaxang(spinmatc[:,:])
#calculate the SU(2) matrix from the angle and spin axis vector
su2.append(self.axangsu2(v,th))
del tmp
#calculating the symmetries in complex harmonics
mat = Elk_tools.symlat_to_complex_harmonics(self,n_symm,n_corr_shells,symmatc,corr_shells)
mat = self.rotsym(n_corr_shells,corr_shells,n_symm,ind,basis,T,mat)
#The reading is done. Some variables may need to change for TRIQS compatibility.
#Alter size of some of the arrays if spin orbit coupling is enabled.
#For SO in Elk, the eigenvalues and eigenvector band indices are in asscending order w.r.t energy
if(SO==1):
[shells,corr_shells,dim_reps,n_orbitals,proj_mat,T,mat]=self.update_so_quatities(n_shells,shells,n_corr_shells,corr_shells,n_inequiv_shells,dim_reps,n_k,n_symm,n_orbitals,proj_mat,T,su2,mat)
#reduce n_spin_blocs
n_spin_blocs = SP + 1 - SO
#put the energy eigenvalues arrays in TRIQS format
hopping = self.sort_dft_eigvalues(n_spin_blocs,SO,n_k,n_orbitals,band_window,en,energy_unit)
#Elk does not use global to local matrix rotation (Rotloc) as is done in Wien2k. However, the projectors
#require a symmetry matrix to rotate from jatom to iatom. Below finds the non inversion
#symmetric matrices which were used in calculating the projectors
use_rotations = 1
rot_mat = [numpy.identity(corr_shells[icrsh]['dim'], complex) for icrsh in range(n_corr_shells)]
for icrsh in range(n_corr_shells):
#return inequivalent index
incrsh = corr_to_inequiv[icrsh]
#return first inequivalent corr_shell index
jcrsh = inequiv_to_corr[incrsh]
#want to rotate atom to first inequivalent atom in list
iatom = corr_shells[jcrsh]['atom']
for isym in range(n_symm):
jatom=perm[isym][corr_shells[icrsh]['atom']-1]
#determinant determines if crystal symmetry matrix has inversion symmetry (=-1)
det = numpy.linalg.det(symmat[isym][:,:])
if((jatom==iatom)&(det>0.0)):
#local rotation which rotates equivalent atom into its local coordinate system
#(inverse of the symmetry operator applied to the projectors in Elk)
rot_mat[icrsh][:,:]=mat[isym][icrsh][:,:].conjugate().transpose()
#used first desired symmetry in crystal symmetry list
break
# Elk does not currently use time inversion symmetry
rot_mat_time_inv = [0 for i in range(n_corr_shells)]
#Check that the charge of all the bands and below the correlated window have been calculated correctly
[density_required, charge_below] = self.check_dens(n_k,nstsv,occ,bz_weights,n_spin_blocs,band_window,SO)
#calculate the required charge (density_required) to remain charge neutral
mpi.report("The total charge of the system = %f" %density_required)
mpi.report("The charge below the correlated window = %f" %charge_below)
mpi.report("The charge within the correlated window = %f" %(density_required - charge_below))
#Elk interface does not calculate theta projectors, hence orbits are the same as Wannier projectors
orbits=[]
#remove the spatom index to avoid errors in the symmetry routines
for ish in range(n_corr_shells):
#remove "spatom"
del corr_shells[ish]['spatom']
orbits.append(corr_shells[ish].copy())
for ish in range(n_shells):
#remove "spatom"
del shells[ish]['spatom']
n_orbits=len(orbits)
#Note that the T numpy array is defined for all shells.
#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 = 'elk'
# Save it to the HDF:
ar = HDFArchive(self.hdf_file, 'a')
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]
del ar
# Save it to the HDF:
ar = HDFArchive(self.hdf_file, 'a')
symm_subgrp=self.symmcorr_subgrp
#Elk does not use time inversion symmetry
time_inv = [0 for j in range(n_symm)]
mat_tinv = [numpy.identity(orbits[orb]['dim'], complex)
for orb in range(n_orbits)]
#Save all the symmetry data
if not (symm_subgrp in ar):
ar.create_group(symm_subgrp)
things_to_save_sym = ['n_symm', 'n_atoms', 'perm',
'orbits', 'SO', 'SP', 'time_inv', 'mat', 'mat_tinv']
for it in things_to_save_sym:
ar[symm_subgrp][it] = locals()[it]
del ar
#Save misc info
things_to_save_misc = ['band_window','vkl','nstsv']
# Save it to the HDF:
ar = HDFArchive(self.hdf_file, 'a')
if not (self.misc_subgrp in ar):
ar.create_group(self.misc_subgrp)
for it in things_to_save_misc:
ar[self.misc_subgrp][it] = locals()[it]
del ar
mpi.report('Converted the Elk ground state data')
def convert_bands_input(self):
"""
Reads the appropriate files and stores the data for the bands_subgrp in the hdf5 archive.
"""
# Read and write only on the master node
if not (mpi.is_master_node()):
return
filext='_WANBAND.OUT'
dft_file='PROJ'+filext
mpi.report("Reading %s" % dft_file)
#Energy conversion - Elk uses Hartrees
energy_unit = 27.2113850560 # Elk uses hartrees
shells=[]
#read information about projectors calculated in the Elk calculation
[gen_info,n_corr_shells,n_inequiv_shells,corr_to_inequiv,inequiv_to_corr,corr_shells,n_reps,dim_reps,ind,basis,T] = read_Elk.read_proj(self,dft_file)
#get info for HDF5 file from gen_info
n_k=gen_info['n_k']
SP=gen_info['spinpol']-1
#Elk uses spinor wavefunctions. Therefore these two spinor wavefunctions have spin-orbit coupling incorporated in them. Here we read in the spinors
n_spin_blocs = SP + 1
SO=gen_info['SO']
#Elk only calculates Wannier projectors (no theta projectors generated):
n_shells=n_corr_shells
for ish in range(n_shells):
shells.append(corr_shells[ish].copy())
#remove last 2 entries from corr_shlls
del shells[ish]['SO']
del shells[ish]['irep']
shells[ish]['dim'] = 2*shells[ish]['l']+1
#read in the band eigenvalues
mpi.report("Reading BAND.OUT")
en=numpy.loadtxt('BAND.OUT')
nstsv=int(len(en[:,1])/n_k)
#convert the en array into a workable format
entmp = numpy.zeros([n_k,nstsv], complex)
enj=0
for ist in range(nstsv):
for ik in range(n_k):
entmp[ik,ist]=en[enj,1]
enj+=1
del en
#read projectors
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), nstsv], complex)
mpi.report("Reading projector(s)")
for ish in range(n_corr_shells):
[n_orbitals,band_window,rep,proj_mat]=read_Elk.read_projector(self,corr_shells,n_spin_blocs,ish,proj_mat,ind,T,basis,filext)
#alter arrays for spin-orbit coupling
if(SO==1):
mat=[]
su2=[]
n_symm=1
[shells,corr_shells,dim_reps,n_orbitals,proj_mat,T,mat]=self.update_so_quatities(n_shells,shells,n_corr_shells,corr_shells,n_inequiv_shells,dim_reps,n_k,n_symm,n_orbitals,proj_mat,T,su2,mat,sym=False)
#reduce n_spin_blocs
n_spin_blocs = SP + 1 - SO
#put the energy eigenvalues arrays in TRIQS format
hopping = self.sort_dft_eigvalues(n_spin_blocs,SO,n_k,n_orbitals,band_window,entmp,energy_unit)
# No partial projectors generated here, so set to 0:
n_parproj = numpy.array([0])
proj_mat_all = numpy.array([0])
# Save it to the HDF:
ar = HDFArchive(self.hdf_file, 'a')
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]
del ar
mpi.report('Converted the band data')
def convert_contours_input(self,kgrid=None,ngrid=None):
r"""
Reads the appropriate files and stores the data for the cont_subgrp in the hdf5 archive.
Parameters
----------
kgrid : size (4,3) double numpy array, optional
Numpy array defining the reciprocal lattice vertices used in the Elk Fermi
surface calculation. Each row has the following meaning:
grid3d[0,:] - origin lattice vertex
grid3d[1,:] - b1 lattice vertex
grid3d[2,:] - b2 lattice vertex
grid3d[3,:] - b3 lattice vertex
ngrid : size (3) integer numpy array, optional
Numpy array for the number of points along each (b1,b2,b3) lattice vertices
Note that these inputs relate to the plot3d input of Elk.
"""
if not (mpi.is_master_node()):
return
filext='_FS.OUT'
dft_file='PROJ'+filext
mpi.report("Reading %s" % dft_file)
#read the symmetries and k-points first
#read kpoints calculated in the Elk FS calculation
mpi.report("Reading KPOINT_FS.OUT")
[bz_weights,vkl]=read_Elk.read_kpoints(self,filext=filext)
n_k=vkl[:,0].size
#Need lattice symmetries to unfold the irreducible BZ
#Read symmetry files
mpi.report("Reading SYMCRYS.OUT")
[n_symm,spinmat,symlat,tr] = read_Elk.readsym(self)
#generate full vectors for Fermi surface plotting along with index mapping
#to irreducible vector set.
if (ngrid is not None and kgrid is not None):
mpi.report('Using User defined k-mesh')
#check variables are in correct format
if ngrid.size != 3:
assert 0, "The input numpy ngrid is not the required size of 3!"
elif ngrid.dtype != int:
assert 0, "The input numpy ngrid is not an array of integers."
elif kgrid.shape != (4,3):
assert 0, "The input numpy kgrid is not the required size of (4x3)!"
#generate full set of k-points with mapping to reduced set
[BZ_vkl, BZ_iknr, BZ_n_k] = Elk_tools.plotpt3d(self,n_k,vkl,n_symm,symlat,kgrid,ngrid)
elif (ngrid is None and kgrid is None):
mpi.report('No grid dimension input for Fermi surface.')
mpi.report('Calculating k-points by folding out irreducible vectors instead if using symmetries.')
mpi.report('Warning! This may not equate to the same set of vectors used to generate the Fermi surface data.')
[BZ_vkl, BZ_iknr, BZ_n_k] = Elk_tools.bzfoldout(self,n_k,vkl,n_symm,symlat)
else:
assert 0, "Either input both ngrid and kgrid numpy arrays or neither."
#return all threads apart from master
if not (mpi.is_master_node()):
return
# Read and write the following only on the master thread
#Energy conversion - Elk uses Hartrees
energy_unit = 27.2113850560 # Elk uses hartrees
shells=[]
#read information about projectors calculated in the Elk calculation
[gen_info,n_corr_shells,n_inequiv_shells,corr_to_inequiv,inequiv_to_corr,corr_shells,n_reps,dim_reps,ind,basis,T] = read_Elk.read_proj(self,dft_file)
#get info for HDF5 file from gen_info
n_k=gen_info['n_k']
SP=gen_info['spinpol']-1
#Elk uses spinor wavefunctions. Therefore these two spinor wavefunctions have spin-orbit coupling incorporated in them. Here we read in the spinors
n_spin_blocs = SP + 1
SO=gen_info['SO']
#Elk only calculates Wannier projectors (no theta projectors generated):
n_shells=n_corr_shells
for ish in range(n_shells):
shells.append(corr_shells[ish].copy())
#remove last 2 entries from corr_shlls
del shells[ish]['SO']
del shells[ish]['irep']
shells[ish]['dim'] = 2*shells[ish]['l']+1
#read in the eigenvalues used for the FS calculation
mpi.report("Reading EIGVAL_FS.OUT and EFERMI.OUT")
[en,occ,nstsv]=read_Elk.read_eig(self,filext=filext)
#read projectors
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), nstsv], complex)
mpi.report("Reading projector(s)")
for ish in range(n_corr_shells):
[n_orbitals,band_window,rep,proj_mat]=read_Elk.read_projector(self,corr_shells,n_spin_blocs,ish,proj_mat,ind,T,basis,filext)
mpi.report("Reading LATTICE.OUT")
[amat,amatinv,bmat,bmatinv,cell_vol] = read_Elk.readlat(self)
#Put eigenvalues into array of eigenvalues for the correlated window
#alter arrays for spin-orbit coupling
if(SO==1):
mat=[]
su2=[]
[shells,corr_shells,dim_reps,n_orbitals,proj_mat,T,mat]=self.update_so_quatities(n_shells,shells,n_corr_shells,corr_shells,n_inequiv_shells,dim_reps,n_k,n_symm,n_orbitals,proj_mat,T,su2,mat,sym=False)
#reduce n_spin_blocs
n_spin_blocs = SP + 1 - SO
#put the energy eigenvalues arrays in TRIQS format
hopping = self.sort_dft_eigvalues(n_spin_blocs,SO,n_k,n_orbitals,band_window,en,energy_unit)
# Save it to the HDF:
ar = HDFArchive(self.hdf_file, 'a')
if not (self.cont_subgrp in ar):
ar.create_group(self.cont_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','bmat',
'BZ_n_k','BZ_iknr','BZ_vkl','hopping']
for it in things_to_save:
ar[self.cont_subgrp][it] = locals()[it]
del ar
mpi.report('Converted the Contours data')
# commented out for now - unsure using this produces DFT+DMFT PDOS.
# The data from BC.OUT are the band-resolved diagonal muffin-tin DFT density matrix elements used in Elk to calculate PDOS
# (the PDOS is calculated from the Trace over the bands indices). Although this is equivalent to using using projectors in DFT and is likely valid for DFT+DMFT,
# the equivalence needs to be thoroughly checked for DFT+DMFT, but would require theta (or similar) projectors from Elk to do so.
# code left here just in case.
# def dft_band_characters(self):
# """
# Reads in the band-resolved muffin-tin density matrix (band characters) generated in Elk
# to be used for PDOS plots.
# """
# #determine file extension
# fileext='.OUT'
# #read number of k-points and eigenstates
# things_to_read = ['n_k','n_orbitals']
# ar = HDFArchive(self.hdf_file, 'r')
# for it in things_to_read:
# setattr(self, it, ar[self.dft_subgrp][it])
# del ar
# if not (mpi.is_master_node()):
# return
# mpi.report("Reading BC%s"%(fileext))
# # get needed data from hdf file
# # from general info
# ar = HDFArchive(self.hdf_file, 'a')
# things_to_read = ['SP', 'SO']
# for it in things_to_read:
# if not hasattr(self, it):
# setattr(self, it, ar[self.dft_subgrp][it])
# #from misc info
# things_to_read = ['nstsv','band_window']
# for it in things_to_read:
# if not hasattr(self, it):
# setattr(self, it, ar[self.misc_subgrp][it])
# #from sym info
# things_to_read = ['n_atoms','perm']
# symm_subgrp=self.symmcorr_subgrp
# for it in things_to_read:
# if not hasattr(self, it):
# setattr(self, it, ar[symm_subgrp][it])
# del ar
# #read in band characters
# [bc,maxlm] = read_Elk.read_bc(self,fileext)
#note that bc is the band resolved inner product of the Elk muffin-tin wave functions in
#a diagonal lm basis (by default). These are used in Elk to calculate the DOS in a diagonal
#irreducible lm basis. This band resolved density matrix in the lm basis will be used
#to project to spectral funtion to get the muffin-tin contributions. There will be an
#interstitial contribution (as Elk uses an APW+lo basis) which is the difference between
#the total and summed muffin-tin contributions. Also note that the bc array should be
#symmetrised within Elk.
#general variables
# lmax = int(numpy.sqrt(maxlm)-1)
# n_spin_blocs = self.SP + 1 - self.SO
# so = self.SO + 1
# #get the sort entry which is just the species index for Elk
# [ns, na, atpos]=read_Elk.read_geometry(self)
# isrt=0
# sort=numpy.zeros([self.n_atoms],int)
# #arrange sort(species) order
# for i in range(ns):
# for ia in range(na[i]):
# sort[isrt]=i
# isrt+=1
# #updating n_shells to include all the atoms and l used in the Elk calculation.
# n_shells = self.n_atoms * (lmax+1)
# shells = []
# shell_entries = ['atom', 'sort', 'l', 'dim']
# for iat in range(self.n_atoms):
# for l in range(lmax+1):
# #sort is not known from Elk outputs
# tmp = [iat+1, sort[iat]+1, l, so*(2*l+1)]
# shells.append({name: int(val) for name, val in zip(shell_entries, tmp)})
# del tmp, ns, na, atpos, isrt, shell_entries
# #overwrite n_shells and shells
# things_to_save = ['n_shells', 'shells']
# ar = HDFArchive(self.hdf_file, 'a')
# for it in things_to_save:
# ar[self.dft_subgrp][it] = locals()[it]
# # Initialise P, here a double list of matrices:
# band_dens_muffin = numpy.zeros([self.n_k, n_spin_blocs, n_shells, so*(2*lmax+1), numpy.max(self.n_orbitals)], float)
# for ik in range(self.n_k):
# for isp in range(n_spin_blocs):
# ish=0
# for iat in range(self.n_atoms):
# for l in range(lmax+1):
# #variables for putting subset of bc in proj_mat_all
# lm_min=l**2
# lm_max=(l+1)**2
# nst=self.n_orbitals[ik,isp]
# ibot=self.band_window[isp][ik, 0]-1
# itop=ibot+nst
# dim=l*2+1
#check use of abs (negative values should be close to 0)
# band_dens_muffin[ik,isp,ish,0:dim,0:nst] = \
# bc[lm_min:lm_max,isp,iat,ibot:itop,ik]
# if(self.SO==1):
# band_dens_muffin[ik,isp,ish,dim:2*dim,0:nst] = \
# bc[lm_min:lm_max,1,iat,ibot:itop,ik]
# ish+=1
# things_to_save = ['band_dens_muffin']
# # Save it all to the HDF:
# with HDFArchive(self.hdf_file, 'a') as ar:
# if not (self.bc_subgrp in ar):
# ar.create_group(self.bc_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_dens_muffin']
# for it in things_to_save:
# ar[self.bc_subgrp][it] = locals()[it]
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
# get needed data from hdf file
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)
things_to_read = ['SP', 'SO','n_k','n_orbitals']
for it in things_to_read:
if not hasattr(self, it):
setattr(self, it, ar[self.dft_subgrp][it])
#from misc info
things_to_read = ['band_window','vkl','nstsv']
for it in things_to_read:
if not hasattr(self, it):
setattr(self, it, ar[self.misc_subgrp][it])
#unlike in WIEN2k, Elk writes the velocities (momentum) matrix elements for all bands.
#Therefore, we can use the indices in the n_orbitals array to extract the desired elements.
#However, the PMAT.OUT file is in Fortran-binary, so the file is read in by python wrappers
#around the reading fortran code.
# Read relevant data from PMAT.OUT binary file
###########################################
# band_window_optics: same as Elk converter's band_window, but rearranged to be compatible
# for the transport calculations.
# velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics
# and each k-point.
#load fortran wrapper module
import triqs_dft_tools.converters.elktools.elkwrappers.getpmatelk as et
#elk velocities for all bands
pmat=numpy.zeros([self.nstsv,self.nstsv,3],dtype=complex)
n_spin_blocks = self.SP + 1 - self.SO
#TRIQS' velocities array used in its transport routines
velocities_k = [[] for isp in range(n_spin_blocks)]
#TRIQS' band_window array used in its transport routines
band_window_optics = []
mpi.report("Reading PMAT.OUT")
#read velocities for each k-point
for ik in range(self.n_k):
#need to use a fortran array for wrapper
f_vkl = numpy.asfortranarray(self.vkl[ik,:])
#read the ik velocity using the wrapper
pmat[:,:,:]=et.getpmatelk(ik+1,self.nstsv,f_vkl)
#loop over spin
for isp in range(n_spin_blocks):
#no. correlated bands at ik
nu1=self.band_window[isp][ik,0]-1
nu2=self.band_window[isp][ik,1]-1
n_bands=nu2-nu1+1
#put into velocity array (code similar to that in wien.py.
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)
#CHECK these lines
velocity_xyz[:,:,:]=pmat[nu1:nu2+1,nu1:nu2+1,:]
velocities_k[isp].append(velocity_xyz)
#rearrange Elk's band_window array into band_window_optics array format
for isp in range(n_spin_blocks):
band_window_optics_isp = []
for ik in range(self.n_k):
nu1=self.band_window[isp][ik,0]
nu2=self.band_window[isp][ik,1]
band_window_optics_isp.append((nu1, nu2))
n_bands=nu2-nu1+1
band_window_optics.append(numpy.array(band_window_optics_isp))
#read in the cell volume from LATTICE.OUT
mpi.report("Reading LATTICE.OUT")
[amat,amatinv,bmat,bmatinv,cell_vol] = read_Elk.readlat(self)
#read in the crystal symmetries
mpi.report("Reading SYMCRYS.OUT")
[n_symmetries,spinmat,rot_symmetries,tr] = read_Elk.readsym(self)
# 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]
things_to_save_misc = ['n_symmetries', 'rot_symmetries','cell_vol']
for it in things_to_save_misc:
ar[self.misc_subgrp][it] = locals()[it]
mpi.report("Reading complete!")
Reference in New Issue View Git Blame Copy Permalink