mirror of
https://github.com/triqs/dft_tools
synced 2024-12-23 12:55:17 +01:00
ad8c4e75fe
Adding Elk-TRIQS interface (first iteration) This interface reads in Elk's ground-state files / projectors generated by a specific Elk interface code version (https://github.com/AlynJ/Elk_interface-TRIQS). The interface can perform charge-self consistent DFT+DMFT calculations using the aforementioned Elk code version, including spin orbit-coupling. Hence, this is the first open source interfaced DFT code to triqs with FCSC support. The commit includes detailed documentation and tutorials on how to use this interface. Moreover, further new post-processing routines are added for Fermi surface plots and spectral functions (A(w)) from the elk inputs. The interface was tested by A. James and A. Hampel. However, this is the first iteration of the interface and should be used with care. Please report all bugs. List of changes: --------------- - sumk.py: added cacluation of charge density correction for elk (dm_type='elk'). - sumk_dft_tools.py: added new post-processing functions to calculate the Fermi surface and A(w) from the Elk inputs. - documentation and tutorial files amended for this interface. - added python tests for the Elk converter.
688 lines
26 KiB
Python
688 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.
|
|
|
|
Yields
|
|
------
|
|
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.
|
|
|
|
Yields
|
|
------
|
|
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'):
|
|
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], numpy.int)
|
|
proj=[]
|
|
shells=[]#numpy.zeros([n_shells], numpy.int)
|
|
corr_shells=[]#numpy.zeros([n_shells], numpy.int)
|
|
prjtype=[]
|
|
wan=[]
|
|
proj_info=[]#numpy.zeros([n_shells], numpy.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[ip]
|
|
#appends the mapping array to index of inequivalent atom in shells
|
|
inequiv_to_corr.append(icorr)
|
|
for ia in range(neqatom[ip]):
|
|
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=numpy.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=numpy.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 numer of inequivalent atoms
|
|
n_inequiv_shells+=1
|
|
#increase the inequiv_to_corr value
|
|
icorr+=neqatom[ip]
|
|
#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], numpy.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]], numpy.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], numpy.float_) / float(n_k)
|
|
#array for lattice vectors
|
|
vkl = numpy.ones([n_k,3], numpy.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=[]
|
|
#maximum symmetries
|
|
nsym = 48
|
|
#set up symmetry matrices
|
|
for isym in range(nsym):
|
|
symmat.append(numpy.zeros([3, 3], numpy.float_))
|
|
spinmat.append(numpy.zeros([3, 3], numpy.float_))
|
|
tr.append(numpy.zeros([3], numpy.float_))
|
|
#read the number of crystal symmetries
|
|
x = next(R)
|
|
nsym = int(atof(x[0]))
|
|
#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 the symmetries in lattice coordinates
|
|
"""
|
|
dft_file='LATTICE.OUT'
|
|
R = self.read_elk_file2( dft_file, self.fortran_to_replace)
|
|
try:
|
|
amat = numpy.zeros([3, 3], numpy.float_)
|
|
amatinv = numpy.zeros([3, 3], numpy.float_)
|
|
bmat = numpy.zeros([3, 3], numpy.float_)
|
|
bmatinv = numpy.zeros([3, 3], numpy.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])
|
|
#reciprocal lattice matrices
|
|
#cycling through information which is not needed
|
|
for i in range(5):
|
|
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
|
|
|
|
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, numpy.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
|
|
|
|
#band character dependent calculations
|
|
def read_bc(self):
|
|
"""
|
|
Read in the ELK generated band characters from BC.OUT
|
|
"""
|
|
|
|
#import string
|
|
file = 'BC.OUT'
|
|
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], numpy.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)
|
|
|