Version 1.3. Added plotting utility

BUILD

@@ -6,7 +6,8 @@ source files. They can be bundled together into an
 executable archive by running the build.sh script.
 Once vasp_unfold is generated by build sript, it can
 be run independently and you can place it wherever
-it fits you best to use it.
+it fits you best to use it. Since version 1.3 this
+script also builds plotting utility fatplot.
 
 If Python is invoked in a non-standard way on your
 system, please change the PYTHON_COMMAND variable

CHANGELOG

@@ -0,0 +1,49 @@
+Jan 23th 2015
+=============
+
+v1.0
+
+* Added the possiblity to parse the PROCAR file from the
+  non-collinear spin-polarized calculation.
+
+* Added the possibility to unfold both spin up and spin
+  down components.
+
+
+Jan 28th 2015
+=============
+
+v1.1
+
+* Refactored code into separate files. App can be packed 
+  into single executable with pack.sh script. 
+
+* Spin polarized and non-collinear PROCARS are now fully 
+  parsed.
+
+
+Aug 18th 2015
+=============
+
+v1.2
+
+* Can handle vacancies and excess atoms not. One to one
+  mapping check is now optional.
+
+
+Aug 31st 2015
+=============
+
+v1.3
+
+* Changed the way orbital weights are unfolded. Orbital
+  weights from the non-collinear spin-polarized calculation
+  are now fully unfolded.
+
+* Added a small plotting utility to plot fatbands.
+
+* File structure reorganized again.
+
+* Added extended error message regarding broken PROCAR
+  formatting.
+

README.md

