1#!/usr/local/bin/python3.8
2
3# B. Malone, 10/18/09
4# Usage: qp_shifts.py sigma_hp.log
5
6# Script takes a sigma_hp.log as input and produces two files, 'cond' and
7# 'val'. These files contain the quasiparticle shifts for each band. 'cond'
8# contains all conduction bands and 'val' valence bands.
9
10# Output format:
11# E_lda E_qp-E_lda
12#       .
13#       .
14#       .
15
16# You can use the files to determine the QP scissor shift parameters
17# ecs, evs, ecdel, evdel, ev0, and ec0. ec0, ev0, ecs, and evs can simply be read
18# off from the files as the lowest conduction band energy, highest valence
19# band energy, the shift of the lowest conduction band, and the shift of the
20# highest valence band. To determine ecdel and evdel (the slopes) you will
21# need to fit a line to the points with a program like gnuplot.
22
23# For example:
24# gnuplot > f(x)=a*x+b;a=1;b=1;     # a,b initial guesses, doesn't matter
25# gnuplot > fit f(x) 'cond' via a,b
26# gnuplot > [output giving you the fitted parameters]
27# gnuplot > plot f(x), 'cond'       # let you see the fitted line along
28                                    # with the data points
29
30
31import sys
32
33
34try:
35    infilename=sys.argv[1]
36except:
37    sys.exit('Error: The sigma_hp.log file should be the first argument')
38
39print('Which band is the top of the valence band?')
40valtop=int(input())
41
42infile=open(infilename,'r')
43
44condlist=[]
45vallist=[]
46
47flag=0
48kpcounter=0
49
50
51while True:
52    line=infile.readline()
53    if line=='': # EOF
54        break
55    if line.find('k = ')!=-1 and line.find('spin =')!=-1: # found kpoint
56        kpcounter=kpcounter+1
57        bandcounter=0
58        for i in range(1,3):   # read in two lines
59            infile.readline()
60        while True:        # reading in bands
61            line=infile.readline()
62            if line.find('.')!=-1: # there are energy values (not blank)
63                bandcounter=bandcounter+1
64                theline=line.split()
65                evalue=float(theline[1])
66                GWvalue=float(theline[9])
67                diff=GWvalue-evalue
68                addition=(evalue,diff)
69                if int(theline[0])<=valtop:
70                    vallist.append(addition)
71                else:
72                    condlist.append(addition)
73            else:          # done reading all bands for this kpoint
74                break
75
76
77infile.close()
78
79print 'There are',bandcounter,'bands and', kpcounter, 'kpoints in the file'
80
81
82vallist.sort()
83condlist.sort()
84
85# now write to outfiles
86
87outfile1=open('cond','w')
88outfile2=open('val','w')
89
90for i in condlist:
91    print >>outfile1, i[0],i[1]
92
93for i in vallist:
94    print >>outfile2, i[0],i[1]
95
96outfile1.close()
97outfile2.close()
98
99print('Done, info written to \'cond\' and \'val\'')
100