3
0
mirror of https://github.com/triqs/dft_tools synced 2024-12-23 12:55:17 +01:00
dft_tools/python/triqs_dft_tools/converters/elktools/readElkfiles.py
Xiao Jiang fafd625879 Update readElkfiles.py
During the Elk conversion process, two files, SYMCRYS.OUT and LATTICE.OUT, need to be read. However, after the Elk calculation is completed, there are some blank lines in these files, which leads to recognition failure of readElkfiles.py.  By using line.strip(), non-blank lines can be directly read.
2024-09-09 20:02:13 +02:00

690 lines
26 KiB
Python

##########################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2019 by A. D. N. James, M. Zingl 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/>.
#
##########################################################################
from types import *
import numpy
from triqs_dft_tools.converters.converter_tools import *
import os.path
from locale import atof
class readElkfiles:
"""
Read in the Elk output files for the elk_converter.py routines.
"""
def __init__(self):
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'
def read_elk_file(self, filename, to_replace):
#this function is almost identical to read_fortran_file, but it removes words after ':'
"""
Returns a generator that yields all numbers in the Fortran file as float, with possible replacements.
Parameters
----------
filename : string
Name of Fortran-produced file.
to_replace : dict of str:str
Dictionary defining old_char:new_char.
Returns
-------
string
The next number in file.
"""
#import os.path
#import string
import locale
from locale import atof
if not(os.path.exists(filename)):
raise IOError("File %s does not exist." % filename)
for line in open(filename, 'r'):
for old, new in to_replace.items():
line = line.replace(old, new)
# removes the substring after ':' character - if it exists
if ":" in line:
line = self.split_string(line)
if "(" in line:
line = self.split_string3(line)
for x in line.split():
yield atof(x)
def read_elk_file2(self, filename, to_replace):
#this function is almost identical to read_fortran_file, but it removes words after '('
"""
Returns a generator that yields all numbers in the Fortran file as float, with possible replacements.
Parameters
----------
filename : string
Name of Fortran-produced file.
to_replace : dict of str:str
Dictionary defining old_char:new_char.
Returns
-------
string
The next number in file.
"""
import os.path
import string
if not(os.path.exists(filename)):
raise IOError("File %s does not exist." % filename)
for line in open(filename, 'r'):
if line.strip():
for old, new in to_replace.items():
line = line.replace(old, new)
# removes the substring after '(' character - if it exists
if "(" in line:
line = self.split_string2(line,'(')
# removes the substring after '+' character - if it exists
if "+" in line:
line = self.split_string2(line,'+')
# removes the substring after '|' character - if it exists
if "|" in line:
line = self.split_string2(line,'|')
#only include lines which have multiple characters
if(len(line)>1):
yield line.split()
def split_string(self,line):
"""
This removes the excess information after ':' of the string read in from the file
"""
temp = line.split(':')
line = temp[0]
return line
def split_string2(self,line,str):
"""
This removes the excess information after 'str' of the string read in from the file
"""
temp = line.split(str)
line = temp[0]
return line
def split_string3(self,line):
"""
This removes the excess information after '(' of the string read in from the file
"""
temp = line.split('(')
line = temp[0]
return line
def read_proj(self,dft_file):#,things_to_set):
"""
This function reads the contents of PROJ.OUT and returns them for general use.
"""
#read in PROJ.OUT file
R = self.read_elk_file( dft_file, self.fortran_to_replace)
try:
#arrays of entry names
gen_info_entries = ['nproj', 'n_k', 'spinpol', 'SO','natm']
proj_entries = ['sort', 'natom', 'l', 'dim']
at_entries = ['spatom','atom']
#Read in the No. of projectors, no. of k-points, spin and spinorb info
gen_info = {name: int(val) for name, val in zip(gen_info_entries, R)}
#mpi.report(gen_info)
#read in the information of all the projectors
#initialized variables
icorr=0
lmax=3 #maximum l value that will be used
n_shells=0
n_inequiv_shells=0
#local arrays
neqatom=[]#numpy.zeros([n_shells], int)
proj=[]
shells=[]#numpy.zeros([n_shells], int)
corr_shells=[]#numpy.zeros([n_shells], int)
prjtype=[]
wan=[]
proj_info=[]#numpy.zeros([n_shells], int)
T=[]
basis=[]
inequiv_to_corr=[]
corr_to_inequiv=[]
ind=[]
proj_idx=[]
at=[]
idxlm=[]
#assign global variables
n_k=gen_info['n_k']
SO=gen_info['SO']
#loop over projector flags
for ip in range(gen_info['nproj']):
#projector index
proj_idx.append(int(next(R)))
proj.append({name: int(val) for name, val in zip(proj_entries, R)})
na=proj[ip]['natom'] # integer which reduces when reading in atom info
while(na>0):
neqatom.append(int(next(R)))
na-=neqatom[n_inequiv_shells]
#appends the mapping array to index of inequivalent atom in shells
inequiv_to_corr.append(icorr)
for ia in range(neqatom[n_inequiv_shells]):
corr_to_inequiv.append(n_inequiv_shells)
at.append({name: int(val) for name, val in zip(at_entries, R)})
shells.append(proj[ip].copy())
shells[n_shells].update({'atom':at[n_shells]['atom'], 'spatom':at[n_shells]['spatom']})
corr_shells.append(shells[n_shells].copy())
n_orb=2*shells[n_shells]['l']+1
#lm submatrix indices
idxlm.append(numpy.zeros(2*lmax+1, dtype=int))
nrep=proj[ip]['dim']
for i in range(nrep):
idxlm[n_shells][i]=next(R)-1
ind.append(idxlm[n_shells][0:nrep])
#determine basis type and transformation matrix
basis.append(int(next(R)))
#determine whether which basis the projectors where generated in
#spherical harmonics
T.append(numpy.zeros([n_orb, n_orb], dtype=complex))
#Elk generated unitary basis
if (basis[n_shells]==2):
#reads the transformation matrix
for i in range(n_orb):
for j in range(n_orb):
T[n_shells][i, j] = next(R)
for i in range(n_orb):
for j in range(n_orb):
T[n_shells][i, j] += 1j * next(R)
#transformation matrix from spherical harmonics to cubic basis
elif (basis[n_shells]==1):
T[n_shells]=self.determine_T(shells[n_shells]['l'])
else:
for i in range(n_orb):
T[n_shells][i,i] = 1.0
#index for the next inequivalent atom (+1 to index is not needed as this is incorporated in
#neqatom[ish]
n_shells+=1
#increase the inequiv_to_corr value
icorr+=neqatom[n_inequiv_shells]
#increase the numer of inequivalent atoms
n_inequiv_shells+=1
#end reading the file if read the last projector
if (ip+1==gen_info['nproj']):
break
#determine the irreducible representation
dim_reps=[]
n_reps=[]
irep=[]
for ish in range(n_inequiv_shells):
isheq=inequiv_to_corr[ish]
[n_reps,dim_reps,irep] = self.determine_rep(isheq,ish,corr_shells,basis,ind,n_reps,dim_reps,irep)
for ish in range(n_shells):
ishin=corr_to_inequiv[ish]
corr_shells[ish].update({'SO':SO, 'irep':irep[ishin]})
except StopIteration: # a more explicit error if the file is corrupted.
#raise "Elk_converter : reading PROJ.OUT file failed!"
raise IOError("Elk_converter : reading PROJ.OUT file failed!")
R.close()
#output desired information
return (gen_info,n_shells,n_inequiv_shells,corr_to_inequiv,inequiv_to_corr,corr_shells,n_reps,dim_reps,ind,basis,T)
def determine_T(self,l):
"""
Current version calculates the transformation matrix T to convert the inputs from spherical
harmonics to the cubic basis (as used in TRIQS and Wien2k).
This routine is currently very similar to spherical_to_cubic in the TRIQS library.
This routine can be extended to include other unitary rotation matrices
"""
from math import sqrt
size=2*l+1
T = numpy.zeros((size,size),dtype=complex)
if(l == 0):
T[0,0]= 1.0
elif(l == 1):
T[0,0] = 1.0/sqrt(2); T[0,2] = -1.0/sqrt(2)
T[1,0] = 1j/sqrt(2); T[1,2] = 1j/sqrt(2)
T[2,1] = 1.0
elif l == 2:
T[0,2] = 1.0
T[1,0] = 1.0/sqrt(2); T[1,4] = 1.0/sqrt(2)
T[2,0] =-1.0/sqrt(2); T[2,4] = 1.0/sqrt(2)
T[3,1] = 1.0/sqrt(2); T[3,3] =-1.0/sqrt(2)
T[4,1] = 1.0/sqrt(2); T[4,3] = 1.0/sqrt(2)
elif l == 3:
#This needs to be checked
T[0,0] = 1.0/sqrt(2); T[0,6] = -1.0/sqrt(2)
T[1,1] = 1.0/sqrt(2); T[1,5] = 1.0/sqrt(2)
T[2,2] = 1.0/sqrt(2); T[2,4] = -1.0/sqrt(2)
T[3,3] = 1.0
T[4,2] = 1j/sqrt(2); T[4,4] = 1j/sqrt(2)
T[5,1] = 1j/sqrt(2); T[5,5] = -1j/sqrt(2)
T[6,0] = 1j/sqrt(2); T[6,6] = 1j/sqrt(2)
else: raise ValueError("determine: implemented only for l=0,1,2,3")
return numpy.matrix(T)
def determine_rep(self,ish,ishin,corr_shells,basis,ind,n_reps,dim_reps,irep):
"""
Determines the irreducible representation used for projection calculation.
Only for Cubic harmonics at the moment
"""
#if all of the orbitals were used to construct the projectors
rep=corr_shells[ish]['dim']
if(rep==2*corr_shells[ish]['l']+1):
n_reps.append(1)
dim_reps.append([corr_shells[ish]['dim']])
irep.append(1)
#if a subset of obitals were used for generating the projectors (l>1)
#if Elk generated T matrix is used
elif((basis[ish]==2)&(corr_shells[ish]['l']>1)):
if (corr_shells[ish]['l']==2):
n_reps.append(2)
dim_reps.append([3, 2])
irep.append(dim_reps[ishin].index(rep)+1)
elif (corr_shells[ish]['l']==3):
n_reps.append(3)
dim_reps.append([1, 3, 3])
#cubic T matrix
elif((basis[ish]==1)&(corr_shells[ish]['l']>1)):
if (corr_shells[ish]['l']==2):
n_reps.append(2)
dim_reps.append([2, 3])
irep.append(dim_reps[ishin].index(rep)+1)
elif (corr_shells[ish]['l']==3):
n_reps.append(3)
dim_reps.append([2, 2, 3])
#determine the dim_reps from the lm indices in ind
if(rep==3):
irep.append(3)
else:
for i in range(2):
#look for ind[i]==0
if(ind[i]==0):
irep.append(1)
break
#default to irep=2
elif(i==1):
irep.append(2)
else:
raise ValueError("Elk_converter (determine_rep) : irreducible representations were not found!")
return (n_reps,dim_reps,irep)
def read_projector(self,shell,n_spin_blocks,ish,proj_mat,ind,T,basis,filext):
"""
This function reads the contents of WANPROJ_L**_S**_A****.OUT
"""
import string
l = str(shell[ish]['l']).zfill(2)
s = str(shell[ish]['sort']).zfill(2)
a = str(shell[ish]['spatom']).zfill(4)
#construct file name from atoms info
wan_shells_file = 'WANPROJ_L'+l+'_S'+s+'_A'+a+filext
R = self.read_elk_file(wan_shells_file, self.fortran_to_replace)
try:
#no. of kpts and number of orbital
gen_entries = ['n_k', 'lmmax','dim_rep']
gen = {name: int(val) for name, val in zip(gen_entries, R)}
#projector lm size
dim=gen['lmmax']
#no. of lm used to construct projector
dim_rep=gen['dim_rep']
lat=[]
n_k=gen['n_k']
n_orbitals = numpy.zeros([n_k, n_spin_blocks], int)
band_window = [None for isp in range(n_spin_blocks)]
for isp in range(n_spin_blocks):
band_window[isp] = numpy.zeros([n_k, 2], dtype=int)
for ik in range(0,n_k):
#k-point index and correlated band window indices
lat_entries = ['ik','vklx','vkly','vklz']
lat.append({name: int(val) for name, val in zip(lat_entries, R)})
#loop over each spin
for isp in range(0,n_spin_blocks):
#band window indices
proj_entries = ['isp', 'ist_min','ist_max']
proj_dim={name: int(val) for name, val in zip(proj_entries, R)}
if(isp+1!=proj_dim['isp']):
raise IOError("Elk_converter (",wan_shells_file,") : reading spin projecotrs failed!")
return
#setting band window arrays
n_orbitals[ik,isp]=proj_dim['ist_max']-proj_dim['ist_min']+1
band_window[isp][ik, 0] = proj_dim['ist_min']
band_window[isp][ik, 1] = proj_dim['ist_max']
#define temporary matrix for reading in the projectors
mat = numpy.zeros([dim, n_orbitals[ik,isp]], complex)
# Real part
for j in range(dim):
for i in range(n_orbitals[ik,isp]):
mat[j, i] = next(R)
# Imag part:
for j in range(dim):
for i in range(n_orbitals[ik,isp]):
mat[j, i] += 1j * next(R)
#rotate projectors into basis T
if(basis[ish]!=0):
mat[:,:]=numpy.matmul(T[ish],mat[:,:])
#put desired projector subset into master array
proj_mat[ik,isp,ish,0:dim_rep,0:n_orbitals[ik,isp]]=mat[ind[ish],0:n_orbitals[ik,isp]]
#delete temporary index
del mat
#resize array length of band indices to maximum number used
proj_mat=proj_mat[:,:,:,:,:numpy.max(n_orbitals)]
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter (read projectors): reading file failed!")
R.close()
return(n_orbitals,band_window,dim_rep,proj_mat)
def read_eig(self,filext=".OUT"):
"""
This function reads the contents of EIGVAL.OUT and EFERMI.OUT
"""
import string
#construct file name from atoms info
R = self.read_elk_file( self.efermi_file, self.fortran_to_replace)
try:
efermi=next(R)
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter (EFERMI.OUT): reading file failed!")
eval_file = 'EIGVAL'+filext
R = self.read_elk_file( eval_file, self.fortran_to_replace)
try:
#no. of kpts
n_k=int(next(R))
#no. second variation states (no. bands)
nstsv=int(next(R))
en=[]
occ=[]
kp2=[]
for ik in range(0,n_k):
#k index and corresponing lattice coordinates
k_entries = ['ik', 'vklx','vkly','vklz']
kp2.append([{name: val for name, val in zip(k_entries, R)}])
en.append([])
occ.append([])
for ist in range(0,nstsv):
#stores the band index, energy eigenvalues and occupancies
en_entries = ['ist', 'hopping','occ']
read_en={name: val for name, val in zip(en_entries, R)}
#Energy eigenvalues and occupations
en[ik].append(read_en['hopping']-efermi)
occ[ik].append(read_en['occ'])
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter (EIGVAL.OUT): reading file failed!")
R.close()
return(en,occ,nstsv)
def read_kpoints(self,filext=".OUT"):
"""
This function reads the contents of KPOINTS.OUT
"""
import string
kp_file = 'KPOINTS'+filext
R = self.read_elk_file( kp_file, self.fortran_to_replace)
try:
#no. of kpts
n_k=int(next(R))
kp=[]
#reads in the k index, lattice vectors, weights and nmat for each kpt
#array for bz weights
bz_weights = numpy.ones([n_k], float) / float(n_k)
#array for lattice vectors
vkl = numpy.ones([n_k,3], float)
for ik in range(n_k):
#k-grid info
k_entries = ['ik', 'vklx','vkly','vklz', 'bz_weights', 'nmat']
kp.append({name: val for name, val in zip(k_entries, R)})
#convert to integers
kp[ik]['ik']=int(kp[ik]['ik'])
#read in BZ weights
bz_weights[ik]=kp[ik]['bz_weights']
j=0
for i in range(1,4):
vkl[ik,j]=kp[ik][k_entries[i]]
j=j+1
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter (KPOINTS.OUT): reading file failed!")
R.close()
return(bz_weights,vkl)
def readsym(self):
"""
Read in the (crystal) symmetries in lattice coordinates
"""
dft_file='SYMCRYS.OUT'
R = self.read_elk_file2( dft_file, self.fortran_to_replace)
try:
symmat=[]
spinmat=[]
tr=[]
#read the number of crystal symmetries
x = next(R)
nsym = int(atof(x[0]))
#set up symmetry matrices
for isym in range(nsym):
symmat.append(numpy.zeros([3, 3], float))
spinmat.append(numpy.zeros([3, 3], float))
tr.append(numpy.zeros([3], float))
#read each symmetry
for isym in range(nsym):
#read the symmetry index and check it
x = next(R)
isymm = int(x[3])
if(isym+1!=isymm):
raise IOError("Elk_converter : reading symmetries failed!")
#next read has no useful infomation
x = next(R)
#read in the translation vector used for symmetries
x = next(R)
for i in range(3):
tr[isym][i]=atof(x[i])
#read in the spatial symmetries
#string with no useful information
x = next(R)
#read in the spatial symmetry
for i in range(3):
x = next(R)
for j in range(3):
symmat[isym][i,j]=int(atof(x[j]))
#string with no useful information
x = next(R)
#read in the spin symmetries
for i in range(3):
x = next(R)
for j in range(3):
spinmat[isym][i,j]=int(atof(x[j]))
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter : reading SYMCRYS.OUT file failed!")
R.close()
return nsym,spinmat,symmat,tr
def readlat(self):
"""
Read in information about the lattice.
"""
dft_file='LATTICE.OUT'
R = self.read_elk_file2( dft_file, self.fortran_to_replace)
try:
amat = numpy.zeros([3, 3], float)
amatinv = numpy.zeros([3, 3], float)
bmat = numpy.zeros([3, 3], float)
bmatinv = numpy.zeros([3, 3], float)
#real space lattice matrices
#cycling through information which is not needed
for i in range(4):
x = next(R)
#reading in the lattice vectors as matrix
for i in range(3):
x = next(R)
for j in range(3):
amat[i,j] = atof(x[j])
#cycling through information which is not needed
x = next(R)
#reading in the inverse lattice matrix
for i in range(3):
x = next(R)
for j in range(3):
amatinv[i,j] = atof(x[j])
#read in cell volume (for transport)
x = next(R)
cell_vol = atof(x[-1])
#cycling through information which is not needed
for i in range(4):
x = next(R)
#reading in the reciprocal lattice vectors as matrix
for i in range(3):
x = next(R)
for j in range(3):
bmat[i,j] = atof(x[j])
#cycling through information which is not needed
x = next(R)
#reading in the inverse reciprocal lattice matrix
for i in range(3):
x = next(R)
for j in range(3):
bmatinv[i,j] = atof(x[j])
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter : reading PROJ.OUT file failed!")
R.close()
return amat,amatinv,bmat,bmatinv,cell_vol
def read_geometry(self):
"""
This function reads the contents of GEOMETRY.OUT
"""
import string
#read in the Geometery file
dft_file='GEOMETRY.OUT'
#atom positions in lattice coordinates
R = self.read_elk_file2( dft_file, self.fortran_to_replace)
nspecies=0
na=[]
atpos_entries = ['kx', 'ky','kz','Bkx', 'Bky', 'Bkz']
atpos=[]
try:
i=0
#cycle through file until "atoms" is found
while(i<1000):
x = next(R)
if "atoms" in x:
break
elif(i==1000):
raise IOError("Elk_converter : Could not find 'atoms' in GEOMETRY.OUT!")
i+=1
#read in the number of species
x = next(R)
ns = int(atof(x[0]))
#loop over species
for js in range(ns):
#Skip species name
x = next(R)
#read in the number of atoms per species
x = next(R)
na.append(int(atof(x[0])))
#loop over atomss pre species
atpos.append([])
for ia in range(na[js]):
atpos[js].append(numpy.zeros(6, float))
x = next(R)
for j in range(6):
atpos[js][ia][j]=atof(x[j])
except StopIteration: # a more explicit error if the file is corrupted.
raise IOError("Elk_converter : Could not open GEOMETRY.OUT!")
R.close()
return ns, na, atpos
#commented out for now - unsure this will produce DFT+DMFT PDOS
##band character dependent calculations
# def read_bc(self,fileext):
# """
# Read in the ELK generated band characters from BC.OUT
# """
# #import string
# file = 'BC'+fileext
# R = self.read_elk_file(file, self.fortran_to_replace)
# try:
# #no. of kpts and number of orbital
# gen_entries = ['maxlm', 'nspinor','natmtot','nstsv','nkpt','irep']
# gen = {name: int(val) for name, val in zip(gen_entries, R)}
# #projector lm size
# #check the read in information complies with previous read in data
# nspinor=self.SP+1
# if(gen['nspinor'] != nspinor):
# mpi.report('HDF file nspinor = %s'%nspinor)
# mpi.report('BC.OUT nspinor = %s'%gen['nspinor'])
# raise IOError("Elk_converter (",file,") : reading nspinor failed!")
# return
# if(gen['natmtot'] != self.n_atoms):
# raise IOError("Elk_converter (",file,") : reading no. of atoms failed!")
# return
# if(gen['nstsv'] != self.nstsv):
# raise IOError("Elk_converter (",file,") : reading all states failed!")
# return
# if(gen['nkpt'] != self.n_k):
# raise IOError("Elk_converter (",file,") : reading kpoints failed failed!")
# return
# if(gen['irep'] == 0):
# raise IOError("Elk_converter (",file,") : Band characters are in spherical hamonics, may have issues with the PDOS!")
# return
# dim=gen['maxlm']
# lmax=numpy.sqrt(dim)-1
# bc = numpy.zeros([dim,nspinor,self.n_atoms,self.nstsv,self.n_k], float)
# for ik in range(0,self.n_k):
# for iatom in range(0,self.n_atoms):
# for ispn in range(0,nspinor):
# entry = ['ispn','ias','is','ia','ik']
# ent = {name: int(val) for name, val in zip(entry, R)}
# #k-point index and correlated band window indices
# #check read in values
# if(ent['ispn'] != ispn+1):
# raise IOError("Elk_converter (",file,") : reading ispn failed!")
# return
# if(ent['ias'] != iatom+1):
# raise IOError("Elk_converter (",file,") : reading iatom failed!")
# return
# if(ent['ik'] != ik+1):
# raise IOError("Elk_converter (",file,") : reading ik failed!")
# return
# for ist in range(self.nstsv):
# for lm in range(dim):
# bc[lm,ispn,iatom,ist,ik] = next(R)
# except StopIteration: # a more explicit error if the file is corrupted.
# raise IOError("Elk_converter (read BC.OUT): reading file failed!")
# R.close()
# return(bc,dim)