#!/usr/bin/env python
# Copyright 2013 National Renewable Energy Laboratory, Golden CO, USA
# This file is part of NREL MatDB.
#
# NREL MatDB 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.
#
# NREL MatDB 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 NREL MatDB. If not, see <http://www.gnu.org/licenses/>.
import datetime, re, sys, traceback, os.path
import xml.etree.cElementTree as etree
import numpy as np
np.seterr( all='raise', under='ignore')
import pylada.vasp # used by parsePylada tor parse OUTCAR files
import ScanXml # used to parse vasprun.xml files
import ScanOutcar # used to parse OUTCAR files
import wrapUpload # used for parseBoolean
#====================================================================
[docs]class ResClass:
'''
An empty class used as a container for parseDir results.
The parseDir function will call either parsePylada or parseXml,
and they will save the VASP results as attributes of
an instance of ResClass.
'''
def __str__(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
#====================================================================
def badparms( msg):
'''
Prints an error msg and usage info, and exits with rc 1.
'''
print '\nError: %s' % (msg,)
print 'Parms:'
print ' -bugLev <int> debug level'
print ' -readType <string> outcar / xml / compare'
print ' -allowExc <boolean> false/true: continue after error.'
print ' -inDir <string> dir containing input OUTCAR or vasprun.xml'
print ' -maxLev <int> max levels to print for xml'
print ''
print 'Examples:'
print './readVasp.py -bugLev 5 -readType xml -allowExc false -inDir tda/testlada.2013.04.15.fe.len.3.20/icsd_044729/icsd_044729.cif/hs-anti-ferro-0/relax_cellshape/0 -maxLev 0'
sys.exit(1)
#====================================================================
[docs]def main():
'''
Test driver: Extracts info from the output of a VASP run.
Command line parameters:
================ ========= ==============================================
Parameter Type Description
================ ========= ==============================================
**-bugLev** integer Debug level. Normally 0.
**-readType** string If 'outcar', read the OUTCAR file.
If 'xml', read the vasprun.xml file.
If 'compare', read both files and compare.
**-allowExc** boolean false/true: continue after error.
**-inDir** string Input directory containing OUTCAR
and/or vasprun.xml.
**-maxLev** int Max number of levels to print for xml
================ ========= ==============================================
'''
bugLev = 0
readType = None
allowExc = None
inDir = None
maxLev = 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 == '-readType': readType = val
elif key == '-allowExc': allowExc = wrapUpload.parseBoolean( val)
elif key == '-inDir': inDir = val
elif key == '-maxLev': maxLev = int( val)
else: badparms('unknown key: "%s"' % (key,))
if bugLev == None: badparms('parm not specified: -bugLev')
if readType == None: badparms('parm not specified: -readType')
if allowExc == None: badparms('parm not specified: -allowExc')
if inDir == None: badparms('parm not specified: -inDir')
if maxLev == None: badparms('parm not specified: -maxLev')
##np.set_printoptions( threshold=10000)
resObj = parseDir( bugLev, readType, allowExc, inDir, maxLev)
print '\nmain: resObj:\n%s' % (resObj,)
np.set_printoptions(threshold='nan')
np.set_printoptions(linewidth=80)
print '\n===== main: final resObj =====\n'
print ''
print 'import datetime'
print 'import numpy as np'
print ''
keys = resObj.__dict__.keys()
keys.sort()
msg = ''
for key in keys:
val = resObj.__dict__[key]
stg = repr(val)
if type(val).__name__ == 'ndarray':
print ''
print '# %s shape: %s' % (key, val.shape,)
stg = 'np.' + stg
print '%s = %s' % (key, stg,)
#====================================================================
# Returns ResClass instance.
# If all is ok, we return a ResClass instance with
# resObj.excMsg = None
# resObj.excTrace = None
# Else we return a ResClass instance with
# resObj.excMsg = exception message
# resObj.excTrace = exception traceback
[docs]def parseDir(
bugLev,
readType,
allowExc,
inDir,
maxLev):
'''
Extracts info from the output of a VASP run.
**Parameters**:
* bugLev (int): Debug level. Normally 0.
* readType (str):
If 'outcar', calls :class:`nrelmat.ScanOutcar` to
read the OUTCAR file.
If 'xml', calls :class:`nrelmat.ScanXml` to read the vasprun.xml file.
If 'compare', read both files and compare.
* allowExc (bool): continue after error.
* inDir (str): Input directory containing OUTCAR
and/or vasprun.xml.
* maxLev (int) Max number of levels to print for xml (normally 0)
**Returns**:
* resObj (class ResClass): data object with attributes set.
'''
resObj = ResClass()
resObj.excMsg = None
resObj.excTrace = None
try:
if not os.path.isdir(inDir):
throwerr('inDir is not a dir: "%s"' % (inDir,))
if readType == 'outcar':
scanner = ScanOutcar.ScanOutcar( bugLev, inDir, resObj) # fills resObj
#xxx del calcMisc( bugLev, inDir, resObj)
#xxx del calcBandgaps( bugLev, inDir, resObj)
elif readType == 'xml':
inFile = os.path.join( inDir, 'vasprun.xml')
if not os.path.isfile(inFile):
throwerr('inFile is not a file: "%s"' % (inFile,))
ScanXml.parseXml( bugLev, inFile, maxLev, resObj) # fills resObj
#xxx del calcMisc( bugLev, inDir, resObj)
#xxx del calcBandgaps( bugLev, inDir, resObj)
elif readType == 'compare':
resOut = ResClass()
scanner = ScanOutcar.ScanOutcar( bugLev, inDir, resOut) # fills resOut
#xxx del calcMisc( bugLev, inDir, resOut)
#xxx del calcBandgaps( bugLev, inDir, resOut)
resXml = ResClass()
inFile = os.path.join( inDir, 'vasprun.xml')
if not os.path.isfile(inFile):
throwerr('inFile is not a file: "%s"' % (inFile,))
ScanXml.parseXml( bugLev, inFile, maxLev, resXml) # fills resXml
#xxx del calcMisc( bugLev, inDir, resXml)
#xxx del calcBandgaps( bugLev, inDir, resXml)
cmsg = wrapUpload.deepCompare(
bugLev, 1.e-15, 'outcar', resOut, 'xml', resXml)
if cmsg == None:
print 'readVasp: deepCompare says identical'
else:
print '\n===== readVasp: begin deepCompare diff msg =====\n'
print '%s' % (cmsg,)
print '\n===== readVasp: end deepCompare diff msg =====\n'
# The pylada scanner is obsolete: it doesn't handle
# some files, and has bugs.
#elif readType == 'pylada':
# parsePyladaUnused( bugLev, inFile, resObj) # fills resObj
else: throwerr('unknown readType: %s' % (readType,))
except Exception, exc:
resObj.excMsg = repr(exc)
resObj.excTrace = traceback.format_exc( limit=None)
if bugLev >= 1:
print 'readVasp.py. caught exc: %s' % (resObj.excMsg,)
print ' readType: "%s"' % (readType,)
print ' inDir: "%s"' % (inDir,)
print '===== traceback start ====='
print resObj.excTrace
print '===== traceback end ====='
if not allowExc: throwerr('caught: %s' % (exc,))
return resObj
#====================================================================
# Too unreliable; too many types of VASP runs
#def setRunType( bugLev, resObj):
# '''
# Sets resObj.runType = 'typeGw' or 'typeChi' or 'typeRelax',
# depending on resObj.algo.
# Called by ScanOutcar and ScanXml (a callback).
# '''
#
# # xxx if not resObj.__dict__.has_key('algo'):
# # xxx resObj.__dict__['algo'] = 'Normal'
#
# algo = getattr( resObj, 'algo', 'normal') # if not found, assume 'normal'
# algo = algo.lower()
# if algo == 'a': algo = 'all'
#
# # See http://cms.mpi.univie.ac.at/wiki/index.php/ALGO
# if algo in ['gw', 'gw0', 'scgw', 'scgw0']:
# resObj.runType = 'typeGw'
# elif algo in ['chi']:
# resObj.runType = 'typeChi'
# elif algo in [ 'all', 'conjugate',
# 'damp', 'damped', 'diag',
# 'eigenval', 'fast', 'none', 'norm', 'normal',
# 'nothing', 'subrot', 'veryfast', 'n/a']:
# resObj.runType = 'typeRelax'
# else: throwerr('unknown algo: %s' % (algo,))
# if bugLev >= 5: print 'algo: %s runType: %s' % (algo, resObj.runType,)
#====================================================================
# xxx del def calcMisc(
# xxx del bugLev,
# xxx del inDir,
# xxx del resObj):
# xxx del '''
# xxx del Sets resObj.atomTypes, atomNames, etc, based on
# xxx del resObj.typeNames, etc. The atom* attributes are one per atom;
# xxx del the type* attributes are one per unique element. So
# xxx del H2O has 3 atomNames ['H', 'H', 'O'], and 2 typeNames ['H', 'O']
# xxx del with typeNums [2, 1].
# xxx del '''
# xxx del
# xxx del if resObj.iterRealTimes == None: resObj.iterTotalTime = None
# xxx del else: resObj.iterTotalTime = sum( resObj.iterRealTimes)
# xxx del if bugLev >= 5:
# xxx del print 'calcMisc: iterTotalTime: %s' % (resObj.iterTotalTime,)
# xxx del
# xxx del numType = len( resObj.typeNames)
# xxx del if len( resObj.typeNums) != numType: throwerr('numType mismatch')
# xxx del
# xxx del atomTypes = [] # indices into typeNames, typeNums
# xxx del atomNames = []
# xxx del atomMasses_amu = []
# xxx del atomPseudos = []
# xxx del atomValences = []
# xxx del for ii in range( numType):
# xxx del numEle = resObj.typeNums[ii]
# xxx del atomTypes += numEle * [ ii ]
# xxx del atomNames += numEle * [ resObj.typeNames[ii] ]
# xxx del atomMasses_amu += numEle * [ resObj.typeMasses_amu[ii] ]
# xxx del atomPseudos += numEle * [ resObj.typePseudos[ii] ]
# xxx del atomValences += numEle * [ resObj.typeValences[ii] ]
# xxx del if bugLev >= 5:
# xxx del print 'calcMisc: atomTypes: %s' % (atomTypes,)
# xxx del print 'calcMisc: atomNames: %s' % (atomNames,)
# xxx del print 'calcMisc: atomMasses_amu: %s' % (atomMasses_amu,)
# xxx del print 'calcMisc: atomPseudos: %s' % (atomPseudos,)
# xxx del print 'calcMisc: atomValences: %s' % (atomValences,)
# xxx del
# xxx del if getattr( resObj, 'atomTypes', None) == None:
# xxx del resObj.atomTypes = atomTypes
# xxx del else:
# xxx del if resObj.atomTypes != atomTypes:
# xxx del throwerr('atomTypes mismatch:\n old: %s\n new: %s\n'
# xxx del % (resObj.atomTypes, atomTypes,))
# xxx del
# xxx del if getattr( resObj, 'atomNames', None) == None:
# xxx del resObj.atomNames = atomNames
# xxx del else:
# xxx del if resObj.atomNames != atomNames:
# xxx del throwerr('atomNames mismatch:\n old: %s\n new: %s\n'
# xxx del % (resObj.atomNames, atomNames,))
# xxx del
# xxx del if getattr( resObj, 'atomMasses_amu', None) == None:
# xxx del resObj.atomMasses_amu = atomMasses_amu
# xxx del else:
# xxx del if resObj.atomMasses_amu != atomMasses_amu:
# xxx del throwerr('atomMasses_amu mismatch:\n old: %s\n new: %s\n'
# xxx del % (resObj.atomMasses_amu, atomMasses_amu,))
# xxx del
# xxx del if getattr( resObj, 'atomPseudos', None) == None:
# xxx del resObj.atomPseudos = atomPseudos
# xxx del else:
# xxx del if resObj.atomPseudos != atomPseudos:
# xxx del throwerr('atomPseudos mismatch:\n old: %s\n new: %s\n'
# xxx del % (resObj.atomPseudos, atomPseudos,))
# xxx del
# xxx del if getattr( resObj, 'atomValences', None) == None:
# xxx del resObj.atomValences = atomValences
# xxx del else:
# xxx del if resObj.atomValences != atomValences:
# xxx del throwerr('atomValences mismatch:\n old: %s\n new: %s\n'
# xxx del % (resObj.atomValences, atomValences,))
# xxx del
# xxx del
# xxx del
# xxx del
# xxx del # totalValence = sum( count[i] * valence[i])
# xxx del # PyLada calls this valence.
# xxx del resObj.totalValence = np.dot( resObj.typeNums, resObj.typeValences)
# xxx del if bugLev >= 5:
# xxx del print 'calcMisc: totalValence: %s' % (resObj.totalValence,)
# xxx del
# xxx del
# xxx del
# xxx del # In defect calcs of Lany, we can have numElectron==288
# xxx del # and totalValence==289, leaving the supercell with
# xxx del # positive charge, for a donor.
# xxx del resObj.netCharge = resObj.totalValence - resObj.numElectron
# xxx del
# xxx del
# xxx del
# xxx del # volume, Angstrom^3
# xxx del # The scale is hard coded as 1.0 in PyLada crystal/read.py,
# xxx del # in both icsd_cif_a and icsd_cif_b.
# xxx del if getattr( resObj, 'finalBasisMat', None) == None:
# xxx del resObj.finalVolumeCalc_ang3 = None
# xxx del resObj.finalDensity_g_cm3 = None
# xxx del else:
# xxx del volScale = 1.0
# xxx del resObj.finalVolumeCalc_ang3 = abs( np.linalg.det(
# xxx del volScale * resObj.finalBasisMat))
# xxx del # Density
# xxx del # xxx better: get atomic weights from periodic table
# xxx del volCm3 = resObj.finalVolumeCalc_ang3 / (1.e8)**3 # 10**8 Angstrom per cm
# xxx del totMass = np.dot( resObj.typeNums, resObj.typeMasses_amu)
# xxx del totMassGm = totMass * 1.660538921e-24 # 1.660538921e-24 g / amu
# xxx del resObj.finalDensity_g_cm3 = totMassGm / volCm3
# xxx del if bugLev >= 5:
# xxx del print 'finalVolumeCalc_ang3: %s' % (resObj.finalVolumeCalc_ang3,)
# xxx del print 'volCm3: %g' % (volCm3,)
# xxx del print 'totMass: %g' % (totMass,)
# xxx del print 'totMassGm: %g' % (totMassGm,)
# xxx del print 'finalDensity_g_cm3: %g' % (resObj.finalDensity_g_cm3,)
# xxx del
# xxx del
# xxx del
# xxx del # reciprocal space volume, * (2*pi)**3
# xxx del # As in PyLada.
# xxx del if getattr( resObj, 'finalBasisMat', None) == None:
# xxx del resObj.recipVolume = None
# xxx del else:
# xxx del invMat = np.linalg.inv( volScale * resObj.finalBasisMat)
# xxx del resObj.recipVolume = abs( np.linalg.det( invMat)) * (2 * np.pi)**3
# xxx del if bugLev >= 5:
# xxx del print 'recipVolume: origMat:\n%s' \
# xxx del % (repr(volScale * resObj.finalBasisMat),)
# xxx del print 'recipVolume: invMat:\n%s' % (repr(invMat),)
# xxx del print 'recipVolume: det:\n%s' % (repr(np.linalg.det( invMat)),)
# xxx del print 'recipVolume: %g' % (resObj.recipVolume,)
# xxx del
# xxx del
# xxx del # Calc pressure
# xxx del # xxx Caution: we do not include the non-diag terms in:
# xxx del # VASP: force.F: FORCE_AND_STRESS: line 1410:
# xxx del # PRESS=(TSIF(1,1)+TSIF(2,2)+TSIF(3,3))/3._q &
# xxx del # & -DYN%PSTRESS/(EVTOJ*1E22_q)*LATT_CUR%OMEGA
# xxx del
# xxx del if getattr( resObj, 'finalStressMat_kbar', None) == None:
# xxx del resObj.finalPressure_kbar = None
# xxx del else:
# xxx del resObj.finalPressure_kbar \
# xxx del = resObj.finalStressMat_kbar.diagonal().sum() / 3.0
# xxx del if bugLev >= 5:
# xxx del print 'getStressMat: finalPressure_kbar: %s' \
# xxx del % (resObj.finalPressure_kbar,)
# xxx del
# xxx del
# xxx del #====================================================================
# xxx del
# xxx del def calcBandgaps(
# xxx del bugLev,
# xxx del inDir,
# xxx del resObj):
# xxx del '''
# xxx del Using resObj.eigenMat sets resObj.cbMinVals, cbMinIxs, cbMin,
# xxx del and similarly for Max,
# xxx del and bandgapDirects, bandgapIndirects, and bandgap.
# xxx del '''
# xxx del
# xxx del if resObj.eigenMat == None:
# xxx del resObj.cbMinVals = None
# xxx del resObj.cbMinIxs = None
# xxx del resObj.vbMaxVals = None
# xxx del resObj.vbMaxIxs = None
# xxx del resObj.bandgapDirects = None
# xxx del resObj.bandgapIndirects = None
# xxx del resObj.bandgapDirect = None
# xxx del resObj.bandgapIndirect = None
# xxx del
# xxx del else:
# xxx del # Conduction band min
# xxx del # cbMinMat[isp][ikp] = min of eigvals > efermi0
# xxx del #
# xxx del # Valence band max
# xxx del # vbMaxMat[isp][ikp] = max of eigvals <= efermi0
# xxx del
# xxx del print 'calcBandGaps: eig shape: ', resObj.eigenMat.shape
# xxx del print 'calcBandGaps: numSpin: %d numKpoint: %d numBand: %d' \
# xxx del % (resObj.numSpin, resObj.numKpoint, resObj.numBand,)
# xxx del
# xxx del numSpin = resObj.numSpin
# xxx del numKpoint = resObj.numKpoint
# xxx del # Sometimes numBand == 128 but eigenMat.shape == [x,x,64]
# xxx del #xxx numBand = resObj.numBand
# xxx del numBand = resObj.eigenMat.shape[2]
# xxx del
# xxx del # For each isp and kpoint, what is min eigvalue > efermi0
# xxx del cbMinMat = np.empty( [numSpin, numKpoint], dtype=float)
# xxx del vbMaxMat = np.empty( [numSpin, numKpoint], dtype=float)
# xxx del cbMinMat.fill( np.inf)
# xxx del vbMaxMat.fill( -np.inf)
# xxx del
# xxx del # For each isp, what is min or max across all kpoints
# xxx del cbMinVals = np.empty( [numSpin], dtype=float)
# xxx del vbMaxVals = np.empty( [numSpin], dtype=float)
# xxx del cbMinVals.fill( np.inf)
# xxx del vbMaxVals.fill( -np.inf)
# xxx del
# xxx del # For each isp, which kpoint has the min or max
# xxx del cbMinIxs = np.empty( [numSpin], dtype=int)
# xxx del vbMaxIxs = np.empty( [numSpin], dtype=int)
# xxx del cbMinIxs.fill( -1)
# xxx del vbMaxIxs.fill( -1)
# xxx del
# xxx del for isp in range( numSpin):
# xxx del for ikp in range( numKpoint):
# xxx del for iband in range( numBand):
# xxx del val = resObj.eigenMat[isp][ikp][iband]
# xxx del
# xxx del algorithm = 'occupMat'
# xxx del
# xxx del if algorithm == 'fermi0':
# xxx del if val <= resObj.efermi0:
# xxx del if val > vbMaxMat[isp][ikp]: vbMaxMat[isp][ikp] = val
# xxx del if val > vbMaxVals[isp]:
# xxx del vbMaxVals[isp] = val
# xxx del vbMaxIxs[isp] = ikp
# xxx del else: # val > resObj.efermi0
# xxx del if val < cbMinMat[isp][ikp]: cbMinMat[isp][ikp] = val
# xxx del if val < cbMinVals[isp]:
# xxx del cbMinVals[isp] = val
# xxx del cbMinIxs[isp] = ikp
# xxx del
# xxx del elif algorithm == 'occupMat':
# xxx del if numSpin == 1:
# xxx del occLo = 0.02
# xxx del occHi = 1.98
# xxx del elif numSpin == 2:
# xxx del occLo = 0.01
# xxx del occHi = 0.99
# xxx del else: throwerr('unknown numSpin')
# xxx del
# xxx del occ = resObj.occupMat[isp][ikp][iband]
# xxx del if occ >= occHi:
# xxx del if val > vbMaxMat[isp][ikp]: vbMaxMat[isp][ikp] = val
# xxx del if val > vbMaxVals[isp]:
# xxx del vbMaxVals[isp] = val
# xxx del vbMaxIxs[isp] = ikp
# xxx del if occ <= occLo:
# xxx del if val < cbMinMat[isp][ikp]: cbMinMat[isp][ikp] = val
# xxx del if val < cbMinVals[isp]:
# xxx del cbMinVals[isp] = val
# xxx del cbMinIxs[isp] = ikp
# xxx del
# xxx del else: throwerr('unknown algorithm')
# xxx del
# xxx del # Find bandgapDirects[isp] = min gap for the same kpoint
# xxx del # Find bandgapIndirects[isp] = min gap across kpoints
# xxx del bandgapDirects = numSpin * [np.inf]
# xxx del bandgapIndirects = numSpin * [np.inf]
# xxx del for isp in range( numSpin):
# xxx del for ikp in range( numKpoint):
# xxx del gap = max( 0, cbMinMat[isp][ikp] - vbMaxMat[isp][ikp])
# xxx del if gap < bandgapDirects[isp]: bandgapDirects[isp] = gap
# xxx del bandgapIndirects[isp] = max( 0, cbMinVals[isp] - vbMaxVals[isp])
# xxx del
# xxx del resObj.cbMinVals = cbMinVals
# xxx del resObj.cbMinIxs = cbMinIxs
# xxx del resObj.vbMaxVals = vbMaxVals
# xxx del resObj.vbMaxIxs = vbMaxIxs
# xxx del resObj.bandgapDirects = bandgapDirects
# xxx del resObj.bandgapIndirects = bandgapIndirects
# xxx del resObj.bandgapDirect = min( bandgapDirects)
# xxx del resObj.bandgapIndirect = min( bandgapIndirects)
# xxx del
# xxx del # This is the PyLada version
# xxx del #resObj.cbMin = min( cbMinVals)
# xxx del #resObj.vbMax = max( vbMaxVals)
# xxx del #resObj.bandgap = resObj.cbMin - resObj.vbMax
# xxx del #if resObj.bandgap < 0: resObj.bandgap = 0
# xxx del
# xxx del if bugLev >= 5:
# xxx del
# xxx del print '\ncalcBandgaps: efermi0: %s' % (resObj.efermi0,)
# xxx del for isp in range( numSpin):
# xxx del print '\ncalcBandgaps: start isp: %d' % (isp,)
# xxx del for ikp in range( numKpoint):
# xxx del vval = vbMaxMat[isp][ikp]
# xxx del cval = cbMinMat[isp][ikp]
# xxx del vmsg = '%8.4f' % (vval,)
# xxx del cmsg = '%8.4f' % (cval,)
# xxx del if ikp == vbMaxIxs[isp]: vmsg += ' ***'
# xxx del if ikp == cbMinIxs[isp]: cmsg += ' ***'
# xxx del print ' isp: %d ikp: %3d vbMaxMat: %-14s cbMinMat: %-14s' \
# xxx del % ( isp, ikp, vmsg, cmsg,)
# xxx del print ''
# xxx del
# xxx del print 'calcBandgaps: vbMaxIxs: %s' % (resObj.vbMaxIxs,)
# xxx del print 'calcBandgaps: vbMaxVals: %s' % (resObj.vbMaxVals,)
# xxx del print 'calcBandgaps: cbMinIxs: %s' % (resObj.cbMinIxs,)
# xxx del print 'calcBandgaps: cbMinVals: %s' % (resObj.cbMinVals,)
# xxx del
# xxx del print 'calcBandgaps: bandgapDirects: %s' % (resObj.bandgapDirects,)
# xxx del print 'calcBandgaps: bandgapIndirects: %s' % (resObj.bandgapIndirects,)
# xxx del print 'calcBandgaps: bandgapDirect: %s' % (resObj.bandgapDirect,)
# xxx del print 'calcBandgaps: bandgapIndirect: %s' % (resObj.bandgapIndirect,)
# xxx del
# xxx del # Write plot data to tempplota%ispin
# xxx del for isp in range( numSpin):
# xxx del with open('tempplota%d' % (isp,), 'w') as fout:
# xxx del for ikp in range( numKpoint):
# xxx del print >> fout, '%g %g' % (vbMaxMat[isp,ikp], cbMinMat[isp,ikp],)
# xxx del
# xxx del # Write plot data to tempplotb%ispin
# xxx del for isp in range( numSpin):
# xxx del with open('tempplotb%d' % (isp,), 'w') as fout:
# xxx del for ikp in range( numKpoint):
# xxx del for iband in range( numBand):
# xxx del print >> fout, '%d %g' \
# xxx del % (ikp, resObj.eigenMat[isp][ikp][iband],)
# xxx del
# xxx del
# xxx del # Obsolete way to calc bandgaps
# xxx del ###if bugLev >= 5: print '\n===== cbMin, vbMax, bandgap =====\n'
# xxx del
# xxx del #### Find cbm = min of eigs > efermi0
# xxx del #### Find vbm = max of eigs <= efermi0
# xxx del
# xxx del ###cbms = resObj.numSpin * [np.inf]
# xxx del ###vbms = resObj.numSpin * [-np.inf]
# xxx del ###cbmKpis = resObj.numSpin * [None]
# xxx del ###vbmKpis = resObj.numSpin * [None]
# xxx del
# xxx del ###for isp in range( resObj.numSpin):
# xxx del ### eigs = resObj.eigenMat[isp]
# xxx del ### for ikp in range( resObj.numKpoint):
# xxx del ### for iband in range( resObj.numBand):
# xxx del ### val = eigs[ikp][iband]
# xxx del ### if val > resObj.efermi0:
# xxx del ### cbms[isp] = min( cbms[isp], val)
# xxx del ### cbmKpis[isp] = ikp
# xxx del ### if val <= resObj.efermi0:
# xxx del ### vbms[isp] = max( vbms[isp], val)
# xxx del ### vbmKpis[isp] = ikp
# xxx del
# xxx del ###cbms = map( float, cbms) # change type from numpy.float64 to float
# xxx del ###vbms = map( float, vbms) # change type from numpy.float64 to float
# xxx del
# xxx del ###resObj.cbms = cbms
# xxx del ###resObj.vbms = vbms
# xxx del ###resObj.cbmKpis = cbmKpis
# xxx del ###resObj.vbmKpis = vbmKpis
# xxx del ###resObj.cbMin = min( cbms) # This is PyLada's cbm
# xxx del ###resObj.vbMax = max( vbms) # This is PyLada's vbm
# xxx del
# xxx del ###resObj.bandgaps = [ (cbms[ii] - vbms[ii]) for ii in range( resObj.numSpin)]
# xxx del ###resObj.bandgapa = min( resObj.bandgaps)
# xxx del ###resObj.bandgap = resObj.cbMin - resObj.vbMax # This is PyLada version
# xxx del
# xxx del ###if bugLev >= 5:
# xxx del ### print 'cbmKpis: %s cbms: %s' % (cbmKpis, cbms,)
# xxx del ### print 'vbmKpis: %s vbms: %s' % (vbmKpis, vbms,)
# xxx del ### print 'cbMin: %g' % (resObj.cbMin,)
# xxx del ### print 'vbMax: %g' % (resObj.vbMax,)
# xxx del ### print 'bandgaps: %s' % (resObj.bandgaps,)
# xxx del ### print 'bandgapa: %g' % (resObj.bandgapa,)
# xxx del ### print 'bandgap: %g' % (resObj.bandgap,)
#====================================================================
# No longer used:
# Calc the efermi0 by adding up occupancies,
# as in PyLada.
def calcEfermi0Unused( bugLev, resObj):
if getattr( resObj, 'eigenMat', None) == None:
resObj.efermi0Calc = None
else:
# efermi0Calc
# Make an array allEigs of all the eigenvalues for all the
# kpoints, with each eigenvalue replicated by the associated
# kpoint weight.
# Further, if ispin == 1 (non spin polarized),
# replicate each eigenvalue again.
#
# Sort allEigs by increasing value.
# Let factor =
# if numSpin == 1: 2/sum(weights) # non spin polarized
# if numSpin == 2: 1/sum(weights) # spin polarized
# Start summing: occ = factor * (index into allEigs).
# When occ = totalValence, the corresponding eigVal is fermi0K.
numSpin = resObj.numSpin
numKpoint = resObj.numKpoint
numBand = resObj.numBand
# xxx check all this
numEig = 0
for ikp in range( numKpoint):
numEig += numBand * resObj.kpointMults[ikp]
numEig *= 2 # for two spin levels, whether numSpin is 1 or 2
numEig = int( round( numEig))
if bugLev >= 5:
print 'calcEfermi0: numSpin: %d numKpoint: %d numBand: %d' \
% (numSpin, numKpoint, numBand,)
print 'calcEfermi0: numEig: %d' % (numEig,)
avgWt = sum( resObj.kpointMults) / float( numKpoint)
allEigs = np.zeros( [numEig])
kk = 0
for isp in range( numSpin):
for ikp in range( numKpoint):
for iband in range( numBand):
# Replicate by kpoint weight
for irepl in range( int( round( resObj.kpointMults[ ikp]))):
allEigs[kk] = resObj.eigenMat[ isp, ikp, iband]
kk += 1
if numSpin == 1:
# Replicate for both spins
allEigs[kk] = allEigs[kk-1]
kk += 1
if kk != numEig: throwerr('numEig mismatch')
allEigs.sort()
if bugLev >= 5:
print 'calcEfermi0: numEig: %d' % (numEig,)
print 'calcEfermi0: allEigs.shape: %s' % (allEigs.shape,)
# Find indx such that
# occupancy = indx / sum(kpWts) == valence
# Set dindx = num eigenvalues used.
# The array index of the last used eigenvalue is dindx-1.
dindx = resObj.totalValence * sum( resObj.kpointMults)
indx = int( round( dindx)) - 1
if indx < 0 or indx >= allEigs.shape[0]:
throwerr('bad fermi. dindx: %g indx: %d shape: %s' \
% (dindx, indx, allEigs.shape))
resObj.efermi0Calc = allEigs[indx]
if bugLev >= 5:
print 'calcEfermi0: totalValence: %g' % (resObj.totalValence,)
print 'calcEfermi0: sum( resObj.kpointMults): %g' % (sum( resObj.kpointMults),)
print 'calcEfermi0: len(allEigs): %g' % (len(allEigs),)
print 'calcEfermi0: dindx: %g indx: %d' % (dindx, indx,)
print 'calcEfermi0: efermi0: %g' % (resObj.efermi0,)
print 'calcEfermi0: efermi0Calc: %g' % (resObj.efermi0Calc,)
for ii in range( max( 0, indx-10), min( numEig, indx+10)):
msg = 'calcEfermi0: allEigs[%d]: %g' % (ii, allEigs[ii],)
if ii == indx: msg += ' ***'
print msg
#====================================================================
#====================================================================
# The pylada scanner is obsolete: it doesn't handle
# some files, and has bugs.
def parsePyladaUnused( bugLev, inFile, resObj):
'''
Extracts info from the OUTCAR file from a VASP run,
using the PyLada vasp.Extract API.
**Parameters**:
* bugLev (int): Debug level. Normally 0.
* inFile (str): Path of the input OUTCAR file.
* resObj (class ResClass): data object: we set attributes here.
**Returns**:
* None
'''
ex = pylada.vasp.Extract( inFile)
if not ex.success: throwerr('file %s is not complete' % (inFile,))
#===== program, version, date etc =====
resObj.runDate = ex.datetime
# iterTimes are pairs: [cpuTime, realTime]
resObj.iterCpuTimes = [pair[0] for pair in ex.iterTimes]
resObj.iterRealTimes = [pair[1] for pair in ex.iterTimes]
resObj.iterTotalTime = np.sum( resObj.iterRealTimes)
#===== incar parameters =====
# Use float() to get rid of Quantities units
resObj.algo = ex.algo
resObj.ediff = float( ex.ediff)
resObj.encut_ev = float( ex.encut) # eV
resObj.ibrion = ex.ibrion
resObj.isif = ex.isif
resObj.systemName = ex.system
#===== general parameters =====
#===== electronic parameters =====
resObj.ialgo = ex.ialgo
resObj.numBand = ex.nbands
resObj.numElectron = ex.nelect
resObj.icharg = ex.icharg
resObj.numSpin = ex.ispin
#===== atom info =====
resObj.typeNames = np.array( ex.species, dtype=str)
resObj.typeNums = np.array( ex.stoichiometry, dtype=int)
resObj.typeValences = np.array( ex.ionic_charges, dtype=float)
# totalValence = sum( count[i] * valence[i])
# PyLada calls this valence.
resObj.totalValence = np.dot( resObj.typeNums, resObj.typeValences)
if bugLev >= 5: print 'totalValence: %g' % (resObj.totalValence,)
if resObj.numElectron != resObj.totalValence:
# xxx should this be an error?
print('%g == numElectron != totalValence == %g' \
% (resObj.numElectron, resObj.totalValence,))
struct = ex.structure
natom = len( struct)
resObj.atomNames = [ struct[ii].type for ii in range( natom)]
resObj.atomMasses_amu = None # not available
resObj.atomPseudos = None # not available
resObj.atomValences = None # not available
# atomPseudos: only the first title is available
# atomMass, Valence: not available
# partial_charges and magnetization are not available in xml
# resObj.partialChargeMat = ex.partial_charges
# resObj.magnetizationMat = ex.magnetization
#===== initial structure =====
# In structure.cell, the basis vectors are columns.
# Change to rows.
resObj.initialBasisMat = ex.initial_structure.cell.T
resObj.initialRecipBasisMat = None # not available
struct = ex.initial_structure
natom = len( struct)
cartPos = natom * [None]
for ii in range( natom):
cartPos[ii] = struct[ii].pos # no units
resObj.initialCartPosMat = np.array( cartPos)
resObj.initialFracPosMat = np.dot(
resObj.initialCartPosMat, np.linalg.inv( resObj.initialBasisMat))
#===== final structure =====
# In structure.cell, the basis vectors are columns.
# Change to rows.
resObj.finalBasisMat = ex.structure.cell.T
struct = ex.structure
natom = len( struct)
cartPos = natom * [None]
for ii in range( natom):
cartPos[ii] = struct[ii].pos
resObj.finalCartPosMat = np.array( cartPos)
resObj.finalFracPosMat = np.dot(
resObj.finalCartPosMat, np.linalg.inv( resObj.finalBasisMat))
#===== kpoints =====
resObj.kpointCartMat = ex.kpoints # transform
resObj.numKpoint = resObj.kpointCartMat.shape[0]
resObj.kpointWeights = None # xxx not available
#resObj.kpointFracMat = np.dot(
# resObj.kpointCartMat, np.linalg.inv( resObj.initialRecipBasisMat))
# initialRecipBasisMat is not available, so
resObj.kpointFracMat = None
if bugLev >= 5:
print 'numKpoint: %g' % (resObj.numKpoint,)
print 'kpointFracMat:\n%s' % (repr(resObj.kpointFracMat),)
print 'kpointCartMat:\n%s' % (repr(resObj.kpointCartMat),)
#===== final volume and density =====
resObj.finalVolume_ang3 = float( ex.volume) # Angstrom^3
resObj.recipVolume = float( ex.reciprocal_volume)
resObj.finalDensity_g_cm3 = float( ex.density) # g/cm3
#===== last calc forces =====
resObj.finalForceMat_ev_ang = np.array( ex.forces) # eV/angstrom
resObj.finalStressMat_kbar = np.array( ex.stress) # kbar
resObj.finalPressure_kbar = float( ex.pressure) # kbar
#===== eigenvalues and occupancies =====
resObj.eigenMat = np.array( ex.eigenvalues)
# Caution: for non-magnetic (numSpin==1),
# OUTCAR has occupMat values = 2, while vasprun.xml has values = 1.
# For magnetic (numSpin==2), both OUTCAR and vasprun.xml have 1.
resObj.occupMat = np.array( ex.occupations)
#===== misc junk =====
#===== energy, efermi0 =====
resObj.energyNoEntrps = [float( ex.total_energy)] # eV
resObj.efermi0 = float( ex.fermi0K) # eV
#===== cbMin, vbMax, bandgap =====
resObj.cbMin = float( ex.cbm) # eV
resObj.vbMax = float( ex.vbm) # eV
resObj.bandgap = resObj.cbMin - resObj.vbMax # eV
if bugLev >= 5:
print 'cbMin: %g' % (resObj.cbMin,)
print 'vbMax: %g' % (resObj.vbMax,)
print 'bandgap: %g' % (resObj.bandgap,)
resObj.cbms = None # not available
resObj.vbms = None # not available
resObj.cbmKpis = None # not available
resObj.vbmKpis = None # not available
resObj.bandgapa = None # not available
resObj.bandgaps = None # not available
return
# xxx to do:
# energy == totalEnergy
# energy_sigma0 totalEnergySigma0
# fermi0K fermi0
# fermi_energy fermiEnergy
# total_energy totalEnergy
# array(-13.916343) * eV
#====================================================================
def throwerr( msg):
'''
Prints an error message and raises Exception.
**Parameters**:
* msg (str): Error message.
**Returns**
* (Never returns)
**Raises**
* Exception
'''
# print msg
# print >> sys.stderr, msg
raise Exception( msg)
#====================================================================
if __name__ == '__main__': main()
#====================================================================