Commit ae6a31fb authored by Pardini, Lorenzo (lopa)'s avatar Pardini, Lorenzo (lopa)
Browse files

fixed dos units and gw

parent cc0eb5aa
......@@ -6,7 +6,7 @@ from nomadcore.unit_conversion.unit_conversion import convert_unit
from nomadcore.unit_conversion import unit_conversion
class DosHandler(xml.sax.handler.ContentHandler):
def __init__(self, backend, spinTreat):
def __init__(self, backend, spinTreat, unitCellVol):
self.backend = backend
self.dosSectionGIndex = -1
self.inDos = False
......@@ -24,6 +24,7 @@ class DosHandler(xml.sax.handler.ContentHandler):
self.numDosVal = 0
self.dosProjDummy = []
self.dosProjDummy2 = []
self.unitCellVol = float(unitCellVol[0])
def endDocument(self):
# self.inDos = False
......@@ -32,10 +33,13 @@ class DosHandler(xml.sax.handler.ContentHandler):
# self.backend.closeSection("section_atom_projected_dos",self.dosProjSectionGIndex)
self.dosSectionGIndex = -1
# self.dosProjSectionGIndex = -1
# print("self.unitCellVol=",self.unitCellVol)
def startElement(self, name, attrs):
ha_per_joule = convert_unit(1, "hartree", "J")
bohr_cube_to_m_cube = convert_unit(1, "bohr^3", "m^3")
# bohr_cube_to_m_cube = convert_unit(1, "bohr^3", "m^3")
# print("bohr_cube_to_m_cube = ", bohr_cube_to_m_cube )
# print("ha_per_joule= ",ha_per_joule)
fromH = unit_conversion.convert_unit_function("hartree", "J")
if name == "totaldos":
self.dosSectionGIndex = self.backend.openSection("section_dos")
......@@ -50,11 +54,11 @@ class DosHandler(xml.sax.handler.ContentHandler):
self.inDosProj = True
elif name == "point":
if self.inDos:
self.totDos.append(ha_per_joule*bohr_cube_to_m_cube*float(attrs.getValue('dos')))
self.totDos.append(ha_per_joule*self.unitCellVol*float(attrs.getValue('dos')))
# self.energy.append(float(attrs.getValue('e')))
self.energy.append(fromH(float(attrs.getValue('e'))))
elif self.inDosProj:
self.dosProj.append(ha_per_joule*bohr_cube_to_m_cube*float(attrs.getValue('dos')))
self.dosProj.append(ha_per_joule*self.unitCellVol*float(attrs.getValue('dos')))
# self.energy.append(float(attrs.getValue('e')))
self.energy.append(fromH(float(attrs.getValue('e'))))
elif name == "diagram":
......@@ -154,8 +158,8 @@ class DosHandler(xml.sax.handler.ContentHandler):
def characters(self, content):
pass
def parseDos(inF, backend, spinTreat):
handler = DosHandler(backend, spinTreat)
def parseDos(inF, backend, spinTreat, unitCellVol):
handler = DosHandler(backend, spinTreat, unitCellVol)
logging.error("will parse")
xml.sax.parse(inF, handler)
logging.error("did parse")
......@@ -21,8 +21,10 @@ class GWParser(object):
self.vertexLabels = []
self.vertexNum = 0
def parseGW(self, gwFile, backend, dftMethodSectionGindex, dftSingleConfigurationGindex):
def parseGW(self, gwFile, backend, dftMethodSectionGindex, dftSingleConfigurationGindex, xcName, unitCellVol):
# logging.error("GW onClose_section_single_configuration_calculation")
# print("xcNameGW=", xcName)
self.unitCellVol = float(unitCellVol[0])
backend.openNonOverlappingSection("section_single_configuration_calculation")
if dftSingleConfigurationGindex is not None:
backend.openNonOverlappingSection("section_calculation_to_calculation_refs")
......@@ -48,6 +50,7 @@ class GWParser(object):
if os.path.exists(inputFile):
selfGWSetGIndex = backend.openSection("section_method")
backend.addValue('electronic_structure_method', "G0W0")
backend.addValue('x_exciting_gw_starting_point', xcName)
if dftMethodSectionGindex is not None:
m2mGindex = backend.openNonOverlappingSection("section_method_to_method_refs")
backend.addValue("method_to_method_ref", dftMethodSectionGindex)
......@@ -58,20 +61,25 @@ class GWParser(object):
npol = 0
scrtype = "rpa"
snempty = 0
pnempty = 0
coreflag = "all"
fgrid = "gaule2"
lmaxmb = 3
epsmb = 0.0001
# gmb = 1.0
sciavtype = "isotropic"
k1 = 0
k2 = 0
# f1 = 0
# f2 = 0
s1 = 0
s2 = 0
# m1 = 0
# m2 = 0
# bc1 = 0
# bc2 = 0
# sc1 = 0
# sc2 = 0
m1 = 0
m2 = 0
bc1 = 0
bc2 = 0
sc1 = 0
sc2 = 0
with open(inputFile) as g:
i = 0
while 1:
......@@ -86,12 +94,12 @@ class GWParser(object):
# if s[0] == "</freqgrid>": f2 = i
if s[0] == "<selfenergy": s1 = i
if s[0] == "</selfenergy>": s2 = i
# if s[0] == "<mixbasis": m1 = i
# if s[0] == "</mixbasis>": m2 = i
# if s[0] == "<barecoul": bc1 = i
# if s[0] == "</barecoul>": bc2 = i
# if s[0] == "<scrcoul": sc1 = i
# if s[0] == "</scrcoul>": sc2 = i
if s[0] == "<mixbasis": m1 = i
if s[0] == "</mixbasis>": m2 = i
if s[0] == "<barecoul": bc1 = i
if s[0] == "</barecoul>": bc2 = i
if s[0] == "<scrcoul": sc1 = i
if s[0] == "</scrcoul>": sc2 = i
with open(inputFile) as g:
i = 0
while 1:
......@@ -108,19 +116,35 @@ class GWParser(object):
actype = s[1][1:-1]
if (s[0] == "npol") and (i >= k1):
npol = s[1][1:-1]
if (s[0] == "nempty") and (i >= k1) and (i <= k2):
pnempty = s[1][1:-1]
if (s[0] == "scrtype") and (i >= k1):
scrtype = s[1][1:-1]
if (s[0] == "nempty") and (i >= s1) and (i <= s2):
snempty = s[1][1:-1]
if (s[0] == "nempty") and (i >= m1) and (i <= m2):
lmaxmb = s[1][1:-1]
if (s[0] == "nempty") and (i >= m1) and (i <= m2):
epsmb = s[1][1:-1]
# if (s[0] == "nempty") and (i >= m1) and (i <= m2):
# gmb = s[1][1:-1]
if (s[0] == "sciavtype") and (i >= sc11) and (i <= sc2):
sciavtype = s[1][1:-1]
if (s[0] == "fgrid") and (i >= k1):
fgrid = s[1][1:-1]
backend.addValue("x_exciting_GW_frequency_grid_type", fgrid)
backend.addValue("x_exciting_GW_self_energy_c_empty_states", int(snempty))
backend.addValue("x_exciting_GW_core_treatment", coreflag)
backend.addValue("x_exciting_GW_self_energy_singularity_treatment", singularity)
backend.addValue("x_exciting_GW_self_energy_c_analytical_continuation", actype)
backend.addValue("x_exciting_GW_self_energy_c_number_of_poles", int(npol))
backend.addValue("x_exciting_GW_screened_Coulomb", scrtype)
backend.addValue("x_exciting_gw_frequency_grid_type", fgrid)
backend.addValue("x_exciting_gw_self_energy_c_empty_states", int(snempty))
backend.addValue("x_exciting_gw_core_treatment", coreflag)
backend.addValue("x_exciting_gw_self_energy_singularity_treatment", singularity)
backend.addValue("x_exciting_gw_self_energy_c_analytical_continuation", actype)
backend.addValue("x_exciting_gw_self_energy_c_number_of_poles", int(npol))
backend.addValue("x_exciting_gw_screened_Coulomb", scrtype)
backend.addValue("x_exciting_gw_polarizability_empty_states", int(pnempty))
backend.addValue("x_exciting_gw_basis_set", "mixed")
backend.addValue("x_exciting_gw_mixed_basis_lmax", lmaxmb)
backend.addValue("x_exciting_gw_mixed_basis_tolerance", epsmb)
# backend.addValue("x_exciting_gw_mixed_basis_gmax", gmb)
backend.addValue("x_exciting_gw_screened_coulomb_volume_average",sciavtype)
backend.closeSection("section_method",selfGWSetGIndex)
if os.path.exists(vertexGWFile):
......@@ -192,17 +216,17 @@ class GWParser(object):
backend.addValue("number_of_eigenvalues", len(qpE[0]))
backend.addValue("number_of_eigenvalues_kpoints", len(qpGWKpoint))
backend.addValue("eigenvalues_values", qpE)
backend.addValue("x_exciting_GW_qp_linearization_prefactor", Znk)
backend.addValue("x_exciting_gw_qp_linearization_prefactor", Znk)
backend.closeSection("section_eigenvalues",eigvalGWGIndex)
backend.addValue("x_exciting_GW_self_energy_x", Sx)
backend.addValue("x_exciting_GW_self_energy_c", Sc)
backend.addValue("x_exciting_gw_self_energy_x", Sx)
backend.addValue("x_exciting_gw_self_energy_c", Sc)
####################DOS######################
if os.path.exists(dosGWFile):
dosGWGIndex = backend.openSection("section_dos")
ha_per_joule = unit_conversion.convert_unit(1, "hartree", "J")
bohr_cube_to_m_cube = unit_conversion.convert_unit(1, "bohr^3", "m^3")
# bohr_cube_to_m_cube = unit_conversion.convert_unit(1, "bohr^3", "m^3")
fromH = unit_conversion.convert_unit_function("hartree", "J")
with open(dosGWFile) as g:
dosValues = [[],[]]
......@@ -212,7 +236,7 @@ class GWParser(object):
if not s: break
s = s.strip()
s = s.split()
ene, value = fromH(float(s[0])), ha_per_joule*bohr_cube_to_m_cube*float(s[1])
ene, value = fromH(float(s[0])), ha_per_joule*self.unitCellVol*float(s[1])
dosEnergies.append(ene)
if not self.spinTreat:
for i in range(0,2):
......
......@@ -30,6 +30,8 @@ class ExcitingParserContext(object):
self.sim_cell = []
self.cell_format = ''
self.secRunIndex = None
self.unit_cell_vol = 0
self.xcName = None
def onOpen_section_system(self, backend, gIndex, section):
self.secSystemIndex = gIndex
......@@ -50,14 +52,16 @@ class ExcitingParserContext(object):
dirPath = os.path.dirname(self.parser.fIn.name)
gw_File = os.path.join(dirPath, "GW_INFO.OUT")
gwFile = os.path.join(dirPath, "GWINFO.OUT")
# print("xcName= ",self.xcName)
for gFile in [gw_File, gwFile]:
if os.path.exists(gFile):
# logging.error("Starting GW")
gwParser = exciting_parser_gw.GWParser()
gwParser.parseGW(gFile, backend,
dftMethodSectionGindex = self.secMethodIndex,
dftSingleConfigurationGindex = self.secSingleConfIndex)
dftSingleConfigurationGindex = self.secSingleConfIndex,
xcName = self.xcName,
unitCellVol = self.unit_cell_vol)
# logging.error("Finished GW")
break
......@@ -74,6 +78,8 @@ class ExcitingParserContext(object):
backend.addValue("simulation_cell", cell)
def onClose_x_exciting_section_reciprocal_lattice_vectors(self, backend, gIndex, section):
# self.unit_cell_vol = section["x_exciting_unit_cell_volume"]
# print("self.unit_cell_vol= ",self.unit_cell_vol)
recLatticeX = section["x_exciting_geometry_reciprocal_lattice_vector_x"]
recLatticeY = section["x_exciting_geometry_reciprocal_lattice_vector_y"]
recLatticeZ = section["x_exciting_geometry_reciprocal_lattice_vector_z"]
......@@ -84,6 +90,7 @@ class ExcitingParserContext(object):
def onClose_x_exciting_section_xc(self, backend, gIndex, section):
xcNr = section["x_exciting_xc_functional"][0]
# print("xcNr= ",xcNr)
xc_internal_map = {
2: ['LDA_C_PZ', 'LDA_X_PZ'],
3: ['LDA_C_PW', 'LDA_X_PZ'],
......@@ -98,6 +105,8 @@ class ExcitingParserContext(object):
406: ['HYB_GGA_XC_PBEH']
}
for xcName in xc_internal_map[xcNr]:
self.xcName = xcName
# print("xcName= ",self.xcName)
gi = backend.openSection("section_XC_functionals")
backend.addValue("XC_functional_name", xcName)
backend.closeSection("section_XC_functionals", gi)
......@@ -115,7 +124,7 @@ class ExcitingParserContext(object):
if os.path.exists(dosFile):
with open(dosFile) as f:
exciting_parser_dos.parseDos(f, backend, self.spinTreat)
exciting_parser_dos.parseDos(f, backend, self.spinTreat, self.unit_cell_vol)
if os.path.exists(bandFile):
with open(bandFile) as g:
exciting_parser_bandstructure.parseBand(g, backend, self.spinTreat)
......@@ -236,6 +245,8 @@ class ExcitingParserContext(object):
def onClose_section_system(self, backend, gIndex, section):
self.unit_cell_vol = section["x_exciting_unit_cell_volume"]
# print("self.unit_cell_vol= ",self.unit_cell_vol)
backend.addArrayValues('configuration_periodic_dimensions', np.asarray([True, True, True]))
self.secSystemDescriptionIndex = gIndex
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment