exciting_parser_gw.py 14.6 KB
Newer Older
1
2
from builtins import object
import setup_paths
3
import xml.sax
4
5
6
7
from nomadcore.simple_parser import mainFunction, CachingLevel
from nomadcore.simple_parser import SimpleMatcher as SM
from nomadcore.local_meta_info import loadJsonFile, InfoKindEl
from nomadcore.unit_conversion import unit_conversion
8
import os, sys, json, logging, exciting_parser_input
9
10

################################################################
11
# This is the subparser for the exciting GW output
12
13
14
################################################################


Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
15
class GWParser(object):
16
17
18
19
    """context for wien2k In2 parser"""

    def __init__(self):
        self.spinTreat = None
20
21
        self.vertexDist = []
        self.vertexLabels = []
22
        self.vertexNum = 0
23

Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
24
    def parseGW(self, gwFile, backend,  dftMethodSectionGindex, dftSingleConfigurationGindex, xcName, unitCellVol,gmaxvr):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
25
#        logging.error("GW onClose_section_single_configuration_calculation")
26
#        print("xcNameGW=", xcName)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
27
28
        self.gmaxvr = float(gmaxvr[0])
#        print("gkmax= ",self.gmaxvr)
29
        self.unitCellVol = float(unitCellVol[0])
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
30
31
32
33
34
35
        backend.openNonOverlappingSection("section_single_configuration_calculation")
        if dftSingleConfigurationGindex is not None:
            backend.openNonOverlappingSection("section_calculation_to_calculation_refs")
            backend.addValue("calculation_to_calculation_ref", dftSingleConfigurationGindex)
            backend.addValue("calculation_to_calculation_kind", "starting_point")
            backend.closeNonOverlappingSection("section_calculation_to_calculation_refs")
36

Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
37
        dirPath = os.path.dirname(gwFile)
38
39
40
41
42
43
        if os.access(os.path.join(dirPath, "EVALQP.DAT"), os.F_OK):
            eigvalGWFile = os.path.join(dirPath, "EVALQP.DAT")
        elif os.access(os.path.join(dirPath, "EVALQP.TXT"), os.F_OK):
            eigvalGWFile = os.path.join(dirPath, "EVALQP.TXT")
        else:
            pass
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
44
        dosGWFile = os.path.join(dirPath, "TDOS-QP.OUT")
45
46
        bandCarbGWFile = os.path.join(dirPath, "bandstructure-qp.dat")
        bandBorGWFile = os.path.join(dirPath, "BAND-QP.OUT")
47
48
        vertexGWFile = os.path.join(dirPath, "BANDLINES.OUT")
        vertexLabGWFile = os.path.join(dirPath, "bandstructure.xml")
49
        selfCorGWFile = os.path.join(dirPath, "SELFC.DAT")
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
50
        inputFile = os.path.join(dirPath, "input.xml")
51
52
        inputgw1File = os.path.join(dirPath, "input-gw.xml")
        inputgw2File = os.path.join(dirPath, "input_gw.xml")
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
53
54

        if os.path.exists(inputFile):
55
56
57
58
59
60
61
62
63
64
65
66
            if os.path.exists(inputgw1File):
                inputgwFile = inputgw1File
            elif os.path.exists(inputgw2File):
                inputgwFile = inputgw2File
            else:
                inputgwFile = inputFile
        elif os.path.exists(inputgw1File):
            inputgwFile = inputgw1File
        elif os.path.exists(inputgw2File):
            inputgwFile = inputgw2File

        if os.path.exists(inputgwFile):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
67
68
            selfGWSetGIndex = backend.openSection("section_method")
            backend.addValue('electronic_structure_method', "G0W0")
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
69
            backend.addValue('gw_starting_point', xcName)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
70
71
72
73
74
            if dftMethodSectionGindex is not None:
                m2mGindex = backend.openNonOverlappingSection("section_method_to_method_refs")
                backend.addValue("method_to_method_ref", dftMethodSectionGindex)
                backend.addValue("method_to_method_kind", "starting_point")
                backend.closeNonOverlappingSection("section_method_to_method_refs")
75
76
                with open(inputgwFile) as f:
                    exciting_parser_input.parseInput(f, backend, self.gmaxvr)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
77
            backend.closeSection("section_method",selfGWSetGIndex)
78

79
80
81
82
83
84
85
86
87
88
89
90
        if os.path.exists(vertexGWFile):
            with open(vertexGWFile) as g:
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    s = s.split()
                    if len(s) > 0:
                        if not self.vertexDist:
                            self.vertexDist.append(float(s[0]))
                        elif float(s[0]) != self.vertexDist[-1]:
                            self.vertexDist.append(float(s[0]))
91
                self.vertexNum = len(self.vertexDist)-1
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
92

93
94
95
96
97
98
99
100
101
102
        if os.path.exists(vertexLabGWFile):
            with open(vertexLabGWFile) as g:
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    s = s.split()
                    if s[0] == "<vertex":
                        f = s[4].split("\"")
                        self.vertexLabels.append(f[1])
103

104
        if os.path.exists(eigvalGWFile):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
105
            eigvalGWGIndex = backend.openSection("section_eigenvalues")
106
107
            with open(eigvalGWFile) as g:
                qpGWKpoint=[]
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
108
                Vxc = [[],[]]
109
110
111
112
113
114
115
116
117
118
119
120
                Sx = [[],[]]
                Sc = [[],[]]
                qpE = [[],[]]
                Znk = [[],[]]
                fromH = unit_conversion.convert_unit_function("hartree", "J")
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    if "k-point" in s.split():
                        qpGWKpoint.append([])
                        for i in range(0,2):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
121
                            Vxc[i].append([])
122
123
124
125
                            Sx[i].append([])
                            Sc[i].append([])
                            qpE[i].append([])
                            Znk[i].append([])
126
                        x,y,z = float(s.split()[3]),float(s.split()[4]),float(s.split()[5])
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
                        qpGWKpoint[-1].append(x)
                        qpGWKpoint[-1].append(y)
                        qpGWKpoint[-1].append(z)
                    else:
                        s=s.split()
                        if len(s) == 0:
                            continue
                        else:
                            try: not int(s[0])
                            except ValueError:
                                continue
                            if self.spinTreat:
                                pass
                            else:
                                for i in range(0,2):
142
143
144
                                    Sx[i][-1].append(fromH(float(s[4])))
                                    Sc[i][-1].append(fromH(float(s[5])))
                                    qpE[i][-1].append(fromH(float(s[3])))
145
                                    Znk[i][-1].append(float(s[9]))
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
146
                                    Vxc[i][-1].append(fromH(float(s[6])))
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
147
148
149
150
        backend.addValue("eigenvalues_kpoints", qpGWKpoint)
        backend.addValue("number_of_eigenvalues", len(qpE[0]))
        backend.addValue("number_of_eigenvalues_kpoints", len(qpGWKpoint))
        backend.addValue("eigenvalues_values", qpE)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
151
        backend.addValue("gw_qp_linearization_prefactor", Znk)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
152
        backend.closeSection("section_eigenvalues",eigvalGWGIndex)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
153
154
155
        backend.addValue("gw_self_energy_x", Sx)
        backend.addValue("gw_self_energy_c", Sc)
        backend.addValue("gw_xc_potential", Vxc)
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
156
157

        ####################DOS######################
158

Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
159
        if os.path.exists(dosGWFile):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
160
161
            dosGWGIndex = backend.openSection("section_dos")
            ha_per_joule = unit_conversion.convert_unit(1, "hartree", "J")
162
#            bohr_cube_to_m_cube = unit_conversion.convert_unit(1, "bohr^3", "m^3")
163
            fromH = unit_conversion.convert_unit_function("hartree", "J")
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
164
165
166
167
168
169
170
171
            with open(dosGWFile) as g:
                dosValues = [[],[]]
                dosEnergies = []
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    s = s.split()
172
                    ene, value = fromH(float(s[0])), float(s[1])/ha_per_joule
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
173
174
175
176
177
178
                    dosEnergies.append(ene)
                    if not self.spinTreat:
                        for i in range(0,2):
                            dosValues[i].append(value)
                    else:
                        pass
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
179
180
181
182
183
184
            backend.addValue("dos_energies", dosEnergies)
            backend.addValue("dos_values", dosValues)
            backend.addValue("number_of_dos_values", len(dosEnergies))
            backend.closeSection("section_dos",dosGWGIndex)        

        ##################BANDSTRUCTURE#####################
185

186
        if os.path.exists(bandCarbGWFile):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
187
            bandGWGIndex = backend.openSection("section_k_band")
188
189
            fromH = unit_conversion.convert_unit_function("hartree", "J")

190
            with open(bandCarbGWFile) as g:
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
                bandEnergies = [[],[]]
                kpoint = []
                dist = []
                Kindex = [0]
                segmK = []
                segmLength = []
                bandEnergiesSegm = []
                bandGWBE = []
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    s = s.split()                
                    if not self.spinTreat:
                        if len(s) == 0:
                            for i in range(0,2):
                                bandEnergies[i].append([])
                        elif s[0] == "#":
                            for i in range(0,2):
                                bandEnergies[i].append([])
                                numBand = int(s[2])
                                numK = int(s[3])          
                        elif len(s) > 0:
                            for i in range(0,2):
                                bandEnergies[i][-1].append(fromH(float(s[6])))
                            if int(s[0]) == 1:
                                kpoint.append([])
                                dist.append(float(s[5]))
                                kpoint[-1].append([float(s[2]),float(s[3]),float(s[4])])
                    else:
                        pass

                for i in range(0,2):
                    bandEnergies[i].pop()

                for i in range(1,numK):
                    if dist[i] == dist[i-1]:
                        Kindex.append(i)
                Kindex.append(numK)
                for i in range(0,len(Kindex)-1): 
                    segmK.append(dist[Kindex[i]:Kindex[i+1]])

                for i in range(0,len(segmK)):
                    segmLength.append(len(segmK[i]))
                for i in range(0,2):
                    bandEnergiesSegm.append([])
                    for j in range(0,numBand):
                         bandEnergiesSegm[i].append([])
                         for k in range (0,len(Kindex)-1):
                             bandEnergiesSegm[i][j].append(bandEnergies[i][j][Kindex[k]:Kindex[k+1]])
            for i in range(0,len(Kindex)-1):
                bandGWBE.append([])
                for j in range(0,2):
                    bandGWBE[i].append([])
                    for k in range(0,segmLength[i]):
                        bandGWBE[i][j].append([])
                        for l in range(0,numBand):
                            bandGWBE[i][j][-1].append(bandEnergiesSegm[j][l][i][k])

250
            for i in range(0,len(Kindex)-1):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
251
252
253
                bandGWSegmGIndex = backend.openSection("section_k_band_segment")
                backend.addValue("band_energies", bandGWBE[i])
                backend.closeSection("section_k_band_segment",bandGWSegmGIndex)
254

Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
255
            backend.closeSection("section_k_band",bandGWGIndex)
256
257

        if os.path.exists(bandBorGWFile) and not os.path.exists(bandCarbGWFile):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
258
            bandGWGIndex = backend.openSection("section_k_band")
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
            fromH = unit_conversion.convert_unit_function("hartree", "J")
            with open(bandBorGWFile) as g:
                bandEnergies = [[[]],[[]]]
                kappa = [[[]],[[]]]
                dist1 = [[]]
                Kindex = [0]
                segmK = []
                segmLength = []
                bandEnergiesSegm = []
                bandGWBE = []
                while 1:
                    s = g.readline()
                    if not s: break
                    s = s.strip()
                    s = s.split()
                    if not self.spinTreat:
                        if len(s) == 0:
                            for i in range(0,2):
                                bandEnergies[i].append([])
                                kappa[i].append([])
                            dist1.append([])
                        elif len(s) > 0:
                            for i in range(0,2):
282
                                bandEnergies[i][-1].append(fromH(float(s[1])))
283
284
285
286
287
288
289
290
291
292
293
                                kappa[i][-1].append(float(s[0]))
                            dist1[-1].append(float(s[0]))
                        numK = len(kappa[0][0])
                        for i in kappa[0][0]:
                            if kappa[0][0].count(i) > 1:
                                kappa[0][0].remove(i)
                    else:
                        pass
                for i in range(0,2):
                    bandEnergies[i].pop()
                numBand = len(bandEnergies[0])
294
                for i in range(1,numK + self.vertexNum-1):
295
296
                    if dist1[0][i] == dist1[0][i-1]:
                        Kindex.append(i)
297
                Kindex.append(numK + self.vertexNum-1)
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
                for i in range(0,len(Kindex)-1):
                    segmK.append(dist1[0][Kindex[i]:Kindex[i+1]])

                for i in range(0,len(segmK)):
                    segmLength.append(len(segmK[i]))
                for i in range(0,2):
                    bandEnergiesSegm.append([])
                    for j in range(0,numBand):
                         bandEnergiesSegm[i].append([])
                         for k in range (0,len(Kindex)-1):
                             bandEnergiesSegm[i][j].append(bandEnergies[i][j][Kindex[k]:Kindex[k+1]])
            for i in range(0,len(Kindex)-1):
                bandGWBE.append([])
                for j in range(0,2):
                    bandGWBE[i].append([])
                    for k in range(0,segmLength[i]):
                        bandGWBE[i][j].append([])
                        for l in range(0,numBand):
                            bandGWBE[i][j][-1].append(bandEnergiesSegm[j][l][i][k])

318
            for i in range(0,len(Kindex)-1):
Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
319
320
321
                bandGWSegmGIndex = backend.openSection("section_k_band_segment")
                backend.addValue("band_energies", bandGWBE[i])
                backend.closeSection("section_k_band_segment",bandGWSegmGIndex)
322

Pardini, Lorenzo (lopa)'s avatar
Pardini, Lorenzo (lopa) committed
323
324
325
            backend.closeSection("section_k_band",bandGWGIndex)
        backend.closeNonOverlappingSection("section_single_configuration_calculation")
#        logging.error("done GW onClose_section_single_configuration_calculation")
326