#!/usr/bin/env python
import datetime, os, re, sys
import numpy as np
np.seterr( all='raise', under='ignore')
import readVasp
import wrapUpload
#====================================================================
def badparms(
msg): # error message
'''
Displays an error message and quits.
Command line parameters for the test driver:
============== ========= ==============================================
Parameter Type Description
============== ========= ==============================================
-bugLev integer Debug level. Normally 0.
-inDir string Input dir containing OUTCAR, INCAR, POSCAR
============== ========= ==============================================
'''
print '\nError: %s' % (msg,)
print '\n%s' % (badparms.__doc__,)
sys.exit(1)
#====================================================================
[docs]def main():
'''
Test driver.
Extracts info from VASP OUTCAR, INCAR, and POSCAR files.
See doc for badparms.
'''
bugLev = 0
inDir = None
if len(sys.argv) % 2 != 1:
badparms('Parms must be key/value pairs')
for iarg in range( 1, len(sys.argv), 2):
key = sys.argv[iarg]
val = sys.argv[iarg+1]
if key == '-bugLev': bugLev = int( val)
elif key == '-inDir': inDir = val
else: badparms('unknown key: "%s"' % (key,))
if bugLev == None: badparms('parm not specified: -bugLev')
if inDir == None: badparms('parm not specified: -inDir')
resObj = ResClass()
scannerUnused = ScanOutcar( bugLev, inDir, resObj)
#====================================================================
#====================================================================
[docs]class ResClass:
'''
A container for the results of scanOutcar.
The results are saved as attributes of the object.
'''
def __str__(
self): # self
'''
Formats self.
'''
keys = self.__dict__.keys()
keys.sort()
msg = ''
for key in keys:
val = self.__dict__[key]
stg = str( val)
if stg.find('\n') >= 0: sep = '\n'
else: sep = ' '
msg += ' %s: type: %s val:%s%s\n' % (key, type(val), sep, val,)
return msg
#====================================================================
#====================================================================
[docs]class Sspec:
'''
Specifies regex, numMin, etc, for one scalar to be found in OUTCAR.
'''
def __init__(
self, # self
tag, # resObj attribute name to set
pat, # regex to match in OUTCAR
which, # 'first' or 'last': which match to use
numMin, # min number of matches
numMax, # max number of matches, or 0 for any number
tp): # result type: int, float, or str
'''
Saves the parameters.
'''
self.tag = tag
self.pat = pat
self.which = which
self.numMin = numMin
self.numMax = numMax
self.tp = tp
if not (self.pat.startswith('^') and self.pat.endswith('$')):
self.throwerr('invalid spec: %s' % (self,))
if which not in ['first', 'last']:
self.throwerr('invalid spec: %s' % (self,))
if tp not in [ int, float, str]:
self.throwerr('invalid spec: %s' % (self,))
self.regex = re.compile( pat)
def __str__(
self): # self
'''
Formats self.
'''
res = 'tag: %s pat: %s which: %s numMin: %d numMax: %d tp: %s' \
% (self.tag, repr( self.pat), self.which,
self.numMin, self.numMax, self.tp,)
return res
def throwerr(
self, # self
msg): # error message
'''
Raises an Exception.
'''
raise Exception( msg)
#====================================================================
#====================================================================
# Fills resObj.
[docs]class ScanOutcar:
'''
Extracts info from OUTCAR, INCAR, and POSCAR files.
Typical usage::
resObj = ResClass()
ScanOutcar( bugLev, inDir, resObj)
print 'ediff: %g' % (resObj.ediff,)
'''
def __init__(
self, # self
bugLev, # debug level
inDir, # input directory
resObj): # resulting object to be filled
'''
Sets variables.
'''
self.bugLev = bugLev
self.inDir = inDir
if bugLev >= 1:
print 'ScanOutcar: inDir: %s' % (inDir,)
# First, read POSCAR and INCAR
if os.path.exists( os.path.join( self.inDir, 'INCAR')):
self.parseIncar( resObj)
# We no longer use the POSCAR file since everything
# in it can be derived from the OUTCAR file.
#if os.path.exists( os.path.join( self.inDir, 'POSCAR')):
# self.parsePoscar( resObj)
# Now read OUTCAR
fname = os.path.join( inDir, 'OUTCAR')
with open( fname) as fin:
self.lines = fin.readlines()
self.numLine = len( self.lines)
# Strip all lines
for ii in range( self.numLine):
self.lines[ii] = self.lines[ii].strip()
self.getScalars( resObj)
self.getDate( resObj)
self.getTimes( resObj)
self.getSystem( resObj)
self.getCategEnergy( resObj)
self.getBasisMats( resObj)
self.getTypeNames( resObj)
self.getTypeNums( resObj)
self.getTypeMassValence( resObj)
self.getTypePseudos( resObj)
numType = len( resObj.typeNames)
self.getLdau( numType, resObj)
if bugLev >= 5:
print '\nunsorted atomTypes:'
print 'numType: %d' % ( numType,)
print 'typeNames: %s' % ( resObj.typeNames,)
print 'typeNums: %s' % ( resObj.typeNums,)
print 'typeMasses_amu: %s' % ( resObj.typeMasses_amu,)
print 'typeValences: %s' % ( resObj.typeValences,)
print 'typePseudos: %s' % ( resObj.typePseudos,)
print 'numAtom: %d' % ( resObj.numAtom,)
if len( resObj.typeNums) != numType \
or len( resObj.typeMasses_amu) != numType \
or len( resObj.typeValences) != numType \
or len( resObj.typePseudos) != numType:
errMsg = 'types mismatch. fname: %s\n typeNames: %s\n typeNums: %s\n typeMasses_amu: %s\n typeValences: %s\n typePseudos: %s\n' \
% ( fname, resObj.typeNames, resObj.typeNums, resObj.typeMasses_amu, resObj.typeValences, resObj.typePseudos,)
print '\n%s\n' % (errMsg,)
# self.throwerr( errMsg, None)
## Sort parallel arrays typeNames, typeNums, etc,
## by typeNames alphabetic order.
## Use an index sort with typeIxs.
## In rare cases like icsd_024360.cif/hs-ferro
## and icsd_060845_GW, the names are out of order.
## Sort to set typeIxs[newIx] == oldIx.
#ntype = len( resObj.typeNames)
#typeIxs = range( ntype)
#def sortFunc(
# ia, # resObj A
# ib): # resObj B
# '''
# Sort function for sorting resObj.typeNames
# '''
# return cmp( resObj.typeNames[ia], resObj.typeNames[ib])
#typeIxs.sort( sortFunc)
#if bugLev >= 5:
# print 'typeIxs: %s' % (typeIxs,)
#resObj.typeIxs = typeIxs
#resObj.typeNames = [resObj.typeNames[ix] for ix in typeIxs]
#resObj.typeNums = [resObj.typeNums[ix] for ix in typeIxs]
#resObj.typeMasses_amu = [resObj.typeMasses_amu[ix] for ix in typeIxs]
#resObj.typeValences = [resObj.typeValences[ix] for ix in typeIxs]
#resObj.typePseudos = [resObj.typePseudos[ix] for ix in typeIxs]
#if bugLev >= 5:
# print '\nsorted atomTypes:'
# print 'typeIxs: %s' % ( resObj.typeIxs,)
# print 'typeNames: %s' % ( resObj.typeNames,)
# print 'typeNums: %s' % ( resObj.typeNums,)
# print 'typeMasses_amu: %s' % ( resObj.typeMasses_amu,)
# print 'typeValences: %s' % ( resObj.typeValences,)
# print 'typePseudos: %s' % ( resObj.typePseudos,)
## Make sure typenames are in alphabetic order
#for ii in range(len(resObj.typeNames) - 1):
# if resObj.typeNames[ii] > resObj.typeNames[ii+1]:
# throwerr('typeNames not in order')
## Make sure atomnames are in alphabetic order
#for ii in range(len(resObj.atomNames) - 1):
# if resObj.atomNames[ii] > resObj.atomNames[ii+1]:
# throwerr('atomNames not in order')
# New way to get atomNames, atomMasses, atomValences, atomPseudos
resObj.atomNames = []
resObj.atomMasses_amu = []
resObj.atomValences = []
resObj.atomPseudos = []
for ii in range( numType):
nn = resObj.typeNums[ ii]
resObj.atomNames += nn * [resObj.typeNames[ii]]
resObj.atomMasses_amu += nn * [resObj.typeMasses_amu[ii]]
resObj.atomValences += nn * [resObj.typeValences[ii]]
resObj.atomPseudos += nn * [resObj.typePseudos[ii]]
numAtom = len( resObj.atomNames)
self.getInitialPositions( resObj)
self.getFinalPositionsForce( resObj)
self.getStressMat( resObj)
self.getMagmoms( resObj)
self.getKpointMat( resObj)
self.getEigenMat( resObj)
self.getDielectric( resObj)
# Check that the POSCAR posScaleFactor agrees
# with the calculation we did in getKpointMat using OUTCAR.
#
#fname = os.path.join( self.inDir, 'POSCAR')
#with open( fname) as fin:
# lines = fin.readlines()
#scalefac = float( lines[1].strip())
#print 'check scaleFactor: poscar: %g outcar: %g out-pos: %g' \
# % (scalefac, resObj.posScaleFactor, resObj.posScaleFactor - scalefac,)
#====================================================================
[docs] def parseIncar(
self, # self
resObj): # resulting object to be filled
'''
Extracts info from a VASP INCAR file, setting resObj attributes.
'''
fname = os.path.join( self.inDir, 'INCAR')
with open( fname) as fin:
lines = fin.readlines()
for iline in range( len( lines)):
line = lines[iline].strip()
if len(line) > 0 and not line.startswith('#'):
mat = re.match('^([a-zA-Z0-9]+) *= *(.*)$', line)
if mat == None:
self.throwerr('invalid line. iline: %d line: %s' % (iline, line,),
None)
key = mat.group(1).lower()
val = mat.group(2).lower()
# Convert key to our nomenclature
keyMap = {
'ispin': 'numSpin',
'encut': 'encut_ev',
}
if keyMap.has_key( key): key = keyMap[key]
ix = val.find('#') # get rid of trailing comment
if ix >= 0: val = val[:ix].strip()
val = val.strip('"\'') # strip possibly unbalanced quotes
# Try to convert to int
isConv = False
try:
val = int( val)
isConv = True
except ValueError, exc:
pass
# Try to convert to float
if not isConv:
try:
val = float( val)
isConv = True
except ValueError, exc:
pass
# Try to convert to bool
if not isConv:
tstg = val.lower()
if tstg in ['f', 'false', '.false.']:
val = False
isConv = True
elif tstg in ['t', 'true', '.true.']:
val = True
isConv = True
resObj.__dict__[key] = val
if self.bugLev >= 5:
print 'parseIncar: %s: %s' % (key, val,)
#====================================================================
# We no longer use the POSCAR file since everything
# in it can be derived from the OUTCAR file.
# Fills only a few values:
# resObj.posScaleFactor
# resObj.typeNames
# resObj.typeNums
#
# We get the initial positions from the OUTCAR
# in getInitialPositions.
[docs] def parsePoscar(
self, # self
resObj): # resulting object to be filled
'''
Extracts info from a VASP POSCAR file, setting resObj attributes.
No longer used.
'''
fname = os.path.join( self.inDir, 'POSCAR')
with open( fname) as fin:
lines = fin.readlines()
for ii in range( len( lines)):
lines[ii] = lines[ii].strip()
iline = 0
sysName = lines[iline].strip('"\'') # strip possibly unbalanced quotes
iline += 1
posScaleFactor = float( lines[iline])
iline += 1
basisList = []
for ii in range(3):
line = lines[iline]
iline += 1
toks = line.split()
if len(toks) != 3:
self.throwerr('invalid basis. iline: %d line: %s' % (iline, line,),
None)
basisList.append( map( float, toks))
basisMat = posScaleFactor * np.array( basisList, dtype=float)
# typeNames
# This is not documented at
# http://cms.mpi.univie.ac.at/vasp/vasp/POSCAR_file.html
# Some POSCAR files have an inserted line with the typeNames.
line = lines[iline]
if re.match('^\d+ .*$', line): # if no typeNames, we hit typeNums
typeNames = None
else: # else get typeNames
iline += 1
typeNames = line.split()
# typeNums
line = lines[iline]
iline += 1
toks = line.split()
typeNums = map( int, toks)
# posType: 'direct' or 'cartesian'
posType = lines[iline].lower()
iline += 1
# fracPosMat, cartPosMat
posList = []
for ii in range( sum( typeNums)):
line = lines[iline]
iline += 1
toks = line.split()
if len(toks) != 3:
self.throwerr('invalid pos. iline: %d line: %s' % (iline, line,),
None)
posList.append( map( float, toks))
posMat = np.array( posList, dtype=float)
if posType in ['d', 'direct']:
fracPosMat = posMat
cartPosMat = np.dot( fracPosMat, basisMat)
elif posType in ['c', 'cartesian']:
cartPosMat = posMat * posScaleFactor
fracPosMat = np.dot( cartPosMat, np.linalg.inv( basisMat))
else: self.throwerr('unknown posType: %s' % (posType,), None)
resObj.posScaleFactor = posScaleFactor
resObj.typeNames = typeNames
resObj.typeNums = typeNums
# Ignore the matrix items in POSCAR, as they sometimes
# differ from OUTCAR and the word is that OUTCAR is correct.
# resObj.initialBasisMat = basisMat
# resObj.initialFracPosMat = directMat
# resObj.initialCartPosMat = cartMat
if self.bugLev >= 5:
print 'parsePoscar: posScaleFactor: %s' % (resObj.posScaleFactor,)
print 'parsePoscar: typeNames: %s' % (typeNames,)
print 'parsePoscar: typeNums: %s' % (typeNums,)
#====================================================================
[docs] def getScalars(
self, # self
resObj): # resulting object to be filled
'''
Extracts scalars from OUTCAR, setting resObj attributes.
'''
# k-Points NKPTS = 51 number of bands NBANDS= 26
# or
# k-points NKPTS = 260 k-points in BZ NKDIM = 260 number of bands NBANDS= 11
nkpointPat = r'^k-.oints +NKPTS *= *(\d+) .* number of bands +NBANDS *= *\d+$'
nbandPat = r'^k-.oints +NKPTS *= *\d+ .* number of bands +NBANDS *= *(\d+)$'
# Specs for grepping some scalars from OUTCAR.
# Each spec has the fields:
# tag key for resObj.__dict__
# pat regular expression to match
# which which matched line to use: 'first' first line, or 'last'.
# numMin min num lines which must match pat
# numMax 0 or max num lines which must match pat
# tp result type: int, float, or str
specs = [
# EDIFF = 0.6E-04 stopping-criterion for ELM
Sspec( 'ediff', r'^EDIFF *= *([-.E0-9]+) *stopping-criterion.*$',
'first', 1, 1, float),
# ENCUT = 340.0 eV 24.99 Ry 5.00 a.u. 4.51 4.51 4.51*2*pi/ulx,y,z
Sspec( 'encut_ev', r'^ENCUT *= *([-.E0-9]+) eV .*$',
'first', 1, 1, float),
# IALGO = 68 algorithm
Sspec( 'ialgo', r'^IALGO *= *(\d+) +algorithm$',
'first', 1, 1, int),
# ALGO = Fast
# ALGO =GW execute GW part
Sspec( 'algo', r'^ALGO *= *([a-zA-Z0-9]+).*$',
'first', 0, 1, str),
# IBRION = 2 ionic relax: 0-MD 1-quasi-New 2-CG
Sspec( 'ibrion', r'^IBRION *= *([-0-9]+) +ionic relax *: .*$',
'first', 1, 1, int),
# ICHARG = 0 charge: 1-file 2-atom 10-const
Sspec( 'icharg', r'^ICHARG *= *(\d+) +charge *: .*$',
'first', 1, 1, int),
# ISPIN = 2 # 1: non spin polarized, 2: spin polarized
Sspec( 'numSpin', r'^ISPIN *= *(\d+) +spin polarized .*$',
'first', 1, 1, int),
# ISIF = 3 stress and relaxation
Sspec( 'isif', r'^ISIF *= *(\d+) +stress and relaxation$',
'first', 1, 1, int),
# ISMEAR and SIGMA from:
# ISMEAR = 0; SIGMA = 0.01 broadening in eV -4-tet -1-fermi 0-gaus
# ISMEAR = -5; SIGMA = 0.20 broadening in eV -4-tet -1-fermi 0-gaus
Sspec( 'ismear',
r'^ *ISMEAR *= *([-0-9]+); +SIGMA *= *[-.E0-9]+ +broadening .*$',
'first', 1, 1, int),
Sspec( 'sigma',
r'^ *ISMEAR *= *[-0-9]+; +SIGMA *= *([-.E0-9]+) +broadening .*$',
'first', 1, 1, float),
# LDA+U is selected, type is set to LDAUTYPE = 2
Sspec( 'ldautype',
r'^ *LDA\+U is selected, type is set to LDAUTYPE *= *(\d+)$',
'first', 0, 1, int),
# NELECT = 14.0000 total number of electrons
Sspec( 'numElectron',
r'^NELECT *= *([-.E0-9]+) +total number of electrons$',
'first', 1, 1, float),
# k-points NKPTS = 260 k-points in BZ NKDIM = 260 \
# number of bands NBANDS= 13
Sspec( 'numKpoint', nkpointPat, 'first', 1, 1, int),
Sspec( 'numBand', nbandPat, 'first', 1, 1, int),
# Volume
# volume of cell : 20.0121
Sspec( 'finalVolume_ang3', r'^volume of cell *: +([-.E0-9]+)$',
'last', 1, 0, float),
] # specs
for spec in specs:
ixs = self.findLines( [spec.pat], spec.numMin, spec.numMax)
if len(ixs) > 0:
if spec.which == 'first': ix = ixs[0]
elif spec.which == 'last': ix = ixs[-1]
else: self.throwerr('invalid which in spec: %s' % (spec,), None)
mat = spec.regex.match( self.lines[ix])
value = spec.tp( mat.group(1).lower())
# We can get a value from both parseIncar and OUTCAR
if resObj.__dict__.has_key( spec.tag):
cmsg = wrapUpload.deepCompare(
self.bugLev, 1.e-15,
'INCAR', resObj.__dict__[ spec.tag],
'OUTCAR', value)
if cmsg != None:
print 'value conflict (using OUTCAR): spec: %s cmsg: %s' \
% (spec, cmsg,)
resObj.__dict__[spec.tag] = value # override INCAR value
else:
resObj.__dict__[spec.tag] = value
if self.bugLev >= 5:
print 'getScalars: %s: %s' % (spec.tag, value,)
#====================================================================
# runDate
# executed on LinuxIFC date 2013.10.18 08:44:21
[docs] def getDate(
self, # self
resObj): # resulting object to be filled
'''
Extracts runDate from OUTCAR, setting resObj.runDate.
'''
pat = r'^executed on .* date +\d{4}\.\d{2}\.\d{2} +\d{2}:\d{2}:\d{2}$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
toks = self.lines[ixs[0]].split()
stg = ' '.join( toks[-2:])
resObj.runDate = datetime.datetime.strptime( stg, '%Y.%m.%d %H:%M:%S')
if self.bugLev >= 5:
print 'getDate: runDate: %s' % (
resObj.runDate.strftime('%Y.%m.%d %H:%M:%S'),)
#====================================================================
# iterCpuTimes (not in GW calcs)
# iterRealTimes (not in GW calcs)
# Lines are like either of:
# LOOP+: cpu time 22.49: real time 24.43
# LOOP+: VPU time 42.93: CPU time 43.04
# LOOP+: cpu time16935.65: real time16944.21
# LOOP+: cpu time********: real time********
# In the second line, "VPU" means CPU, and "CPU" means wall time.
# See: http://cms.mpi.univie.ac.at/vasp-forum/forum_viewtopic.php?3.336
#
# In all calcs get:
# Total CPU time used (sec): 11.725
# User time (sec): 9.392
# System time (sec): 2.334
# Elapsed time (sec): 14.468
[docs] def getTimes(
self, # self
resObj): # resulting object to be filled
'''
Extracts timing info from OUTCAR, setting attrs in resObj.
'''
pat = r'^LOOP\+ *: +[a-zA-Z]+ +time *([.0-9]+): +[a-zA-Z]+' \
+ r' +time *([.0-9]+)$'
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
if len(ixs) == 0:
resObj.iterCpuTimes = None
resObj.iterRealTimes = None
else:
cpuTimes = []
realTimes = []
for ix in ixs:
mat = re.match( pat, self.lines[ix])
cpuTimes.append( float( mat.group(1)))
realTimes.append( float( mat.group(2)))
resObj.iterCpuTimes = cpuTimes
resObj.iterRealTimes = realTimes
pat = r'^Total CPU time used \(sec\): +([.0-9]+)$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
mat = re.match( pat, self.lines[ixs[0]])
resObj.totalCpuTimeSec = mat.group(0)
pat = r'^User time \(sec\): +([.0-9]+)$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
mat = re.match( pat, self.lines[ixs[0]])
resObj.userTimeSec = mat.group(0)
pat = r'^System time \(sec\): +([.0-9]+)$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
mat = re.match( pat, self.lines[ixs[0]])
resObj.systemTimeSec = mat.group(0)
pat = r'^Elapsed time \(sec\): +([.0-9]+)$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
mat = re.match( pat, self.lines[ixs[0]])
resObj.elapsedTimeSec = mat.group(0)
if self.bugLev >= 5:
print 'getTimes: iterCpuTimes: %s' % (resObj.iterCpuTimes,)
print 'getTimes: iterRealTimes: %s' % (resObj.iterRealTimes,)
print 'getTimes: totalCpuTimeSec: %s' % (resObj.totalCpuTimeSec,)
print 'getTimes: userTimeSec: %s' % (resObj.userTimeSec,)
print 'getTimes: systemTimeSec: %s' % (resObj.systemTimeSec,)
print 'getTimes: elapsedTimeSec: %s' % (resObj.elapsedTimeSec,)
#====================================================================
# angular momentum for each species, LDAU# 1 : L = 2 -1
# U (eV) for each species, LDAU# 1 : U = -2.40 0.00
# J (eV) for each species, LDAU# 1 : J = 0.00 0.00
# nlep for each species, LDAU# 1 : O = 2 1
# angular momentum for each species, LDAU# 2 : L = -1 -1
# U (eV) for each species, LDAU# 2 : U = 0.00 0.00
# J (eV) for each species, LDAU# 2 : J = 0.00 0.00
# nlep for each species, LDAU# 2 : O = 1 1
# angular momentum for each species, LDAU# 3 : L = -1 -1
# U (eV) for each species, LDAU# 3 : U = 0.00 0.00
# J (eV) for each species, LDAU# 3 : J = 0.00 0.00
# nlep for each species, LDAU# 3 : O = 1 1
# angular momentum for each species, LDAU# 4 : L = -1 -1
# U (eV) for each species, LDAU# 4 : U = 0.00 0.00
# J (eV) for each species, LDAU# 4 : J = 0.00 0.00
# nlep for each species, LDAU# 4 : O = 1 1
def getLdau(
self, # self
numType, # num atomic types
resObj): # resulting object to be filled
'''
Extracts LDAU info from OUTCAR, setting attrs in resObj.
'''
pat = r' *angular momentum +for each species, LDAU# *(\d+) *: *L *=.*$'
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
numLdau = len( ixs)
if numLdau == 0:
ldaul_mat = None # numLdau * numType
ldauu_mat = None # numLdau * numType
ldauj_mat = None # numLdau * numType
ldauo_mat = None # numLdau * numType
else:
ldaul_mat = numLdau * [None] # numLdau * numType
ldauu_mat = numLdau * [None] # numLdau * numType
ldauj_mat = numLdau * [None] # numLdau * numType
ldauo_mat = numLdau * [None] # numLdau * numType
for ii in range(len(ixs)):
# Get LDAUL line
ix = ixs[ii]
line = self.lines[ix]
if not re.match(r'^ *angular momentum +for each species, LDAU# +\d+ *: *L *= .*$',
line):
self.throwerr('ldaul pat mismatch', ix)
toks = line.strip().split()
if len(toks) != 10 + numType:
self.throwerr('ldaul len mismatch', ix)
if int(toks[6]) != ii+1:
self.throwerr('ldaul grp mismatch', ix)
vals = map( int, toks[10:])
ldaul_mat[ii] = vals
# Get LDAUU line
ix = ixs[ii] + 1
line = self.lines[ix]
if not re.match(r'^ *U \(eV\) +for each species, LDAU# +\d+ *: *U *= .*$',
line):
self.throwerr('ldauu pat mismatch', ix)
toks = line.strip().split()
if len(toks) != 10 + numType:
self.throwerr('ldauu len mismatch', ix)
if int(toks[6]) != ii+1:
self.throwerr('ldauu grp mismatch', ix)
vals = map( float, toks[10:])
ldauu_mat[ii] = vals
# Get LDAUJ line
ix = ixs[ii] + 2
line = self.lines[ix]
if not re.match(r'^ *J \(eV\) +for each species, LDAU# +\d+ *: *J *= .*$',
line):
self.throwerr('ldauj pat mismatch', ix)
toks = line.strip().split()
if len(toks) != 10 + numType:
self.throwerr('ldauj len mismatch', ix)
if int(toks[6]) != ii+1:
self.throwerr('ldauj grp mismatch', ix)
vals = map( float, toks[10:])
ldauj_mat[ii] = vals
# Get LDAUO line, if any
ix = ixs[ii] + 3
line = self.lines[ix]
if not re.match(r'^ *nlep +for each species, LDAU# +\d+ *: *O *= .*$',
line):
ldauo_mat[ii] = None
else:
toks = line.strip().split()
if len(toks) != 9 + numType:
self.throwerr('ldauo len mismatch', ix)
if int(toks[5]) != ii+1:
self.throwerr('ldauo grp mismatch', ix)
vals = map( int, toks[9:])
ldauo_mat[ii] = vals
if self.bugLev >= 5:
print 'getLdau: raw ldaul_mat:\n%s' % (ldaul_mat,)
print 'getLdau: raw ldauu_mat:\n%s' % (ldauu_mat,)
print 'getLdau: raw ldauj_mat:\n%s' % (ldauj_mat,)
print 'getLdau: raw ldauo_mat:\n%s' % (ldauo_mat,)
ldaul_mat = np.array( ldaul_mat)
ldauu_mat = np.array( ldauu_mat)
ldauj_mat = np.array( ldauj_mat)
ldauo_mat = np.array( ldauo_mat)
resObj.ldaul_mat = ldaul_mat
resObj.ldauu_mat = ldauu_mat
resObj.ldauj_mat = ldauj_mat
resObj.ldauo_mat = ldauo_mat
if self.bugLev >= 5:
print 'getLdau: ldaul_mat:\n%s' % (resObj.ldaul_mat,)
print 'getLdau: ldauu_mat:\n%s' % (resObj.ldauu_mat,)
print 'getLdau: ldauj_mat:\n%s' % (resObj.ldauj_mat,)
print 'getLdau: ldauo_mat:\n%s' % (resObj.ldauo_mat,)
#====================================================================
# system
# SYSTEM = icsd_633029 in icsd_633029.cif, hs-ferr
[docs] def getSystem(
self, # self
resObj): # resulting object to be filled
'''
Extracts systemName from OUTCAR, setting resObj.systemName
'''
pat = r'^SYSTEM *= '
ixs = self.findLines( [pat], 1, 0) # pats, numMin, numMax
toks = self.lines[ixs[0]].split()
nm = ' '.join(toks[2:])
nm = nm.strip('\"\'') # strip unpaired surrounding quotes
resObj.systemName = nm
if self.bugLev >= 5:
print 'getSystem: systemName: %s' % (resObj.systemName,)
#====================================================================
# For doc on ALGO see:
# http://cms.mpi.univie.ac.at/vasp/vasp/Recipe_selfconsistent_GW_calculations.html
# http://cms.mpi.univie.ac.at/vasp/vasp/Recipe_GW_calculations.html
# http://cms.mpi.univie.ac.at/wiki/index.php/ALGO
[docs] def getCategEnergy(
self, # self
resObj): # resulting object to be filled
'''
Extracts energyNoEntrps, efermi0 from OUTCAR, setting resObj attributes.
'''
# Spacing varies:
# energy without entropy = -60.501550 energy(sigma->0) = -60.501640
# energy without entropy= -15.051005 energy(sigma->0) = -15.051046
pat = r'^energy +without +entropy *= *([-.E0-9]+)' \
+ r' +energy\(sigma->0\) *= *[-.E0-9]+$'
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
if len(ixs) == 0:
resObj.energyNoEntrps = None
resObj.efermi0 = None
else:
resObj.energyNoEntrps = []
for ix in ixs:
mat = re.match( pat, self.lines[ix])
resObj.energyNoEntrps.append( float( mat.group(1)))
# E-fermi : 6.5043 XC(G=0): -12.1737 alpha+bet :-11.7851
pat = r'^E-fermi *: *([-.E0-9]+) +XC.*alpha\+bet.*$'
ixs = self.findLines( [pat], 1, 0) # pats, numMin, numMax
ix = ixs[-1]
mat = re.match( pat, self.lines[ix])
resObj.efermi0 = float( mat.group(1))
#====================================================================
# initialBasisMat, finalBasisMat == cell
# direct lattice vectors reciprocal lattice vectors
# 0.0004124 2.1551750 2.1553184 -0.2320833 0.2320088 0.2320195
# 2.1551376 0.0004322 2.1552987 0.2320129 -0.2320836 0.2320237
# 2.1548955 2.1549131 0.0006743 0.2320650 0.2320652 -0.2320942
[docs] def getBasisMats(
self, # self
resObj): # resulting object to be filled
'''
Extracts basis matrices from OUTCAR, setting resObj attributes.
'''
pat = r'^direct lattice vectors +reciprocal lattice vectors$'
ixs = self.findLines( [pat], 2, 0) # pats, numMin, numMax
ix = ixs[0] # initial
initialBasisMat = self.parseMatrix(
6, ix+1, ix+4, 0, 3) # ntok, rowBeg, rowEnd, colBeg, colEnd
initialRecipBasisMat = self.parseMatrix(
6, ix+1, ix+4, 3, 6) # ntok, rowBeg, rowEnd, colBeg, colEnd
# The basisMats can differ from POSCAR by 1.e-1 or so. Why?
# Use the OUTCAR version.
resObj.initialBasisMat = initialBasisMat
resObj.initialRecipBasisMat = initialRecipBasisMat
ix = ixs[-1] # final
# In vasprun.xml and OUTCAR, the basis vectors are rows.
# In PyLada, Mayeul uses the transpose.
resObj.finalBasisMat = self.parseMatrix(
6, ix+1, ix+4, 0, 3) # ntok, rowBeg, rowEnd, colBeg, colEnd
resObj.finalRecipBasisMat = self.parseMatrix(
6, ix+1, ix+4, 3, 6) # ntok, rowBeg, rowEnd, colBeg, colEnd
# Test the linear algebra
test_initialRecip = np.linalg.inv( resObj.initialBasisMat).T
test_finalRecip = np.linalg.inv( resObj.finalBasisMat).T
if self.bugLev >= 5:
print 'getBasisMats: initialBasisMat:\n%s' \
% (repr( resObj.initialBasisMat),)
print 'getBasisMats: initialRecipBasisMat:\n%s' \
% (repr( resObj.initialRecipBasisMat),)
print 'getBasisMats: test_initialRecip:\n%s' \
% (repr( test_initialRecip),)
print 'getBasisMats: finalBasisMat:\n%s' \
% (repr( resObj.finalBasisMat),)
print 'getBasisMats: finalRecipBasisMat:\n%s' \
% (repr( resObj.finalRecipBasisMat),)
print 'getBasisMats: test_finalRecip:\n%s' \
% (repr( test_finalRecip),)
# Test the linear algebra
self.checkClose( resObj.initialRecipBasisMat, test_initialRecip, 1.e-3)
self.checkClose( resObj.finalRecipBasisMat, test_finalRecip, 1.e-3)
#====================================================================
# typeNames == species
# All lines like:
# VRHFIN =Fe: d7 s1
# VRHFIN =O: s2p4
# VRHFIN =O : s2p4
# VRHFIN =Ni:
#
# Other possible patterns we could use, but currently do not.
# POTCAR: PAW_PBE Zn 06Sep2000 # no, repeated
# TITEL = PAW_PBE Zn 06Sep2000
[docs] def getTypeNames(
self, # self
resObj): # resulting object to be filled
'''
Extracts typeNames from OUTCAR, setting resObj attribute.
'''
pat = r'^VRHFIN *= *([a-zA-Z]+) *:'
regexa = re.compile( pat)
typeNames = []
for iline in range( self.numLine):
line = self.lines[iline]
mat = regexa.match( line)
if mat != None:
typeNames.append( mat.group( 1))
# If POSCAR didn't have typeNames
if getattr( resObj, 'typeNames', None) == None:
resObj.typeNames = typeNames
else:
if typeNames != resObj.typeNames:
self.throwerr('typeNames mismatch. POSCAR: %s OUTCAR: %s' \
% (resObj.typeNames, typeNames,), None)
if self.bugLev >= 5:
print 'getTypeNames: typeNames: %s' % (resObj.typeNames,)
#====================================================================
# typeNums == stoichiometry
# ions per type = 1 1
[docs] def getTypeNums(
self, # self
resObj): # resulting object to be filled
'''
Extracts typeNums (ions per type) from OUTCAR, setting resObj attributes.
'''
pat = r'^ions per type *= '
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
toks = self.lines[ixs[0]].split()
typeNums = map( int, toks[4:])
# If POSCAR didn't have typeNums
if getattr( resObj, 'typeNums', None) == None:
resObj.typeNums = typeNums
else:
if typeNums != resObj.typeNums:
self.throwerr('typeNums mismatch. POSCAR: %s OUTCAR: %s' \
% (resObj.typeNums, typeNums,), None)
resObj.numAtom = sum( resObj.typeNums)
if self.bugLev >= 5:
print 'getTypeNums: typeNums: %s' % (resObj.typeNums,)
print 'getTypeNums: numAtom: %d' % (resObj.numAtom,)
#====================================================================
# initialFracPosMat, initialCartPosMat
# position of ions in fractional coordinates (direct lattice)
# 0.00000005 0.99999999 0.99999997
# 0.49999995 0.50000001 0.50000003
#
# position of ions in cartesian coordinates (Angst):
# 4.24342506 2.12183357 2.12263079
# 2.12184258 2.12186596 2.12233823
[docs] def getInitialPositions(
self, # self
resObj): # resulting object to be filled
'''
Extracts initial position matrices from OUTCAR, setting resObj attributes.
'''
patf = r'^position of ions in fractional coordinates \(direct lattice\)$'
patc = r'^position of ions in cartesian coordinates +\(Angst\):$'
# Get initialFracPosMat
ixs = self.findLines( [patf], 1, 1) # pats, numMin, numMax
ix = ixs[0]
resObj.initialFracPosMat = self.parseMatrix(
3, ix+1, ix+1+resObj.numAtom, 0, 3)
# ntok, rowBeg, rowEnd, colBeg, colEnd
# The posMats in OUTCAR and POSCAR can differ by 1.e-1 or so. Why?
# Use the OUTCAR version.
# Get initialCartPosMat
ixs = self.findLines( [patc], 1, 1) # pats, numMin, numMax
ix = ixs[0]
resObj.initialCartPosMat = self.parseMatrix(
3, ix+1, ix+1+resObj.numAtom, 0, 3)
# ntok, rowBeg, rowEnd, colBeg, colEnd
# The posMats in OUTCAR and POSCAR can differ by 1.e-1 or so. Why?
# Use the OUTCAR version.
# Test the linear algebra
test_initialFracPosMat = np.dot(
resObj.initialCartPosMat, np.linalg.inv( resObj.initialBasisMat))
test_initialCartPosMat = np.dot(
resObj.initialFracPosMat, resObj.initialBasisMat)
self.checkClose( resObj.initialFracPosMat, test_initialFracPosMat, 1.e-3)
self.checkClose( resObj.initialCartPosMat, test_initialCartPosMat, 1.e-3)
if self.bugLev >= 5:
print 'getInitialPositions: initialFracPosMat:\n%s' \
% (self.formatMatrix( resObj.initialFracPosMat),)
print 'getInitialPositions: initialCartPosMat:\n%s' \
% (self.formatMatrix( resObj.initialCartPosMat),)
#====================================================================
# finalCartPosMat, forceMat
# POSITION TOTAL-FORCE (eV/Angst)
# ------------------------------------------------------------------
# 0.59671 5.05062 2.21963 0.003680 -0.001241 -0.004039
# 4.98024 2.21990 0.58287 -0.000814 -0.015091 -0.012495
# 2.23430 0.66221 5.03287 0.017324 0.003318 -0.014394
# ...
# ------------------------------------------------------------------
#
# Not found in gw calcs.
[docs] def getFinalPositionsForce(
self, # self
resObj): # resulting object to be filled
'''
Extracts final positions and forces from OUTCAR, setting resObj attributes.
'''
pat = r'^POSITION +TOTAL-FORCE \(eV/Angst\)$'
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
if len(ixs) == 0:
resObj.finalFracPosMat = resObj.initialFracPosMat
resObj.finalCartPosMat = resObj.initialCartPosMat
resObj.finalForceMat_ev_ang = None
else:
ix = ixs[-1] # last one
resObj.finalCartPosMat = self.parseMatrix(
6, ix+2, ix+2+resObj.numAtom, 0, 3)
# ntok, rowBeg, rowEnd, colBeg, colEnd
resObj.finalForceMat_ev_ang = self.parseMatrix(
6, ix+2, ix+2+resObj.numAtom, 3, 6)
# ntok, rowBeg, rowEnd, colBeg, colEnd
resObj.finalFracPosMat = np.dot(
resObj.finalCartPosMat, np.linalg.inv( resObj.finalBasisMat))
if self.bugLev >= 5:
print 'getFinalPositionsForce: finalCartPosMat:\n%s' \
% (self.formatMatrix( resObj.finalCartPosMat),)
print 'getFinalPositionsForce: finalFracPosMat:\n%s' \
% (self.formatMatrix( resObj.finalFracPosMat),)
print 'getFinalPositionsForce: finalForceMat_ev_ang:\n%s' \
% (self.formatMatrix( resObj.finalForceMat_ev_ang),)
#====================================================================
# stressMat
# FORCE on cell =-STRESS in cart. coord. units (eV):
# Direction XX YY ZZ XY YZ ZX
# -----------------------------------------------------------------------
# Alpha Z 76.60341 76.60341 76.60341
# Ewald -318.46783 -318.46428 -318.52959 -0.00310 0.03083 0.02783
# Hartree 63.38495 63.38593 63.35616 0.00114 0.01511 0.01368
# E(xc) -70.45894 -70.45888 -70.45878 -0.00052 -0.00008 -0.00008
# Local 3.47019 3.46661 3.56011 -0.00319 -0.04620 -0.04172
# n-local -6.52091 -6.52080 -6.52183 -0.00373 -0.00131 -0.00134
# augment 82.46862 82.46871 82.46910 0.00638 0.00129 0.00151
# Kinetic 170.67238 170.67209 170.67398 0.00710 0.00098 0.00116
# Fock 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
# -----------------------------------------------------------------------
# Total 1.15189 1.15278 1.15256 0.00410 0.00063 0.00105
# in kB 96.62361 96.69864 96.68018 0.34362 0.05255 0.08820
# external pressure = 96.67 kB Pullay stress = 0.00 kB
#
# Not in gw calcs.
[docs] def getStressMat(
self, # self
resObj): # resulting object to be filled
'''
Extracts final stress matrix from OUTCAR, setting resObj attributes.
'''
pat = r'^FORCE on cell *= *-STRESS in cart. coord. +units \(eV\) *:$'
# In some non-GW calcs with ALGO=DIAG there is no stress matrix.
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
if len(ixs) == 0:
resObj.finalStressMat_ev = None
resObj.finalStressMat_kbar = None
else:
ix = ixs[-1] # last one
# Find the "Total" and "in kB" lines
evline = None
kbline = None
for ii in range(ix+1, ix + 16):
line = self.lines[ii]
if re.match(r'^Total( +[-.E0-9]+){6}$', line): evline = ii
if re.match(r'^in kB( +[-.E0-9]+){6}$', line): kbline = ii
if evline == None: self.throwerr('unknown stress ev', ix)
if kbline == None: self.throwerr('unknown stress kb', ix)
# Extract tokens from the "Total" and "in kB" lines
evtoks = self.lines[ evline].split()[1:] # skip 'Total'
kbtoks = self.lines[ kbline].split()[2:] # skip 'in kB'
if len( evtoks) != 6: self.throwerr('unknown stress evtoks', evline)
if len( kbtoks) != 6: self.throwerr('unknown stress kbtoks', kbline)
evs = map( float, evtoks)
kbs = map( float, kbtoks)
# Create 3x3 matrices for ev and kb stress
evstress = np.zeros( [3, 3], dtype=float)
kbstress = np.zeros( [3, 3], dtype=float)
# Add in the terms XX, YY, ZZ
for ii in range(3):
evstress[ii,ii] += evs[ii]
kbstress[ii,ii] += kbs[ii]
evstress[1,0] = evstress[0,1] = evs[3] # XY
kbstress[1,0] = kbstress[0,1] = kbs[3] # XY
evstress[2,0] = evstress[0,2] = evs[5] # XY
kbstress[2,0] = kbstress[0,2] = kbs[5] # XY
evstress[2,1] = evstress[1,2] = evs[4] # YZ
kbstress[2,1] = kbstress[1,2] = kbs[4] # YZ
resObj.finalStressMat_ev = evstress
resObj.finalStressMat_kbar = kbstress
if self.bugLev >= 5:
print 'getStressMat: finalStressMat_ev:\n%s' \
% (self.formatMatrix( resObj.finalStressMat_ev),)
print 'getStressMat: finalStressMat_kbar:\n%s' \
% (self.formatMatrix( resObj.finalStressMat_kbar),)
#====================================================================
## Old typeValences == ionic_charges
## ZVAL = 8.00 6.00
#def getTypeValencesUnused( self, resObj):
# pat = r'^ZVAL *= '
# ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
# toks = self.lines[ixs[0]].split()
# resObj.typeValences = map( float, toks[2:])
#====================================================================
# typeMasses, typeValences
# POMASS = 55.847; ZVAL = 8.000 mass and valenz
# POMASS = 16.000; ZVAL = 6.000 mass and valenz
# POMASS = 55.85 16.00 # ignore this
#
[docs] def getTypeMassValence(
self, # self
resObj): # resulting object to be filled
'''
Extracts masses and valences from OUTCAR, setting resObj attributes.
'''
pat = r'^POMASS *= *([.0-9]+); +ZVAL += +([.0-9]+) +mass and valenz'
ixs = self.findLines( [pat], 1, 0) # pats, numMin, numMax
typeMasses = []
typeValences = []
for ix in ixs:
mat = re.match( pat, self.lines[ix])
typeMasses.append( float( mat.group(1)))
typeValences.append( float( mat.group(2)))
resObj.typeMasses_amu = typeMasses
resObj.typeValences = typeValences
if self.bugLev >= 5:
print 'getTypeMassValence: typeMasses_amu: %s' % (resObj.typeMasses_amu,)
print 'getTypeMassValence: typeValences: %s' % (resObj.typeValences,)
#====================================================================
# typePseudos
# TITEL = PAW_PBE Fe 06Sep2000
# TITEL = PAW_PBE O_s 07Sep2000
[docs] def getTypePseudos(
self, # self
resObj): # resulting object to be filled
'''
Extracts pseudopotential names from OUTCAR, setting resObj attributes.
'''
pat = r'^TITEL *= *(PAW.*)$'
ixs = self.findLines( [pat], 1, 0) # pats, numMin, numMax
typePseudos = []
for ix in ixs:
mat = re.match( pat, self.lines[ix])
typePseudos.append( mat.group(1).strip())
resObj.typePseudos = typePseudos
if self.bugLev >= 5:
print 'getTypePseudos: typePseudos: %s' % (resObj.typePseudos,)
#====================================================================
# Magmom, magnetic moment, magnetization:
# Returns the 'tot' column.
# Not in vasp 4.6.
#
# magnetization (x)
#
# # of ion s p d tot
# ----------------------------------------
# 1 -0.017 -0.008 -1.935 -1.960
# 2 -0.005 -0.006 0.000 -0.011
# ------------------------------------------------
# tot -0.022 -0.014 -1.935 -1.970
#
# But if there is only one atom there is no "tot" line:
# magnetization (x)
#
# # of ion s p d tot
# ----------------------------------------
# 1 -0.008 -0.014 0.549 0.526
# (blank line)
[docs] def getMagmoms(
self, # self
resObj): # resulting object to be filled
'''
Extracts magnetic moments from OUTCAR, setting resObj attributes.
'''
pat = r'^ *magnetization +\(x\)$'
ixs = self.findLines( [pat], 0, 0) # pats, numMin, numMax
if len(ixs) == 0:
resObj.magmoms = None
else:
ix = ixs[-1] # use the last one
ix += 1; # skip over header 'magnetization (x)'
# Skip over blank lines
while ix < self.numLine and self.lines[ix].strip() == '':
ix += 1
# Skip over title '# of ion' and '-----'
if not (self.lines[ix].strip().startswith('# of ion ')
and self.lines[ix+1].strip().startswith('-----')):
self.throwerr('magmom format unknown', ix)
ix += 2
# Get data lines
magmomTots = []
while ix < self.numLine \
and not len( self.lines[ix].strip()) == 0 \
and not self.lines[ix].strip().startswith('-----'):
line = self.lines[ix]
# Get two groups: ionNum, restOfLine
mat = re.match('^ *(\d+)( +[-.0-9]+)+ *$', line)
if not mat: self.throwerr('magmom format unknown', ix)
ionNum = int( mat.group(1))
magStgs = mat.group(2).strip().split()
magvals = map( float, magStgs)
magmomTots.append( magvals[-1]) # just keep the tot column
ix += 1
if len( magmomTots) != resObj.numAtom:
self.throwerr('magmom num mismatch', ix)
resObj.magmoms = magmomTots
if self.bugLev >= 5:
print 'getMagmoms: magmoms: %s' % (resObj.magmoms,)
#====================================================================
# kpointFracMat, kpointCartMat,
# kpointMults, kpointWeights == normalized kpointMults
# Found 260 irreducible k-points:
#
# Following reciprocal coordinates:
# Coordinates Weight
# 0.000000 0.000000 0.000000 1.000000
# 0.125000 0.000000 0.000000 2.000000
# 0.250000 0.000000 0.000000 2.000000
# ...
# 0.500000 0.500000 0.500000 1.000000
#
# Following cartesian coordinates:
# Coordinates Weight
# 0.000000 0.000000 0.000000 1.000000
# -0.029463 0.029454 0.029457 2.000000
# -0.058925 0.058909 0.058914 2.000000
# ...
# 0.117822 0.117821 0.117795 1.000000
[docs] def getKpointMat(
self, # self
resObj): # resulting object to be filled
'''
Extracts kpoints, weights, mults from OUTCAR, setting resObj attributes.
'''
headPat = r'^Found +\d+ +irreducible k-points *:$'
blankPat = r'^$'
recipPat = r'^Following reciprocal coordinates *:$'
cartPat = r'^Following cartesian coordinates *:$'
wtPat = r'^Coordinates +Weight$'
numKpoint = resObj.numKpoint
ixs = self.findLines(
[headPat, blankPat, recipPat, wtPat], 1, 0) # pats, numMin, numMax
# Get kpointFracMat == generated reciprocol fractional coords.
# This is what vasprun.xml calls "kpoints".
# If there are multiple occurances use the first.
ix = ixs[0] + 4 # start of matrix
kpMat = self.parseMatrix( # kpoints and weights (multiplicities)
4, ix, ix+numKpoint, 0, 4) # ntok, rowBeg, rowEnd, colBeg, colEnd
resObj.kpointFracMat = kpMat[ :, 0:3]
kpMults = kpMat[ :, 3]
resObj.kpointMults = kpMults
resObj.kpointWeights = kpMults / float( kpMults.sum())
if self.bugLev >= 5:
print 'getKpointMat: gen frac: kpointMults: %s' \
% (resObj.kpointMults,)
print 'getKpointMat: gen frac: kpointWeights: %s' \
% (resObj.kpointWeights,)
print 'getKpointMat: gen frac: kpointFracMat:\n%s' \
% (resObj.kpointFracMat,)
# Get generated cartesian coords
ixCart = ix + numKpoint + 1
ixs = self.findLines( [cartPat, wtPat], 1, 0) # pats, numMin, numMax
if ixs[0] != ixCart: self.throwerr('kpoint format', ixCart)
ix = ixCart + 2 # start of matrix
kpMat = self.parseMatrix( # kpoints and weights (multiplicities)
4, ix, ix+numKpoint, 0, 4) # ntok, rowBeg, rowEnd, colBeg, colEnd
resObj.kpointCartMat = kpMat[ :, 0:3]
test_kpointMults = kpMat[ :, 3]
self.checkClose( resObj.kpointMults, test_kpointMults, 1.e-3)
if self.bugLev >= 5:
print 'getKpointMat: gen frac: test_kpointMults: %s' \
% (test_kpointMults,)
print 'getKpointMat: gen cart: kpointCartMat:\n%s' \
% (resObj.kpointCartMat,)
# Back-calculate the posScaleFactor, method A,
# using kpointFracMat.
# We use method B, further below, to set resObj.posScaleFactor.
#
# We should have:
# kpointCartMat == posScaleFactor * np.dot(
# kpointFracMat, initialRecipBasisMat)
#
# So posScaleFactor = median( kpointCartMat / dotProd)
dotProd = np.dot( resObj.kpointFracMat, resObj.initialRecipBasisMat)
maskDotProd = np.ma.masked_values( dotProd, 0)
maskKpoint = np.ma.array( resObj.kpointCartMat, mask=maskDotProd.mask)
if self.bugLev >= 5:
print 'getKpointMat: kpointFracMat shape: %s' \
% (resObj.kpointFracMat.shape,)
print 'getKpointMat: kpointCartMat shape: %s' \
% (resObj.kpointCartMat.shape,)
print 'getKpointMat: initialRecipBasisMat shape: %s' \
% (resObj.initialRecipBasisMat.shape,)
print 'getKpointMat: dotProd shape: %s' % (dotProd.shape,)
ratioMat = maskKpoint / maskDotProd
ratioMean = np.ma.mean( ratioMat)
ratioMedian = np.ma.median( ratioMat)
ratioMin = np.ma.min( ratioMat)
ratioMax = np.ma.max( ratioMat)
if self.bugLev >= 5:
print 'getKpointMat: posScaleFrac ratio method A:\n%s' % (ratioMat,)
print ' A: ratioMean: %g' % (ratioMean,)
print ' A: ratioMedian: %g' % (ratioMedian,)
print ' A: ratioMin: %g' % (ratioMin,)
print ' A: ratioMax: %g' % (ratioMax,)
print ' A: ratioDelta: %g' % (ratioMax - ratioMin,)
# Get the actual kpoints
#
# k-points in units of 2pi/SCALE and weight:
# or k-points in units of 2pi/SCALE and weight: Automatic generation
# 0.00000000 0.00000000 0.00000000 0.016
# 0.08286762 -0.04784364 0.00000000 0.094
# ...
# 0.16573524 0.00000000 0.15507665 0.094
#
# k-points in reciprocal lattice and weights:
# or k-points in reciprocal lattice and weights: Automatic generation
# 0.00000000 0.00000000 0.00000000 0.016
# 0.25000000 0.00000000 0.00000000 0.094
# ...
# Read actual_kpointCartMat_a
pat = r'^k-points in units of 2pi/SCALE and weight:' \
+ '( Automatic generation)?$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
ix = ixs[0] + 1
kpMat = self.parseMatrix( # kpoints and weights (multiplicities)
4, ix, ix+numKpoint, 0, 4) # ntok, rowBeg, rowEnd, colBeg, colEnd
actual_kpointCartMat_a = kpMat[ :, 0:3]
wts = kpMat[ :, 3]
# Sometimes the wts are not normalized to 1, for example in
# icsd_633029.cif/non-magnetic/relax_cellshape/1
# the sum is 1.024
actual_kpointWeights_a = wts / sum( wts)
# Back-calculate the posScaleFactor, method B,
# using actual_kpointCartMat_a.
# Method B has narrower ratioDelta values than either
# method A (above) or method C (below), so we use
# method B to set resObj.posScaleFactor.
#
# We should have:
# kpointFracMat = (1. / posScaleFactor) * np.dot(
# kpointCartMat, inv( initialRecipBasisMat))
#
# So posScaleFactor = median( dotProd / kpointFracMat)
# cartMat is numAtom x 3
# recipBasisMat is 3 x 3
# dotProd is numAtom x 3
dotProd = np.dot(
actual_kpointCartMat_a, np.linalg.inv( resObj.initialRecipBasisMat))
# maskKp is numAtom x 3
# maskDot is numAtom x 3
# ratioMat is numAtom x 3
# But for one atom, with position [0,0,0],
# dotProd == maskDot == maskKp == [--,--,--]
maskKp = np.ma.masked_values( resObj.kpointFracMat, 0)
maskDot = np.ma.array( dotProd, mask=maskKp.mask)
ratioMat = maskDot / maskKp
ratioMean = np.ma.mean( ratioMat)
ratioMedian = np.ma.median( ratioMat)
ratioMin = np.ma.min( ratioMat)
ratioMax = np.ma.max( ratioMat)
if self.bugLev >= 5:
print 'getKpointMat: posScaleFrac ratio method B:'
print ' B: actual_kpointCartMat_a:\n%s' % (actual_kpointCartMat_a,)
print ' B: kpointFracMat:\n%s' % (resObj.kpointFracMat,)
print ' B: resObj.initialRecipBasisMat:\n%s' \
% (resObj.initialRecipBasisMat,)
print ' B: inv( resObj.initialRecipBasisMat):\n%s' \
% (np.linalg.inv( resObj.initialRecipBasisMat),)
print ' B: dotProd:\n%s' % (dotProd,)
print ' B: maskDot:\n%s' % (maskDot,)
print ' B: maskKp:\n%s' % (maskKp,)
print ' B: ratioMat:\n%s' % (ratioMat,)
print ' B: ratioMean: %g' % (ratioMean,)
print ' B: ratioMedian: %g' % (ratioMedian,)
print ' B: ratioMin: %g' % (ratioMin,)
print ' B: ratioMax: %g' % (ratioMax,)
print ' B: ratioDelta: %g' % (ratioMax - ratioMin,)
# If there is only one atom, at [0,0,0],
# the median will be masked. But in this case the
# posScaleFactor is meaningless, so call it 1.0.
if type( ratioMedian).__name__ == 'MaskedConstant': ratioMedian = 1.0
resObj.posScaleFactor = ratioMedian
# Read actual_kpointFracMat_b
pat = r'^k-points in reciprocal lattice and weights:' \
+ '( Automatic generation)?$'
ixs = self.findLines( [pat], 1, 1) # pats, numMin, numMax
ix = ixs[0] + 1
kpMat = self.parseMatrix( # kpoints and weights (multiplicities)
4, ix, ix+numKpoint, 0, 4) # ntok, rowBeg, rowEnd, colBeg, colEnd
actual_kpointFracMat_b = kpMat[ :, 0:3]
wts = kpMat[ :, 3]
# Sometimes the wts are not normalized to 1, for example in
# icsd_633029.cif/non-magnetic/relax_cellshape/1
# the sum is 1.024
actual_kpointWeights_b = wts / sum( wts)
# Back-calculate the posScaleFactor, method C,
# using kpointFracMat_b.
# We use method B, further above, to set resObj.posScaleFactor.
#
# We should have:
# kpointCartMat == posScaleFactor * np.dot(
# kpointFracMat, initialRecipBasisMat)
#
# So posScaleFactor = median( kpointCartMat / dotProd)
dotProd = np.dot( actual_kpointFracMat_b, resObj.initialRecipBasisMat)
maskDotProd = np.ma.masked_values( dotProd, 0)
maskKpoint = np.ma.array( resObj.kpointCartMat, mask=maskDotProd.mask)
ratioMat = maskKpoint / maskDotProd
ratioMean = np.ma.mean( ratioMat)
ratioMedian = np.ma.median( ratioMat)
ratioMin = np.ma.min( ratioMat)
ratioMax = np.ma.max( ratioMat)
if self.bugLev >= 5:
print 'getKpointMat: posScaleFrac ratio method C:\n%s' % (ratioMat,)
print ' C: ratioMean: %g' % (ratioMean,)
print ' C: ratioMedian: %g' % (ratioMedian,)
print ' C: ratioMin: %g' % (ratioMin,)
print ' C: ratioMax: %g' % (ratioMax,)
print ' C: ratioDelta: %g' % (ratioMax - ratioMin,)
# This doesn't work. Some sets of weights don't work this way.
# For example, lany.projects.ppmdd.pub_data.2013_PRB_GGAUvsHSE/
# GGAU_defects/Al1N1/Vc_Al/OUTCAR gives
# uniqWts: [0.015952143569292122, 0.030907278165503486,
# 0.046859421734795605, 0.093718843469591209, 0.18743768693918242]
# fac = 1/uniqWts[0] = 62.6875
# uniqWts * fac are not integers
#
## Find kpoint multiplicities by multiplying kpointWeights
## by the min value that makes them all integers
#wts = actual_kpointWeights_a
#uniqWts = list( set( wts))
#uniqWts.sort()
#minWt = uniqWts[0]
#if self.bugLev >= 5:
# print 'getKpointMat: wts: %s' % (wts,)
# print 'getKpointMat: uniqWts: %s' % (uniqWts,)
# print 'getKpointMat: minWt: %s' % (minWt,)
#factor = None
#for ii in range(1, 10):
# tfact = ii / float( minWt)
# allOk = True
# for wt in uniqWts:
# val = tfact * wt
# if abs( val - round(val)) > 1.e-5:
# allOk = False
# break
# if allOk:
# factor = tfact
# break
#if factor == None: self.throwerr('no kpointMults found', None)
#actual_kpointMults_a = wts * factor
#self.checkClose( resObj.kpointMults, actual_kpointMults_a, 1.e-3)
if self.bugLev >= 5:
print 'getKpointMat: numKpoint: %d' % (numKpoint,)
print 'getKpointMat: kpointMults sum: %g' \
% (resObj.kpointMults.sum(),)
print 'getKpointMat: kpointWeights sum: %g' \
% (resObj.kpointWeights.sum(),)
#====================================================================
# eigenMat, occupMat
#
# spin component 1 # this line is not in non-mag calcs
#
# k-point 1 : 0.0000 0.0000 0.0000
# band No. band energies occupation
# 1 -14.2072 1.00000
# 2 2.5938 1.00000
# ...
# 13 23.5652 0.00000
#
# k-point 2 : 0.1250 0.0000 0.0000
# band No. band energies occupation
# 1 -14.0551 1.00000
# 2 1.3327 1.00000
# ...
# 13 23.5652 0.00000
# ...
# k-point 260 : 0.5000 0.5000 0.5000
# band No. band energies occupation
# 1 -13.0330 1.00000
# 2 -2.2402 1.00000
# ...
# 13 19.0910 0.00000
#
# spin component 2
#
# k-point 1 : 0.0000 0.0000 0.0000
# band No. band energies occupation
# 1 -13.6767 1.00000
# 2 3.2576 1.00000
# ...
[docs] def getEigenMat(
self, # self
resObj): # resulting object to be filled
'''
Extracts eigenvalues and occupancies from OUTCAR, setting resObj attributes.
'''
numSpin = resObj.numSpin
numKpoint = resObj.numKpoint
numBand = resObj.numBand
compPat = '^spin component [12]$'
kpFirstPat = '^k-point +1 *: +[-.E0-9]+ +[-.E0-9]+ +[-.E0-9]+$'
eigenMat = None # Wait until we know numEigen to allocate
occupMat = None # Wait until we know numEigen to allocate
# 2014-06-17:In some GW runs, there are no 'spin component' lines.
## Do we have 'spin component' lines?
#compIxs = self.findLines( [compPat], 0, 0) # pats, numMin, numMax
#if len(compIxs) == 0: # if no spin components, or gw calc ...
# # if numSpin != 1: self.throwerr('numSpin mismatch', None)
#else: # else we have spin components ...
# if numSpin != 2: self.throwerr('numSpin mismatch', None)
# Usually DFT OUTCAR files have smallBandPat.
# Usually GW OUTCAR files have bigBandPat.
# But there are exceptions where GW files use smallBandPat.
smallBandPat = '^band No. +band energies +occupation$'
bigBandPat = '^band No. +\w+-energies +QP-energies +sigma\(\w+\) +V_xc\(\w+\) +V\^pw_x\(r,r\'\) +Z +occupation$'
smallBandIxs = self.findLines( [kpFirstPat, smallBandPat], 0, 0)
bigBandIxs = self.findLines( [kpFirstPat, bigBandPat], 0, 0)
if len(smallBandIxs) >= len( bigBandIxs):
# Probably DFT
firstIxs = smallBandIxs
eigCol = 1
occCol = 2
ntok = 3
hdrDelta = 3
else:
# Probably GW
firstIxs = bigBandIxs
eigCol = 2
occCol = 7
ntok = 8
hdrDelta = 4
if len(firstIxs) == 0:
resObj.numBand = None
resObj.eigenMat = None
resObj.occupMat = None
else:
for isp in range( numSpin):
istart = firstIxs[ -numSpin + isp] # use the last
isection = istart - 1 + hdrDelta
for ikp in range( numKpoint):
# Parse one k-point section, for numBand lines
tmat = self.parseMatrix(
ntok, # num tokens on each line
isection, 0, # rBeg, rEnd==0: go to blank line
0, ntok) # cBeg, cEnd
numEigen = tmat.shape[0]
isection += numEigen + hdrDelta
if eigenMat == None:
# Finally allocate matrices
eigenMat = np.zeros([ numSpin, numKpoint, numEigen])
occupMat = np.zeros([ numSpin, numKpoint, numEigen])
eigenMat[ isp, ikp, 0:numEigen] = tmat[ 0:numEigen, eigCol] # eig col
occupMat[ isp, ikp, 0:numEigen] = tmat[ 0:numEigen, occCol] # occ col
resObj.eigenMat = eigenMat
resObj.occupMat = occupMat
if self.bugLev >= 5:
for isp in range( resObj.numSpin):
print 'getEigenMat: isp: %d eigenMat: shape: %s\n%s' \
% (isp, resObj.eigenMat[isp].shape, resObj.eigenMat[isp],)
for isp in range( resObj.numSpin):
print 'getEigenMat: isp: %d occupMat: shape: %s\n%s' \
% (isp, resObj.occupMat[isp].shape, resObj.occupMat[isp],)
#====================================================================
# Get the line:
# number of dos NEDOS = 5001 number of ions NIONS = 2
#
# The imag part is the nedos lines following
#
# frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects)
# E(ev) X Y Z XY YZ ZX
# ----------------------------------------------------------------------
# 0.000000 6.024874 6.024874 6.024874 -0.000000 -0.000000 -0.000000
#
#
# The real part is the nedos lines following
#
# frequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)
# E(ev) X Y Z XY YZ ZX
# ----------------------------------------------------------------------
# 0.000000 6.024874 6.024874 6.024874 -0.000000 -0.000000 -0.000000
def getDielectric(
self, # self
resObj): # resulting object to be filled
'''
Extracts the dielectric function from OUTCAR, setting resObj attributes.
'''
patImag = '^ *frequency dependent *IMAGINARY DIELECTRIC FUNCTION'
patReal = '^ *frequency dependent *REAL DIELECTRIC FUNCTION'
ixs = self.findLines( [patImag], 0, 1) # numMin, numMax
if len(ixs) == 0:
resObj.dielectricImag = None
resObj.dielectricReal = None
if self.bugLev >= 5:
print 'getDielectric: all None'
else:
irow = ixs[0] + 3
resObj.dielectricImag = self.parseMatrix(
7, # num tokens on each line
irow, 0, # rBeg, rEnd==0: go to blank line
0, 7) # cBeg, cEnd
ixs = self.findLines( [patReal], 1, 1) # numMin, numMax
irow = ixs[0] + 3
resObj.dielectricReal = self.parseMatrix(
7, # num tokens on each line
irow, 0, # rBeg, rEnd==0: go to blank line
0, 7) # cBeg, cEnd
if self.bugLev >= 5:
print 'getDielectric: dielectricImag.shape: %s' \
% (resObj.dielectricImag.shape,)
print 'getDielectric: dielectricReal.shape: %s' \
% (resObj.dielectricReal.shape,)
print 'getDielectric: dielectricImag:\n%s\n' % (resObj.dielectricImag,)
print 'getDielectric: dielectricReal:\n%s\n' % (resObj.dielectricReal,)
#====================================================================
# Returns a list of line numbers ixs such that at each ix,
# lines[ix+i] matches pats[i].
[docs] def findLines(
self, # self
pats, # list of patterns to be matched by consecutive lines
numMin, # min num entire matches; else error
numMax): # max num entire matches or 0; else error
'''
Finds a list of line numbers in OUTCAR matching pats.
Each line num is the start of a group of lines
that consecutively match the patterns in pats.
'''
# Insure all pats are '^...$'
regexs = []
for pat in pats:
regexs.append( re.compile( pat))
resIxs = []
for iline in range( self.numLine):
# Do the pats match starting at iline ...
allOk = True
for ii in range( len( pats)):
if iline + ii < 0 or iline + ii >= self.numLine:
allOk = False
break
if not regexs[ii].match( self.lines[iline+ii]):
allOk = False
break
if allOk: resIxs.append( iline)
if len( resIxs) < numMin \
or (numMax > 0 and len( resIxs) > numMax):
self.throwerr(('num found mismatch. numMin: %d numMax: %d'
+ ' numFound: %d pats: %s')
% (numMin, numMax, len( resIxs), pats,), None)
return resIxs
#====================================================================
[docs] def parseMatrix(
self, # self
ntok, # number of tokens per line
begLine, # start line index
endLine, # 1 + last line index
begCol, # first column index, typically 0
endCol): # 1 + last column index, typically ntok
'''
Parses an array line::
0.0004124 2.1551750 2.1553184 -0.2320833 0.2320088 0.2320195
2.1551376 0.0004322 2.1552987 0.2320129 -0.2320836 0.2320237
2.1548955 2.1549131 0.0006743 0.2320650 0.2320652 -0.2320942
Returns an np.ndarray.
'''
rows = []
iline = begLine
while endLine == 0 or iline < endLine:
line = self.lines[iline]
if len(line) == 0 or re.match( '^-----+$', line):
if endLine != 0: self.throwerr('unexpected empty line', iline)
break
# Kluge:
# Fix run-together lines like: ... -2.252011135-14.597316840 ...
while True:
mat = re.search('\d-', line)
if not mat: break
ix = mat.start(0) + 1 # index of '-'
line = line[:ix] + ' ' + line[ix:]
toks = line.split()
if len(toks) != ntok:
self.throwerr(
'bad ntok. Expected: %d found: %d' % (ntok, len(toks,)), iline)
vals = map( float, toks[begCol:endCol])
rows.append( vals)
iline += 1
line = self.lines[iline]
if not (len(line) == 0 or re.match( '^-----+$', line)):
self.throwerr('bad endLine', endLine)
return np.array( rows, dtype=float)
#====================================================================
[docs] def checkClose(
self, # self
mata, # matrix A
matb, # matrix B
absErr): # max absolute error
'''
Checks that two matrices have the same shape and are close in value.
'''
if mata == None or matb == None:
self.throwerr('mata or matb is None', None)
if mata.shape != matb.shape:
self.throwerr('checkClose: shape mismatch: mata: %s matb: %s' \
% (mata.shape, matb.shape,))
if np.max( np.abs( mata - matb)) > absErr:
msg = 'checkClose: mismatch.\n'
msg += 'mata:\n%s\n' % (mata,)
msg += 'matb:\n%s\n' % (matb,)
msg += 'delta = matb - mata:\n%s\n' % (matb - mata,)
msg += 'max abs delta: %g\n' % (np.max( np.abs( matb - mata)))
msg += '\ncheckClose: mismatch. See above.'
self.throwerr( msg, None)
#====================================================================
[docs] def throwerr(
self, # self
msg, # error message
ix): # current line index (origin 0) in OUTCAR
'''
Raises an Exception
'''
fullMsg = '%s\n' % (msg,)
fullMsg += ' inDir: %s' % (self.inDir,)
if ix != None:
fullMsg += ' ix: %d\n line: %s\n' % (ix, self.lines[ix],)
raise Exception( fullMsg)
#====================================================================
# end of class ScanOutcar
#====================================================================
if __name__ == '__main__': main()