source: ESPRESSO/branches/shirley_QE5.0/PP/tools/sum_states.py @ 441

Last change on this file since 441 was 441, checked in by davegp, 5 years ago

Preliminary commit of the official svn repo for QE into our shirley_QE5.0 branch,
code exported from:
http://qeforge.qe-forge.org/svn/q-e/trunk/espresso
revision 9624

* This has not been compiled yet *

The directory structure of SHIRLEY has been modified to reflect the new structure:
CODE
CODE/src
CODE/examples
...

I am committing this to see if anything is broken and to reduce the svn status output

davegp

  • Property svn:executable set to *
File size: 8.0 KB
Line 
1#! /usr/bin/python
2
3######  SUM STATES #######
4# Python script for summing and ploting the data from the Density Of States
5# files obtained from projwfc.x. It can sum also k-solved dos, and make a plot
6# with mathplotlib (if not available, gnuplot, if not avaible, print to file)
7# if there is not X11 forwarding, plots in terminal.
8# It does something very similar to sumpdos.f90, but with
9# some extra features (use "-h" option).
10#
11# it takes two different inputs, the first one is the pw.x output
12# ("-o" option), which is used for parsing the Fermi energy for fitting
13# the PDOS curve to the right energy. The other files are the pDOS files
14# ("-s" option), that can be given with shell syntax, i.e.
15# pdos_atm*Fe*wfc*d* for summing all the d orbitals of Fe.
16# It can also handle k solved dos files.
17#
18# One of the most useful feature, compared to the sumpdos.x, is the
19# fact that it also builds the picture directly, so it can be directly
20# visualized and exported for inclusion in a document.
21# It uses mathplotlib for plotting, but if no mathplotlib is found in
22# the $PYTHONPATH, it tries to use gnuplot, if no gnuplot available,
23# dumps the output data to a file.
24# In the that no X11 forwarding is available (i.e. ssh to the cluster),
25# it shows a rough graph in the terminal, so we get an idea of the shape
26# of the results.
27#
28# Example of usage:
29# cd ....../espresso-5.0/PP/examples/example02/results/
30# ../../../src/sum_states.py -o ni.dos.out -s
31# ni.pdos_atm#1\(Ni\)_wfc#2\(d\) -t "Example PP/02" -xr -6 2
32#
33#
34# The procedure for obtaining the DOS files is explained
35# i.e. in (espresso-dir)/PP/examples/example02/
36#
37# Author: Dr. Julen Larrucea
38#         University of Bremen,
39#         Bremen Centre for Computational Materials Science, HMI Group
40#         julenl [at] gmail.com  or larrucea [at] hmi.uni-bremen.de
41#
42#  This file is distributed under the terms of the GNU General Public
43#  License. See the file `License'
44#  in the root directory of the present distribution,
45#  or http://www.gnu.org/copyleft/gpl.txt .
46#######################
47
48import sys
49import os
50import fnmatch
51import linecache
52
53# Some default variables
54version=0.2
55pwout=""
56selat="*"
57graphtitle=""
58min_x,max_x=-10,3
59min_y,max_y="",""
60output_file_name="sum_dos.out"
61prt="no"
62
63print " #### sum_states.py version "+str(version)+" #### " 
64
65
66# Check if X11, mathplotlib and gnuplot are available
67try:
68 os.popen("gnuplot -V").read()
69 prog_gnuplot="yes" # gnuplot is installed
70except:
71 prog_gnuplot="no"
72
73
74
75# Parse command line options
76if len(sys.argv)>1:
77 for i in sys.argv:
78  if i.startswith('-'):
79   option=i.split('-')[1]
80   if option=="o":
81     pwout= sys.argv[sys.argv.index('-o')+1]
82   if option=="s":
83    selat= sys.argv[sys.argv.index('-s')+1]
84   if option=="p":
85    prt="yes"
86    if len(sys.argv) > sys.argv.index('-p')+1: # if there is a name after "-p" take it as an output name
87     if sys.argv[sys.argv.index('-p')+1] != "-": # otherwise default name sum_dos.out
88        dos_out_name=sys.argv[sys.argv.index('-p')+1]
89   if option=="t":
90    graphtitle= sys.argv[sys.argv.index('-t')+1]
91   if option=="xr":
92    min_x,max_x= float(sys.argv[sys.argv.index('-xr')+1]),float(sys.argv[sys.argv.index('-xr')+2])
93   if option=="yr":
94    min_y,max_y= float(sys.argv[sys.argv.index('-yr')+1]),float(sys.argv[sys.argv.index('-yr')+2])
95   if option=="v":
96    print "sum_dos.py version: "+version
97    sys.exit()
98   if option=="h":
99    print '''
100    -o QE output file name (for grepping Fermi E)
101    -s Selection of atoms for summing the DOSes. "*" for all, *1*Fe*d* for first Fe atom " (def. "*")
102    -p Print output to a file and aditionaly provide an output name (def. no output and "sum_dos.out")
103    -t set title in the head of the graph
104    -xr set min and max x value for the axes in the graph
105    -yr set min and max y value for the axes in the graph
106    -h print this help
107    -v print version
108       Example: sum_states.py --s sys.pdos_atm#4\(Fe2\)_wfc#2\(d\) -t "Wustite LDA+U single Fe" -xr -9 4
109   '''
110    sys.exit()
111
112
113# Check for mathplotlib/gnuplot and import mpl if possible
114if len(os.popen('echo $DISPLAY').read()) > 1:
115  graphic_plot="yes"
116  try:
117   from pylab import *
118   mplplot="yes"
119   print "pylab imported"
120  except:
121   print "There is no mathplotlib installed. Using gnuplot."
122   mplplot="no"
123   prt="yes"
124else:
125  print "No X11. Trying to plot on terminal"
126  graphic_plot="no"
127  if prog_gnuplot=="no":
128   prt="yes"
129
130
131# if not specified, try to find the espresso output, in order to parse the Fermi energy
132if pwout == "":
133 for filen in filter(os.path.isfile, os.listdir('.')):
134  if "Program PWSCF" in linecache.getline(filen, 2):
135   print "Using " + filen + " as pw.x output. You can specify another one with the -o option."
136   pwout=filen
137
138# Parse Fermi energy from the pw.x output
139if pwout!="":
140 try:
141  os.popen("grep -a 'the Fermi energy is' "+pwout ).read()
142  fermi=float(os.popen("grep -a 'the Fermi energy is' "+pwout ).read().split()[4])
143  print "Fermi energy = ", fermi, "a.u."
144 except:
145  print "WARNING: No Fermi energy found. Using 0 e.V. instead"
146  fermi=0
147else:
148 print "WARNING: No pw.x output found. Using E Fermi = 0 e.V."
149 fermi=0
150
151# List of all DOS files to add
152dosfiles=[]
153for dfile in os.listdir('.'):
154   if fnmatch.fnmatch(dfile, selat):
155     dosfiles.append(dfile)
156if len(dosfiles)==0:
157 print "ERROR: Provide a (list of) valid DOS file(s)"
158 sys.exit()
159
160print "dosfiles list: ",
161for dosfile in dosfiles:
162  print dosfile,
163print ""
164
165# Check wetter we have k-solved DOS
166if open(dosfiles[0],'r').readline().split()[1]=="E":
167  ksolved="no"
168  print "no ksolved"
169elif open(dosfiles[0],'r').readline().split()[1]=="ik":
170  ksolved="yes"
171  print "ksolved"
172
173# Sum over all k-points and files
174mat=[]  # matrix with total sum of ldos
175for i in range(len(dosfiles)):
176 mati=[] # temporal matrix for each DOS file "i"
177 k=0
178 for line in open(dosfiles[i],'r'):
179  if len(line) > 10 and line.split()[0] != "#":
180
181   if ksolved=="no":
182      mati.append([float(line.split()[0]),float(line.split()[1]),float(line.split()[2])])
183
184   if ksolved=="yes":
185      ik = int(line.split()[0])
186      if ik > k:  #if it is a different k block
187         k=int(line.split()[0])
188         oldmat=[] # temporal matrix for each k-point
189      if ik == 1:
190         mati.append([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])]) # append: energy, ldosup, ldosdw
191      elif ik == k and k > 1:
192         oldmat.append([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])])
193      elif len(line) < 5 and k > 1:  #if blank line, sum k-frame to the total
194       for j in range(len(oldmat)): 
195         mati[j]=[mati[j][0],mati[j][1]+oldmat[j][1],mati[j][2]+oldmat[j][2]]
196
197 if mat == []: # if it is the first dos file, copy total matrix (mat) = the first dos files's data
198    mat=mati[:]
199 else:
200    for j in range(len(mati)): # if it is not the first file, sum values
201        mat[j]=[mat[j][0],mat[j][1]+mati[j][1],mat[j][2]+mati[j][2]] 
202
203
204print "...ploting..."
205
206
207if prt=="yes":
208 out=open(output_file_name,"w")
209x,y1,y2=[],[],[]   
210for i in mat:
211 x.append(i[0]-fermi)
212 y1.append(i[1])
213 y2.append(-i[2])
214 if prt=="yes":  # print to a file
215  print>>out, i[0]-fermi, i[1], i[2]
216if prt=="yes":
217 out.close()
218
219
220if graphic_plot=="yes":
221  # if there is matplotlib, generate a plot with it
222  if mplplot=="yes":
223    plot(x,y1,linewidth=1.0)
224    plot(x,y2,linewidth=1.0)
225    print min(y2),max(y1)
226    plt.title(graphtitle)
227    plt.xlabel('E (eV)')
228    plt.ylabel('States')
229    plt.grid(True)
230    plt.rcParams.update({'font.size': 22})
231    plt.fill(x,y1,color='0.8')
232    plt.fill(x,y2,color='0.9')
233    if min_x and max_x:
234     fromx,tox=min_x,max_x
235    plt.axis([fromx, tox, min(y2), max(y1)])
236    show()   
237  elif mplplot=="no" and prog_gnuplot=="yes":  # If no mathplotlib available, use gnuplot
238     os.system("echo \"plot '"+ output_file_name + "' using ($1-"+str(fermi)+"):2 w l, '' u ($1"+str(fermi)+"):3 w l\" | gnuplot -persist") 
239elif graphic_plot=="no":  # If no X forwarding available, show graph in terminal
240  if prog_gnuplot=="yes":
241     os.system("echo \"set terminal dumb; plot '"+ output_file_name + "' using ($1-"+str(fermi)+"):2 w l, '' u ($1-"+str(fermi)+"):3 w l\" | gnuplot -persist")
242
243
244
245
Note: See TracBrowser for help on using the repository browser.