@@ -4,6 +4,8 @@ Here you can download vasp_unfold, a Python script which you can use to unfold t
 
 You can use the code in whatever way you see fit, but if it is used to produce the data for the publication, please cite [Ref. 1](#ref_1).
 
+In addition, a small utility called fatplot is provided to quickly plot the bandstructures contained in the VASPs PROCAR files.
+
 The code is given as is, ie. there is no guarantee that it will work. However, if you have some problems, or the code does not behave in a way you expect it to, I encourage you to report the problem to my [e-mail](mailto:[email protected]) and I will try to fix it as soon as possible.
 
 ## Installation
@@ -12,6 +14,7 @@ Dependencies are
 
 * [Python](http://www.python.org/downloads/) ofcourse. The script requires Python2.7 or greater
 * [NumPy](http://www.scipy.org/scipylib/download.html) library for Python
+* [matplotlib](http://matplotlib.org/downloads.html) library for Python (only needed if fatplot is used)
 
 No special installation is required. Just place it wherever it suits you and run it. 
 
@@ -51,9 +54,7 @@ The unfolded bandstructures will be located in PROCAR.irrep.0 file. In case --al
 
 **NOTE 1**: no whitespace is allowed in the fractional translation generator specification. Also, the components can be either 0, or 1/N, where N is an integer. Floating point values are not allowed. 
 
-**NOTE 2**: in case of spin polarized non-collinear calculation only spin-up and spin-down orbital weight totals will be unfolded. Orbital weight components for x,y and z component of spin won't be unfolded.
-
-**NOTE 3**: do not enable --check-mapping flag if your structure has vacancies or excess atoms, since in this case fractional translations do not map every atom onto some other atom.
+**NOTE 2**: do not enable --check-mapping flag if your structure has vacancies or excess atoms, since in this case fractional translations do not map every atom onto some other atom.
 
 ## Resolving the issues with the code
 

build.sh

@@ -1,10 +1,12 @@
 #! /bin/bash
 
 # This scripts packs the Python source files 
-# into an executable ZIP archive
+# into an executable ZIP archive. This is done
+# separatly for the vasp_unfold and 
 
-# Name for the executable
+# Names for the executables
 VASP_UNFOLD="vasp_unfold"
+PLOT="fatplot"
 
 # Command used to invoke Python
 PYTHON_COMMAND="python"
@@ -13,21 +15,26 @@ PYTHON_COMMAND="python"
 ENV_COMMAND="/usr/bin/env"
 
 
-SRC_FILES="__main__.py parse.py unfolding.py utils.py write.py"
+SRC_FILES="__main__.py parse.py unfolding.py utils.py write.py errors.py"
+PLOT_SRC_FILES="__main__.py"
 
 # Change into source directory
 cd src
 
-# Zip source files
-zip ../${VASP_UNFOLD}.zip ${SRC_FILES}
+# Zip the source files
+cd unfolding; zip ../../${VASP_UNFOLD}.zip ${SRC_FILES}; cd ..;
+cd plot;      zip ../../${PLOT}.zip ${PLOT_SRC_FILES};   cd ..;
 
 cd ..
 
 # Prepend shebang to the zip archive
 echo "#! ${ENV_COMMAND} ${PYTHON_COMMAND}" | cat - ${VASP_UNFOLD}.zip > ${VASP_UNFOLD}
+echo "#! ${ENV_COMMAND} ${PYTHON_COMMAND}" | cat - ${PLOT}.zip > ${PLOT}
 
 # Remove zip archive
 rm ${VASP_UNFOLD}.zip
+rm ${PLOT}.zip
 
 # Make it executable
 chmod +x ${VASP_UNFOLD}
+chmod +x ${PLOT}

src/plot/__main__.py

@@ -0,0 +1,159 @@
+#===========================================================
+#
+#  PROJECT: vasp_unfold/fatplot
+#  FILE:    __main__.py
+#  AUTHOR:  Milan Tomic
+#  EMAIL:   [email protected]
+#  VERSION: 1.3
+#  DATE:    Aug 31st 2015
+#
+#===========================================================
+
+
+import os
+import numpy as np
+import argparse
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plot
+
+# Handle color arguments
+def color(string):
+    try:
+        # Try to convert to tuple
+        return eval(string)
+    except:
+        # If cannot convert to tuple just
+        # return the original string
+        return string
+
+desc_str = '''Simple program used to quickly plot the orbital weights
+of the band structure contained in the specified PROCAR file.
+'''
+
+parser = argparse.ArgumentParser(prog='fatplot', description = desc_str)
+
+parser.add_argument('procar', type=str, help='POSCAR file')
+parser.add_argument('output', type=str, help='Filename for the plot. It accepts '
+                    'all image formats supported by matplotlib.')
+
+parser.add_argument('--figsize', type=eval, default=(6, 4),
+                    help='Figure size in inches formatted as width,height '
+                    '(no spaces allowed). Default is 6,4')
+parser.add_argument('--dpi', type=float, default=300,
+                    help='Plot resolution in dots per inch. Default is 300.')
+parser.add_argument('--marker', type=str, default='o', 
+                    help='Marker for the fatband plot. Default is o.')
+parser.add_argument('--markersize', type=float, default=20,
+                    help='Marker size. Default is 20.')
+parser.add_argument('--color', type=color, default='"steelblue"',
+                    help='Color for the marker. It accepts any color specification '
+                    'accepted by matplotlib. Color specified as r,g,b tuple should ' 
+                    'not have any spaces.')
+parser.add_argument('--elim', type=eval, default=(1,-1),
+                    help='Energy range for the plot, specified as emin,emax '
+                    '(no spaces allowed). Default is entire band range.')
+parser.add_argument('--efermi', type=float, default=0,
+                    help='Fermi energy. Default is 0.')
+parser.add_argument('--pow', type=float, default=1,
+                    help='Raise orbital weights to the specified integer power. '
+                    'Powers larger than 1 help to filter out the ghost bands '
+                    'in the unfolded band structures.')
+
+args = parser.parse_args()
+
+
+# Commands to extract the needed information from the PROCAR file
+grep_npoints = 'grep -m 1 "of k\-points" {0} | tr -s " " | cut -d" " -f4'
+grep_kpoints = 'grep -E "^ k\-point " {0} | tr -s " " | cut -d" " -f5-7'
+grep_bands = 'grep -E "^band" {0} | tr -s " " |  cut -d" " -f5'
+grep_weights = 'grep -E "^tot" {0} | tr -s " " |  cut -d" " -f11'
+
+# Extract the number of k-points
+npoints = int(os.popen(grep_npoints.format(args.procar)).read())
+
+# Extract the k-points
+kpoints = os.popen(grep_kpoints.format(args.procar))
+kpoints = np.fromfile(kpoints, count=3*npoints, dtype=float, sep=' ')
+
+# Extract the band energies
+bands = os.popen(grep_bands.format(args.procar))
+bands = np.fromfile(bands, count=-1, dtype=float, sep=' ')
+
+# Extract the total orbital weights for each band
+weights = os.popen(grep_weights.format(args.procar))
+weights = np.fromfile(weights, count=-1, dtype=float, sep=' ')
+
+# Figure out the number of bands
+nbands = len(bands)/npoints
+
+# Figure out the number of orbital blocks per band
+# 1 for collinear calculation, 4 for non-collinear
+# In non-collinear case we just need the first one
+wdim = len(weights)/(npoints*nbands)
+
+# Reshape the arrays into their proper shapes
+kpoints = kpoints.reshape((npoints, 3))
+bands = bands.reshape((npoints, nbands))-args.efermi
+weights = weights[::wdim].reshape((npoints, nbands))
+
+# Raise the weights to the specified power
+if args.pow != 1:
+    weights = np.power(weights, args.pow)
+
+
+# Try to guess where the high symmetry points 
+# are by looking for the repeated k-points
+dk = np.zeros((npoints, 3), float)
+
+# Displacement between i+1st and ith k-point
+dk[1:] = kpoints[1:]-kpoints[:-1]
+
+# Get the magnitude of displacements
+dk = np.sqrt(np.sum(dk*dk, axis=1))
+
+# Locate the high-symmetry points
+ipoint = np.where(dk < 1e-6)[0]
+
+# Generate plot's x-axis from the adjacent k-point
+# displacement magnitudes
+x = np.cumsum(dk)
+
+# Get the plot's x-coordinates for high-symmetry k-points
+xsym = x[ipoint]
+
+# Figure out the maximum possible energy boundaries 
+# based on the energy extend of all bands
+e_low = np.min(bands)-1.0
+e_high = np.max(bands)+1.0
+
+plot.figure(figsize=args.figsize)
+
+# Plot horizontal line for the Fermi level
+plot.plot(x[[0,-1]], [0, 0], color='gray', zorder=-1)
+
+# Plot vertical lines for the high-symmetry points
+for xi in xsym:
+    plot.plot([xi, xi], [e_low, e_high], color='gray', zorder=-1)
+
+# Plot the weights
+for bi, wi in zip(bands.T, weights.T):
+    plot.scatter(x, bi, s=args.markersize*wi, marker=args.marker, 
+                 color=args.color, lw=0)
+
+# Fix x and y-axis boundaries
+plot.xlim(x[0], x[-1])
+
+if args.elim[0] < args.elim[1]:
+    plot.ylim(*args.elim)
+else:
+    plot.ylim(e_low, e_high)
+
+# Set x-axis ticks
+plot.xticks(xsym, ['']*len(xsym))
+
+# Save the figure
+plot.savefig(args.output, dpi=args.dpi, bbox_inches='tight')
+
+

src/__main__.py → src/unfolding/__main__.py

@@ -1,11 +1,15 @@
 #! /usr/bin/env python
 
 
-# vasp_unfold
+#============================================================================
 #
-# version 1.2
+#  PROJECT: vasp_unfold
+#  FILE:    __main__.py
+#  AUTHOR:  Milan Tomic
+#  EMAIL:   [email protected]
+#  VERSION: 1.3
+#  DATE:    Aug 31st 2015
 #
-# Author: Tomic Milan
 #
 # vasp_unfold is a code whose purpose is to perform the unfolding
 # of band structures calculated with VASP. It is based on the method
@@ -15,6 +19,8 @@
 # This code is provided as is. There is no guarantee that it will work.
 # However, if you encounter errors or unexpected behavor, please report
 # it to [email protected] and I will try to improve it.
+#
+#============================================================================
 
 
 import numpy as np
@@ -24,6 +30,7 @@
 from unfolding import build_translations, build_operators, build_projectors
 from parse import parse_poscar, parse_procar
 from write import write_procar
+import errors
 
 def main():
     desc_str = 'Unfold bands calculated by VASP. For this, phase '\
@@ -95,15 +102,14 @@ def main():
     try:
         data = parse_procar(args.procar)
     except:
-        post_error('Unable to parse the input PROCAR file. Please '
-            'check if the file exists and is formatted properly.')
-
-    phases = np.copy(data[-1])
+        post_error(errors.poscar_parse_error)
 
-    if phases is None:
+    if data[-1] is None:
         post_error('Phase information has to be present in the PROCAR '
             'file. Please repeat the calculation with LORBIT=12.')
 
+    phases = np.copy(data[-1])
+
     norbs = phases.shape[1]/len(spos)
 
     projs = build_projectors(irreps, ops, norbs)
@@ -128,11 +134,13 @@ def main():
                     'the same crystal structure?')
 
         # We update total absolute weights to correspond to
-        # unfolded phases.
-        # WARNING: In non-collinear calculations, mx, my and
-        # mz weight blocks are not unfolded
-        data[-2][:,:,:,0,:] = np.abs(data[-1])
+        # unfolded phases (by multiplying them by the magnitude
+        # ratio of unfolded and folded phases 
+        phase_ratio = np.abs(data[-1])/(np.abs(phases)+1e-4)
 
+        for idim in xrange(data[-2].shape[3]):
+            data[-2][:,:,:,idim,:] *= phase_ratio
+
         write_procar('{0}.irrep.{1}'.format(output, i), *data)