#!/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, fractions, json, os, re, sys
import psycopg2
import numpy as np
from pyspglib import spglib
import wrapUpload
import fillDbVasp
#====================================================================
class DbObj( object):
'''
Empty class so we can set attributes
'''
pass
#====================================================================
# Creates an ASE bulk look-alike object, as required by spglib,
# using the info retrieved from the DB.
class AseBulk:
def __init__( self, rowCellMat, cartPosMat, atomNames):
self.rowCellMat = rowCellMat # cell matrix, in rows
self.cartPosMat = cartPosMat # cartesian position matrix
self.atomNames = atomNames # names of atoms, like
# ['Fe', 'Fe', 'O', 'O', 'O']
def get_cell( self):
# The spglib python wants rows.
# Internally the spglib python converts to columns
# for the underlying spglib C library.
return np.array( self.rowCellMat, dtype='double', order='C')
def get_scaled_positions( self):
cellInv = np.linalg.inv( self.rowCellMat)
scalePosMat = np.array([np.dot( row, cellInv ) for row in self.cartPosMat])
return np.array( scalePosMat, dtype='double', order='C')
def get_atomic_numbers( self):
from pylada import periodic_table
atomNums = [periodic_table.symbols.index(nm)+1 for nm in self.atomNames]
return atomNums
#====================================================================
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 ' -useCommit <boolean> false/true: do we commit changes to the DB.'
print ' -inSpec <string> inSpecJsonFile'
sys.exit(1)
#====================================================================
[docs]def main():
'''
This program adds additional information to the model database table.
The following functions in this file fill the indicated columns
in the model table. The columns must have been created previously
by the fillDbVasp.py function createTable.
============== ============= ==============================================
Function Column Notes
============== ============= ==============================================
addChemform formula Standard chemical formula: ``"H2 O"``.
addChemform sortedformula Chemical formula: ``"H2 O"``, with
atoms sorted in alphabetic order.
addChemForm chemtext Structured formula, easier for parsing.
Every token is surrounded by spaces, and
single occurance atoms have the explicit
"1" count:
``" H 2 O 1 "``.
addMinenergy minenergyid == mident having min energyperatom
for this formula and finalspacegroupnum.
============== ============= ==============================================
Command line parameters:
============== =========== ===============================================
Parameter Type Description
============== =========== ===============================================
**-bugLev** integer Debug level. Normally 0.
**-useCommit** boolean false/true: do we commit changes to the DB.
**-inSpec** string JSON file containing DB parameters. See below.
============== =========== ===============================================
**inSpec File Parameters:**
=================== ==============================================
Parameter Description
=================== ==============================================
**dbhost** Database hostname.
**dbport** Database port number.
**dbuser** Database user name.
**dbpswd** Database password.
**dbname** Database database name.
**dbschema** Database schema name.
**dbtablemodel** Database name of the "model" table.
=================== ==============================================
**inSpec file example:**::
{
"dbhost" : "scctest",
"dbport" : "6432",
"dbuser" : "x",
"dbpswd" : "x",
"dbname" : "cidlada",
"dbschema" : "satom",
"dbtablemodel" : "model"
}
'''
bugLev = None
useCommit = None
inSpec = 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 == '-useCommit': useCommit = wrapUpload.parseBoolean( val)
elif key == '-inSpec': inSpec = val
else: badparms('unknown key: "%s"' % (key,))
if bugLev == None: badparms('parm not specified: -bugLev')
if useCommit == None: badparms('parm not specified: -useCommit')
if inSpec == None: badparms('parm not specified: -inSpec')
# Register psycopg2 adapter for np.ndarray
# The function formatArray is defined below.
def adaptVal( val):
msg = fillDbVasp.formatArray( val)
# AsIs provides getquoted() which just calls the wrapped object's str().
return psycopg2.extensions.AsIs( msg)
psycopg2.extensions.register_adapter( np.ndarray, adaptVal)
psycopg2.extensions.register_adapter( np.float, adaptVal)
psycopg2.extensions.register_adapter( np.float64, adaptVal)
psycopg2.extensions.register_adapter( np.int64, adaptVal)
psycopg2.extensions.register_adapter( np.string_, adaptVal)
augmentDb( bugLev, useCommit, inSpec)
#====================================================================
[docs]def augmentDb( bugLev, useCommit, inSpec):
'''
Adds additional information to the model database table.
See documentation for the :func:`main` function.
**Parameters**:
* bugLev (int): Debug level. Normally 0.
* useCommit (bool): do we commit changes to the DB.
* inSpec (str): Name of JSON file containing DB parameters.
See description at :func:`main`.
**Returns**
* None
'''
print 'augmentDb: useCommit: %s' % (useCommit,)
with open( inSpec) as fin:
specMap = json.load( fin)
dbhost = specMap.get('dbhost', None)
dbport = specMap.get('dbport', None)
dbuser = specMap.get('dbuser', None)
dbpswd = specMap.get('dbpswd', None)
dbname = specMap.get('dbname', None)
dbschema = specMap.get('dbschema', None)
dbtablemodel = specMap.get('dbtablemodel', None)
if dbhost == None: badparms('inSpec name not found: dbhost')
if dbport == None: badparms('inSpec name not found: dbport')
if dbuser == None: badparms('inSpec name not found: dbuser')
if dbpswd == None: badparms('inSpec name not found: dbpswd')
if dbname == None: badparms('inSpec name not found: dbname')
if dbschema == None: badparms('inSpec name not found: dbschema')
if dbtablemodel == None: badparms('inSpec name not found: dbtablemodel')
dbport = int( dbport)
conn = None
cursor = None
try:
conn = psycopg2.connect(
host=dbhost,
port=dbport,
user=dbuser,
password=dbpswd,
database=dbname)
if bugLev >= 1:
print 'main: got conn. dbhost: %s dbport: %d' % (dbhost, dbport,)
cursor = conn.cursor()
cursor.execute('set search_path to %s', (dbschema,))
db_rows = None
sqlMsg = 'SELECT * FROM %s ORDER BY mident' % (dbtablemodel,)
# For testing use:
#WHERE wrapid = '@2013.12.04@09.42.21.719000@ciduser@data.redmesa.raw.As@'
#WHERE wrapid = '@2014.03.25@10.22.24.696353@ciduser@data.peregrine.upload.lany.2013.10.31.raw.fixed.2013_PRB_GGAUvsHSE@'
#WHERE wrapid = '@2014.04.08@19.58.58.011410@ssulliva@scratch.ssulliva.slany.gw@'
#WHERE model.mident = ANY( ARRAY[
# 2739, 2740, 2741, 2755, 2817,
# 2818, 2820, 2835, 2836, 2848,
# 3046, 3047, 3048, 3568, 3880,
# 3881, 3882, 4173, 4592, 4593,
# 4594, 4772, 17126, 17193, 23445,
# 27102, 30971, 31004, 31037, 31070,
# 34729, 39530, 39597, 39634, 39653,
# 39668, 59421, 61882, 64391, 66165])
cursor.execute( sqlMsg)
dbRows = cursor.fetchall()
dbCols = [desc[0] for desc in cursor.description]
nrow = len( dbRows)
ncol = len( dbCols)
objPairs = [] # one pair [curObj, newObj] per row
for row in dbRows:
if len(row) != ncol: throwerr('len mismatch')
# Fill curObj with the values of this row
curObj = DbObj()
for icol in range( ncol):
setattr( curObj, dbCols[icol], row[icol])
objPairs.append( (curObj, DbObj()))
#print '\nxxxxxxxxxxxxxx start objPairs[0][0] xxxxxxxxxxxxxxxxx'
#wrapUpload.printMap('objPairs[0][0]', objPairs[0][0].__dict__, 0)
#print 'xxxxxxxxxxxxxx end objPairs[0][0] xxxxxxxxxxxxxxxxx\n'
addMisc( bugLev, objPairs) # updates newObjs
addBandgap( bugLev, objPairs) # updates newObjs
addChemform( bugLev, objPairs) # updates newObjs
addSpaceGroups( bugLev, objPairs) # updates newObjs
# addMinenergy must be after addSpaceGroups
addMinenergy( bugLev, objPairs) # updates newObjs
addEnthalpy( bugLev, objPairs) # updates newObjs
addMagmom( bugLev, objPairs) # updates newObjs
addDosMat( bugLev, objPairs) # updates newObjs
addDielectric( bugLev, objPairs) # updates newObjs
addFamily( bugLev, objPairs) # updates newObjs
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
keys = newObj.__dict__.keys()
keys.sort()
colNames = []
colVals = []
for key in keys:
colNames.append( key)
colVals.append( newObj.__dict__[key])
if len( colNames) > 0:
colNameStg = ','.join( colNames)
colTagStg = (len(colNames) * '%s,') [:-1] # del final ','
updateStg = 'UPDATE %s SET (%s) = (%s) WHERE mident = %d' \
% ( dbtablemodel, colNameStg, colTagStg, curObj.mident,)
cursor.execute( updateStg, colVals)
# xxx outdent?
if useCommit:
conn.commit()
if bugLev >= 1:
print 'augmentDb: commit done. mident: %d formula: %s num col: %d' \
% (curObj.mident, newObj.formula, len( newObj.__dict__),)
finally:
if cursor != None: cursor.close()
if conn != None: conn.close()
#====================================================================
def addMisc( bugLev, objPairs): # updates newObjs
'''
Calculates the atomnames, volume, pressure, misc, for all structures.
For example for H2O readVasp.py gives
and 2 typeNames ['H', 'O'] with typeNums [2, 1],
and we compute the 3 atomNames ['H', 'H', 'O'].
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
# xxx all doc
nrow = len( objPairs)
if bugLev >= 1:
print 'addMisc entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
#print '\nxxxxxxxxxxxx start curobj xxxxxxxxxxxxxxx'
#wrapUpload.printMap('curObj', curObj.__dict__, 0)
#print '\nxxxxxxxxxxxx end curobj xxxxxxxxxxxxxxx'
#print '\nxxxxxxxxxxxx start newobj xxxxxxxxxxxxxxx'
#wrapUpload.printMap('newObj', newObj.__dict__, 0)
#print '\nxxxxxxxxxxxx end newobj xxxxxxxxxxxxxxx'
if curObj.iterrealtimes == None: newObj.itertotaltime = None
else: newObj.itertotaltime = sum( curObj.iterrealtimes)
if bugLev >= 5:
print 'addMisc: itertotaltime: %s' % (newObj.itertotaltime,)
numType = len( curObj.typenames)
if len( curObj.typenums) != numType: throwerr('numType mismatch')
atomTypes = [] # indices into typeNames, typeNums
atomNames = []
atomMasses_amu = []
atomPseudos = []
atomValences = []
for ii in range( numType):
numEle = curObj.typenums[ii]
atomTypes += numEle * [ ii ]
atomNames += numEle * [ curObj.typenames[ii] ]
atomMasses_amu += numEle * [ curObj.typemasses_amu[ii] ]
atomPseudos += numEle * [ curObj.typepseudos[ii] ]
atomValences += numEle * [ curObj.typevalences[ii] ]
if bugLev >= 5:
print 'addMisc: atomTypes: %s' % (atomTypes,)
print 'addMisc: atomNames: %s' % (atomNames,)
print 'addMisc: atomMasses_amu: %s' % (atomMasses_amu,)
print 'addMisc: atomPseudos: %s' % (atomPseudos,)
print 'addMisc: atomValences: %s' % (atomValences,)
newObj.atomnames = atomNames
newObj.atommasses_amu = atomMasses_amu
newObj.atompseudos = atomPseudos
newObj.atomvalences = atomValences
# totalValence = sum( count[i] * valence[i])
# PyLada calls this valence.
newObj.totalvalence = np.dot( curObj.typenums, curObj.typevalences)
if bugLev >= 5:
print 'addMisc: totalValence: %s' % (newObj.totalvalence,)
# In defect calcs of Lany, we can have numElectron==288
# and totalValence==289, leaving the supercell with
# positive charge, for a donor.
newObj.netcharge = newObj.totalvalence - curObj.numelectron
# volume, Angstrom^3
# The scale is hard coded as 1.0 in PyLada crystal/read.py,
# in both icsd_cif_a and icsd_cif_b.
if getattr( curObj, 'finalbasismat', None) == None:
newObj.finalvolume_ang3 = None
newObj.finaldensity_g_cm3 = None
else:
volScale = 1.0
newObj.finalvolume_ang3 = abs( np.linalg.det(
volScale * np.array( curObj.finalbasismat)))
# Density
# xxx better: get atomic weights from periodic table
volCm3 = newObj.finalvolume_ang3 / (1.e8)**3 # 10**8 Angstrom per cm
totMass = np.dot( curObj.typenums, curObj.typemasses_amu)
totMassGm = totMass * 1.660538921e-24 # 1.660538921e-24 g / amu
newObj.finaldensity_g_cm3 = totMassGm / volCm3
if bugLev >= 5:
print 'finalVolume_ang3: %s' % (newObj.finalvolume_ang3,)
print 'volCm3: %g' % (volCm3,)
print 'totMass: %g' % (totMass,)
print 'totMassGm: %g' % (totMassGm,)
print 'finalDensity_g_cm3: %g' % (newObj.finaldensity_g_cm3,)
# reciprocal space volume, * (2*pi)**3
# As in PyLada.
# Someday need new db col for recipvolume
if getattr( curObj, 'finalbasismat', None) == None:
recipvolume = None
else:
basisMat = volScale * np.array( curObj.finalbasismat)
invMat = np.linalg.inv( basisMat)
recipVolume = abs( np.linalg.det( invMat)) * (2 * np.pi)**3
if bugLev >= 5:
print 'recipVolume: basisMat:\n%s' % (repr(basisMat),)
print 'recipVolume: invMat:\n%s' % (repr(invMat),)
print 'recipVolume: det:\n%s' % (repr(np.linalg.det( invMat)),)
print 'recipVolume: %g' % (recipVolume,)
# Calc pressure
# xxx Caution: we do not include the non-diag terms in:
# VASP: force.F: FORCE_AND_STRESS: line 1410:
# PRESS=(TSIF(1,1)+TSIF(2,2)+TSIF(3,3))/3._q &
# & -DYN%PSTRESS/(EVTOJ*1E22_q)*LATT_CUR%OMEGA
if getattr( curObj, 'finalstressmat_kbar', None) == None:
newObj.finalpressure_kbar = None
else:
newObj.finalpressure_kbar \
= np.array( curObj.finalstressmat_kbar).diagonal().sum() / 3.0
if bugLev >= 5:
print 'getStressMat: finalPressure_kbar: %s' \
% (newObj.finalpressure_kbar,)
#====================================================================
def addBandgap( bugLev, objPairs): # updates newObjs
'''
Uses eigenMat to set cbMinVals, cbMinIxs, cbMin,
and similarly for Max,
and bandgapDirects, bandgapIndirects, and bandgap.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addBandgap entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
if curObj.eigenmat == None:
# Someday need new db cols:
#newObj.cbminvals = None
#newObj.cbminixs = None
#newObj.vbmaxvals = None
#newObj.vbmaxixs = None
#newObj.bandgapdirects = None
#newObj.bandgapindirects = None
newObj.cbmin = None
newObj.vbmax = None
newObj.bandgapdirect = None
newObj.bandgapindirect = None
else:
# Conduction band min
# cbMinMat[isp][ikp] = min of eigvals > efermi0
#
# Valence band max
# vbMaxMat[isp][ikp] = max of eigvals <= efermi0
eigenMat = np.array( curObj.eigenmat)
if bugLev >= 5:
print 'addBandGap: eig shape: ', eigenMat.shape
print 'addBandGap: numSpin: %d numKpoint: %d numBand: %d' \
% (curObj.numspin, curObj.numkpoint, curObj.numband,)
numSpin = curObj.numspin
numKpoint = curObj.numkpoint
# Sometimes numBand == 128 but eigenMat.shape == [x,x,64]
#xxx numBand = curObj.numband
numBand = eigenMat.shape[2]
# For each isp and kpoint, what is min eigvalue > efermi0
cbMinMat = np.empty( [numSpin, numKpoint], dtype=float)
vbMaxMat = np.empty( [numSpin, numKpoint], dtype=float)
cbMinMat.fill( np.inf)
vbMaxMat.fill( -np.inf)
# For each isp, what is min or max across all kpoints
cbMinVals = np.empty( [numSpin], dtype=float)
vbMaxVals = np.empty( [numSpin], dtype=float)
cbMinVals.fill( np.inf)
vbMaxVals.fill( -np.inf)
# For each isp, which kpoint has the min or max
cbMinIxs = np.empty( [numSpin], dtype=int)
vbMaxIxs = np.empty( [numSpin], dtype=int)
cbMinIxs.fill( -1)
vbMaxIxs.fill( -1)
for isp in range( numSpin):
for ikp in range( numKpoint):
for iband in range( numBand):
val = eigenMat[isp][ikp][iband]
algorithm = 'occupMat'
if algorithm == 'fermi0':
if val <= curObj.efermi0:
if val > vbMaxMat[isp][ikp]: vbMaxMat[isp][ikp] = val
if val > vbMaxVals[isp]:
vbMaxVals[isp] = val
vbMaxIxs[isp] = ikp
else: # val > curObj.efermi0
if val < cbMinMat[isp][ikp]: cbMinMat[isp][ikp] = val
if val < cbMinVals[isp]:
cbMinVals[isp] = val
cbMinIxs[isp] = ikp
elif algorithm == 'occupMat':
if numSpin == 1:
occLo = 0.02
occHi = 1.98
elif numSpin == 2:
occLo = 0.01
occHi = 0.99
else: throwerr('unknown numSpin')
occ = curObj.occupmat[isp][ikp][iband]
if occ >= occHi:
if val > vbMaxMat[isp][ikp]: vbMaxMat[isp][ikp] = val
if val > vbMaxVals[isp]:
vbMaxVals[isp] = val
vbMaxIxs[isp] = ikp
if occ <= occLo:
if val < cbMinMat[isp][ikp]: cbMinMat[isp][ikp] = val
if val < cbMinVals[isp]:
cbMinVals[isp] = val
cbMinIxs[isp] = ikp
else: throwerr('unknown algorithm')
# Find bandgapDirects[isp] = min gap for the same kpoint
# Find bandgapIndirects[isp] = min gap across kpoints
bandgapDirects = numSpin * [np.inf]
bandgapIndirects = numSpin * [np.inf]
for isp in range( numSpin):
for ikp in range( numKpoint):
gap = max( 0, cbMinMat[isp][ikp] - vbMaxMat[isp][ikp])
if gap < bandgapDirects[isp]: bandgapDirects[isp] = gap
bandgapIndirects[isp] = max( 0, cbMinVals[isp] - vbMaxVals[isp])
# Someday need new db cols:
#newObj.cbminvals = cbMinVals
#newObj.cbminixs = cbMinIxs
#newObj.vbmaxvals = vbMaxVals
#newObj.vbmaxixs = vbMaxIxs
#newObj.bandgapdirects = bandgapDirects
#newObj.bandgapindirects = bandgapIndirects
newObj.cbmin = min( cbMinVals)
newObj.vbmax = max( vbMaxVals)
newObj.bandgapdirect = min( bandgapDirects)
newObj.bandgapindirect = min( bandgapIndirects)
# This is the PyLada version
#newObj.cbmin = min( cbMinVals)
#newObj.vbmax = max( vbMaxVals)
#newObj.bandgap = newObj.cbmin - newObj.vbmax
#if newObj.bandgap < 0: newObj.bandgap = 0
if bugLev >= 5:
print '\naddBandGap: efermi0: %s' % (curObj.efermi0,)
for isp in range( numSpin):
print '\naddBandGap: start isp: %d' % (isp,)
for ikp in range( numKpoint):
vval = vbMaxMat[isp][ikp]
cval = cbMinMat[isp][ikp]
vmsg = '%8.4f' % (vval,)
cmsg = '%8.4f' % (cval,)
if ikp == vbMaxIxs[isp]: vmsg += ' ***'
if ikp == cbMinIxs[isp]: cmsg += ' ***'
print ' isp: %d ikp: %3d vbMaxMat: %-14s cbMinMat: %-14s' \
% ( isp, ikp, vmsg, cmsg,)
print ''
print 'addBandGap: vbMaxIxs: %s' % (vbMaxIxs,)
print 'addBandGap: vbMaxVals: %s' % (vbMaxVals,)
print 'addBandGap: cbMinIxs: %s' % (cbMinIxs,)
print 'addBandGap: cbMinVals: %s' % (cbMinVals,)
print 'addBandGap: bandgapDirects: %s' % (bandgapDirects,)
print 'addBandGap: bandgapIndirects: %s' % (bandgapIndirects,)
print 'addBandGap: newObj.cbmin: %s' % (newObj.cbmin,)
print 'addBandGap: newObj.vbmax: %s' % (newObj.vbmax,)
print 'addBandGap: newObj.bandgapdirect: %s' \
% (newObj.bandgapdirect,)
print 'addBandGap: newObj.bandgapindirect: %s' \
% (newObj.bandgapindirect,)
# Write plot data to tempplota%ispin
for isp in range( numSpin):
with open('tempplota%d' % (isp,), 'w') as fout:
for ikp in range( numKpoint):
print >> fout, '%g %g' % (vbMaxMat[isp,ikp], cbMinMat[isp,ikp],)
# Write plot data to tempplotb%ispin
for isp in range( numSpin):
with open('tempplotb%d' % (isp,), 'w') as fout:
for ikp in range( numKpoint):
for iband in range( numBand):
print >> fout, '%d %g' \
% (ikp, eigenMat[isp][ikp][iband],)
#====================================================================
def addChemform( bugLev, objPairs): # updates newObjs
'''
Calculates the formula, sortedformula, and chemtext fields.
Set attributes:
formula:
chemical formula, determined from typenames and typenums,
and factoring by the greatest common divison.
For example, ``'H2 O'``.
sortedformula:
like formula, but the elements are sorted alphabetically.
For example, ``'H2 O'``.
chemtext:
Same as the formula, but every token is surrounded
spaces and singletons have an explicit '1',
for ease in parsing.
For example, ``' H 2 O 1 '``.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addChemform entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
tnames = curObj.typenames
tnums = curObj.typenums
if tnames != None and tnums != None:
tlen = len( tnames)
if len( tnums) != tlen: throwerr('tlen mismatch')
# Replace tnums with tnums / gcd( tnums)
if tlen == 1: tnums = [1]
else:
gcd = tnums[0]
for jj in range(1,tlen):
gcd = fractions.gcd( gcd, tnums[jj])
for jj in range(tlen):
tnums[jj] /= gcd
# For chemtext we insert an initial and final blank,
# so every value is surrounded by spaces,
# for easier handling in SQL.
formula = ''
chemtext = ''
for ii in range(tlen):
# Coordinate formula format with pyramid: views.py: vwQueryStd
if len(formula) > 0: formula += ' '
formula += tnames[ii]
if tnums[ii] != 1:
formula += str(tnums[ii])
chemtext += ' %s %d' % (tnames[ii], tnums[ii],)
chemtext += ' '
# Index sort for parallel arrays tnames, tnums,
# so we can calc sortedformula.
ixs = range(tlen)
ixs.sort( lambda ia, ib: cmp( tnames[ia], tnames[ib]))
sortedformula = ''
for ii in range(tlen):
# Coordinate formula format with pyramid: views.py: vwQueryStd
if len(sortedformula) > 0: sortedformula += ' '
sortedformula += tnames[ixs[ii]]
if tnums[ixs[ii]] != 1:
sortedformula += str(tnums[ixs[ii]])
if bugLev >= 5:
print ('addChemform: mident: %d tnames: %s tnums: %s'
+ ' formula: %s sortedformula: %s chemtext: %s') \
% (mident, tnames, tnums, repr( formula),
repr( sortedformula), repr( chemtext),)
newObj.formula = formula
newObj.sortedformula = sortedformula
newObj.chemtext = chemtext
else:
newObj.formula = None
newObj.sortedformula = None
newObj.chemtext = None
#====================================================================
def addSpaceGroups( bugLev, objPairs): # updates newObjs
'''
Calculates the initial and final space groups, for all structures.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addSpaceGroups entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
atomnames = curObj.atomnames
initialbasismat = curObj.initialbasismat
initialcartposmat = curObj.initialcartposmat
finalbasismat = curObj.finalbasismat
finalcartposmat = curObj.finalcartposmat
if initialbasismat == None \
or initialcartposmat == None \
or finalbasismat == None \
or finalcartposmat == None \
or atomnames == None:
initialSgName = None
initialSgNum = None
finalSgName = None
finalSgNum = None
else:
(initialSgName, initialSgNum) = getSpaceGroup(
initialbasismat, initialcartposmat, atomnames)
(finalSgName, finalSgNum) = getSpaceGroup(
finalbasismat, finalcartposmat, atomnames)
newObj.initialspacegroupname = initialSgName
newObj.initialspacegroupnum = initialSgNum
newObj.finalspacegroupname = finalSgName
newObj.finalspacegroupnum = finalSgNum
#====================================================================
def getSpaceGroup( rowCellMat, cartPosMat, atomnames):
# Create the ASE bulk look-alike object needed by spglib
aseBulk = AseBulk( rowCellMat, cartPosMat, atomnames)
# Call spglib
scaleFactor = 1.0
symprec = 1.e-2
sgStg = spglib.get_spacegroup( aseBulk, symprec)
# Parse the resulting string into name, number.
mat = re.match(r'^(.+)\((\d+)\)$', sgStg)
if mat == None:
throwerr('unknown sgStg. mident: %d sgStg: %s' % (mident, sgStg,))
sgName = mat.group(1).strip()
sgNum = int( mat.group(2).strip())
# spglib returns bogus name when num is 0.
if sgNum == 0: sgName = None
return (sgName, sgNum)
#====================================================================
# addMinenergy must be after addSpaceGroups
[docs]def addMinenergy( bugLev, objPairs): # updates newObjs
'''
Calculates the min energyperatom for each formula.
minenergyid == the mident of the row having the min energyperatom
for this formula.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addMinenergy entry: nrow: %d' % (nrow,)
# Find the min energy for each group
# having standard=fere and an ICSD SG,
# and the same formula, same ICSD-SG, same PP, same author
minMap = {} # (formula,finalspacegroupnum,etc) -> [minEnergy, mident]
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
standards = curObj.meta_standards
sgnum = newObj.finalspacegroupnum
formula = curObj.formula
pseudos = curObj.typepseudos
lastName = curObj.meta_lastname
firstName = curObj.meta_firstname
if 'fere' in standards \
and None not in [sgnum, formula, pseudos, lastName, firstName]:
keyStg = '%s,%s,%s,%s,%s' \
% (formula, sgnum, pseudos, lastName, firstName,)
mident = curObj.mident
energy = curObj.energyperatom
pair = minMap.get( keyStg, None)
if pair == None or energy < pair[0]:
minMap[keyStg] = (energy, mident)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
standards = curObj.meta_standards
sgnum = newObj.finalspacegroupnum
formula = curObj.formula
pseudos = curObj.typepseudos
lastName = curObj.meta_lastname
firstName = curObj.meta_firstname
# mident having the min energy for this formula
minId = None
if 'fere' in standards \
and None not in [sgnum, formula, pseudos, lastName, firstName]:
keyStg = '%s,%s,%s,%s,%s' \
% (formula, sgnum, pseudos, lastName, firstName,)
pair = minMap[keyStg] # [minEnergy, minId]
minId = pair[1]
if bugLev >= 5:
print 'addMinenergy: mident: %d formula: %s sgnum: %s minId: %d' \
% (mident, formula, sgnum, minId,)
newObj.minenergyid = minId
#====================================================================
[docs]def addEnthalpy( bugLev, objPairs): # updates newObjs
'''
Calculates the enthalpy of formation per atom, for all structures.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addEnthalpy entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
formula = curObj.formula
typenames = curObj.typenames
typenums = curObj.typenums
typepseudos = curObj.typepseudos
netcharge = curObj.netcharge
meta_standards = curObj.meta_standards
energy = curObj.energyperatom
enthalpy = None
if 'fere' in meta_standards and netcharge == 0:
enthalpy = calcEnthalpy(
bugLev, mident, typenames, typenums, typepseudos, energy)
if bugLev >= 5:
print 'addEnthalpy: mident: %d formula: %s energy: %s enthalpy: %s' \
% (mident, formula, energy, enthalpy,)
newObj.enthalpy = enthalpy
#====================================================================
def addMagmom( bugLev, objPairs): # updates newObjs
'''
Calculates the magnetic moment for all structures.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addMagmom entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
mident = curObj.mident
efermi0 = curObj.efermi0
eigenmat = np.array( curObj.eigenmat)
magmom = None
if efermi0 != None and eigenmat != None and eigenmat.shape[0] == 2:
eig0 = eigenmat[0].flatten()
eig1 = eigenmat[1].flatten()
occ0 = len([x for x in eig0 if x <= efermi0])
occ1 = len([x for x in eig1 if x <= efermi0])
nkpts = eigenmat.shape[1]
magmom = abs( occ0 - occ1) / float( nkpts)
if bugLev >= 5:
print ('addMagmom: mident: %d occ0: %d occ1: %d' \
+ ' nkpts: %d magmom: %g') \
% ( mident, occ0, occ1, nkpts, magmom,)
newObj.magneticmoment = magmom
#====================================================================
def addDosMat( bugLev, objPairs): # updates newObjs
'''
Removes the dosmat value from entries where it doesn't belong.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addDosMat entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
dosmat = curObj.dosmat
if 'post-dos' in curObj.meta_standards and dosmat != None:
dosmat = np.array( dosmat)
if curObj.numspin == 1:
dosmat = dosmat[ :, :2] # just keep E, DOS
elif curObj.numspin == 2:
dosmat = dosmat[ :, :3] # just keep E, DOSup, DOSdown
# xxx start here. We want:
# dosmat[ :, 1] = dosmat[:,1] - newObj.vbmax
# but each time we run aug, we decrement dosmat.
#
# xxx separate all destructive changes: dosmat, dielectric*
else: dosmat = None
newObj.dosmat = dosmat
#====================================================================
def addDielectric( bugLev, objPairs): # updates newObjs
'''
Removes the dosmat value from entries where it doesn't belong.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addDielectric entry: nrow: %d' % (nrow,)
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj]
if 'post-lopt' not in curObj.meta_standards:
newObj.dielectricimag = None
newObj.dielectricreal = None
#====================================================================
def addFamily( bugLev, objPairs): # updates newObjs
'''
Removes the dosmat value from entries where it doesn't belong.
**Parameters**:
* bugLev (int): Debug level; normally 0.
* curObjs: current objects, one per DB row
* newObjs: objects with new fields only, one per DB row
'''
nrow = len( objPairs)
if bugLev >= 1:
print 'addFamily entry: nrow: %d' % (nrow,)
shaMap = {} # curObj.hashstring -> pair [curObj, newObj]
for iobj in range(nrow):
pair = objPairs[iobj]
shaMap[pair[0].hashstring] = pair
# Set parentids
for iobj in range(nrow):
(curObj, newObj) = objPairs[iobj] # pair [curObj, newObj]
parentids = None
parStgs = curObj.meta_parents
if parStgs != None:
parentids = [shaMap[parStg][0].mident for parStg in parStgs]
newObj.parentids = parentids
# Set familyid = most distant ancestor
for iobj in range(nrow):
(curObjOrig, newObjOrig) = objPairs[iobj] # pair [curObj, newObj]
# Run up the parent chain to the most distant ancestor
(curObjAnc, newObjAnc) = (curObjOrig, newObjOrig)
while curObjAnc.meta_parents != None:
# Only use the first parent
(curObjAnc, newObjAnc) = shaMap[ curObjAnc.meta_parents[0]]
if curObjAnc != curObjOrig:
newObjAnc.familyid = curObjAnc.mident # set ancestor's familyid
newObjOrig.familyid = curObjAnc.mident # set our familyid
#====================================================================
[docs]def calcEnthalpy(
bugLev,
mident,
typenames,
typenums,
typepseudos,
energy):
'''
Calculates the enthalpy of formation per atom, for a single structure.
**Parameters (using example magnetitie, Fe3 O4)**:
* mident (int): The unique DB identifier.
* typenames (str[]): List of unique type names, ex. ['Fe', 'O']
* typenums (int[]): Number of each unique type, ex. [3, 4]
* typepseudos (str[]): Pseudopotential for each type,
['PAW_PBE Fe 06Sep2000', 'PAW_PBE O_s 07Sep2000']
* energy (float): final total energy without entropy, per atom in cell
**Returns**
* enthalpy (float): enthalpy of formation per atom
See:
DOI: 10.1103/PhysRevB.85.115104
http://prb.aps.org/pdf/PRB/v85/i11/e115104
Correcting density functional theory for accurate predictions
of compound enthalpies of formation:
Fitted elemental-phase reference energies
Vladan Stevanovic, Stephan Lany, Xiuwen Zhang, Alex Zunger
page 11, Table V: mu^{FERE}
Email from Stevanovic, Monday, August 12, 2013 5:45 PM
For the enthalpy of formation calculations I will use the
following example. So, if we have a ternary compound with the
chemical formula A_mB_nX_l, where A, B and X represent the
symbols of the elements forming the compound and if N=m+n+l the
total number of atoms in the compound the enthalpy of formation
can be calculated as:
dHf (eV/atom) = [Etot(eV/atom) * N - (m*musMap['A']
+ n*musMap['B'] + l*musMap['X'] ) ] / N
where Etot(eV/atom) is the VASP total energy per atom which is
currently listed in the database and musMap['A'],
musMap['B'], and musMap['X'] you read from the attached
dictionary.
This is the same as:
enthalpy = - sum_i( n_i * mus_i) / sum( n_i)
where n_i = number of element i, and mus_i = mus value for element i.
We use this calculation for ALL structures, not just ternaries.
'''
musList = [
# Element Pseudopotential FERE value
# ------- --------------- ----------
[ 'Ag', 'Ag', -0.82700958541595615, ],
[ 'Al', 'Al', -3.02, ],
[ 'As', 'As', -5.06, ],
[ 'Au', 'Au', -2.2303410086960551, ],
[ 'Ba', 'Ba_sv', -1.3944992462870172, ],
[ 'Be', 'Be', -3.3972092621754264, ],
[ 'Bi', 'Bi_d', -4.3853003286558812, ],
[ 'Ca', 'Ca_pv', -1.64, ],
[ 'Cd', 'Cd', -0.56, ],
# but if LDAUTYPE=2 and U=5eV, use Cd: +0.19 # xxx
[ 'Cl', 'Cl', -1.6262437135301639, ],
[ 'Co', 'Co', -4.7543486260270402, ],
[ 'Cr', 'Cr_pv', -7.2224146752384204, ],
[ 'Cu', 'Cu', -1.9725806522979044, ],
[ 'F', 'F', -1.7037867766570287, ],
[ 'F', 'F_s', -1.68, ],
[ 'Fe', 'Fe', -6.1521343161090325, ],
[ 'Ga', 'Ga_d', -2.37, ],
[ 'Ge', 'Ge_d', -4.137439286830797, ],
[ 'Hf', 'Hf_pv', -7.397695761161847, ],
[ 'Hg', 'Hg', -0.12361566177444684, ],
[ 'In', 'In_d', -2.31, ],
[ 'Ir', 'Ir', -5.964577394407752, ],
[ 'K', 'K_sv', -0.80499202755075006, ],
[ 'La', 'La', -3.6642174822805287, ],
[ 'Li', 'Li_sv', -1.6529591953887741, ],
[ 'Mg', 'Mg', -0.99, ],
[ 'Mg', 'Mg_pv', -0.99, ],
[ 'Mn', 'Mn', -6.9965778258511993, ],
[ 'N', 'N', -8.51, ],
[ 'N', 'N_s', -8.51, ],
[ 'Na', 'Na_pv', -1.0640326227725869, ],
[ 'Nb', 'Nb_pv', -6.6867516375690608, ],
[ 'Ni', 'Ni', -3.5687859474688026, ],
[ 'O', 'O_s', -4.76, ],
[ 'O', 'O', -4.73, ],
[ 'P', 'P', -5.64, ],
[ 'Pd', 'Pd', -3.1174044624888873, ],
[ 'Pt', 'Pt', -3.9527597082085424, ],
[ 'Rb', 'Rb_sv', -0.6750560483522855, ],
[ 'Rh', 'Rh', -4.7622899695820369, ],
[ 'S', 'S', -4.00, ],
[ 'Sb', 'Sb', -4.2862260747305099, ],
[ 'Sc', 'Sc_sv', -4.6302422200922519, ],
[ 'Se', 'Se', -3.55, ],
[ 'Si', 'Si', -4.9927748122726356, ],
[ 'Sn', 'Sn_d', -3.7894939351245469, ],
[ 'Sr', 'Sr_sv', -1.1674559193419329, ],
[ 'Ta', 'Ta_pv', -8.8184831379805324, ],
[ 'Te', 'Te', -3.2503408197224912, ],
[ 'Ti', 'Ti_pv', -5.5167842601434147, ],
[ 'V', 'V_pv', -6.4219725884764864, ],
[ 'Y', 'Y_sv', -4.812621315561298, ],
[ 'Zn', 'Zn', -0.84, ],
# but if LDAUTYPE=2 and U=6eV, use Zn: -0.51 # xxx
[ 'Zr', 'Zr_sv', -5.8747056261113126, ],
]
musMap = {}
for ele in musList:
musMap[ele[1]] = ele[2]
if typenums == None or len(typenums) == 0: enthalpy = None
else:
totNum = sum( typenums)
enthalpy = energy
for ii in range(len(typepseudos)):
pseudoStg = typepseudos[ii] # like 'PAW_PBE Ca_pv 06Sep2000'
num = typenums[ii]
toks = pseudoStg.strip().split()
if len(toks) != 3:
throwerr('bad pseudo. mident: %s names: %s nums: %s pseudos: %s'
% (mident, typenames, typenums, typepseudos,))
pseudoName = toks[1]
mufere = musMap.get( pseudoName, None)
if bugLev >= 5:
print 'calcEnthalpy: mident: %d pseudoStg: %s mufere: %s' \
% (mident, pseudoStg, mufere,)
if mufere == None:
print 'calcEnthalpy: warning: no pseudo "%s". mident: %s names: %s nums: %s pseudos: %s' \
% (pseudoName, mident, typenames, typenums, typepseudos,)
enthalpy = None
if enthalpy != None:
enthalpy -= num * mufere / float(totNum)
if bugLev >= 5:
print 'calcEnthalpy: mident: %d enthalpy: %s' % (mident, enthalpy,)
return enthalpy
#====================================================================
[docs]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()