Handled new PROCAR version and fixed Selective Dynamics bug

CHANGELOG

@@ -79,4 +79,13 @@ v1.33
 
 v1.34
 
-* Fixed bug with numpy fraction issues.
+* Fixed bug with numpy fraction issues.
+
+
+  December 12th 2017
+=============
+
+v1.4
+
+* Added handling of PROCAR files for VASP 5.4.4
+* Fixed problem where POSCARS with Selective Dynamics were not handled

README.md

@@ -40,6 +40,7 @@ where the parameters are
 --eps            Numerical tolerance for position discrimination
 --all-irreps     Write all irreps from the unfolding
 --check-mapping  Verify if fractional translations map atoms one-to-one
+--vasp-version   Which version of VASP was used to produce the PROCAR file
 poscar           Location of POSCAR file
 procar           Location of PROCAR file
 ```

src/unfolding/__main__.py

@@ -7,8 +7,8 @@
 #  FILE:    __main__.py
 #  AUTHOR:  Milan Tomic
 #  EMAIL:   [email protected]
-#  VERSION: 1.32
-#  DATE:    Apr 25th 2017
+#  VERSION: 1.4
+#  DATE:    December 12th 2017
 #
 #
 # vasp_unfold is a code whose purpose is to perform the unfolding
@@ -26,7 +26,7 @@
 import numpy as np
 import argparse
 import sys
-from utils import post_error, translation
+from utils import post_error, translation, version
 from unfolding import build_translations, build_operators, build_projectors
 from parse import parse_poscar, parse_procar
 from write import write_procar
@@ -84,7 +84,11 @@ def main():
                         'one-to-one, ie. map every atom on exactly one other atom '
                         'in the unit cell. This MUST not be enabled for the cases '
                         'where vacancies or excess atoms are present.')
-
+
+    parser.add_argument('--vasp-version', type=version, default='5.2.2', help='Version of VASP'
+                        'that produced the PROCAR file. All versions prior to 5.4.4'
+                        'use the same format. Since 5.4.4 there is a new format.')
+
     args = parser.parse_args()
 
     tgens = args.tgen
@@ -95,19 +99,19 @@ def main():
         cell, spos, symbols = parse_poscar(args.poscar)
     except:
         post_error('Unable to parse the input POSCAR file. Please '
-            'check if the file exists and is formatted properly')
+            'check if the file exists and is formatted properly.', True)
 
     ops = build_operators(spos, trans, args.check_mapping, args.eps)
 
     try:
-        data = parse_procar(args.procar)
-    except:
-        post_error(errors.poscar_parse_error)
+        data = parse_procar(args.procar, args.vasp_version)
+    except Exception as exc:
+        post_error(errors.poscar_parse_error, True)
 
     if data[-1] is None:
         post_error('Phase information has to be present in the PROCAR '
-            'file. Please repeat the calculation with LORBIT=12.')
-
+            'file. Please repeat the calculation with LORBIT=12.', True)
+    
     phases = np.copy(data[-1])
 
     norbs = phases.shape[1]/len(spos)

src/unfolding/errors.py

@@ -4,8 +4,8 @@
 #  FILE:    __errors__.py
 #  AUTHOR:  Milan Tomic
 #  EMAIL:   [email protected]
-#  VERSION: 1.32
-#  DATE:    Apr 25th 2017
+#  VERSION: 1.4
+#  DATE:    December 12th 2017
 #
 #===========================================================
 

src/unfolding/parse.py

@@ -4,8 +4,8 @@
 #  FILE:    parse.py
 #  AUTHOR:  Milan Tomic
 #  EMAIL:   [email protected]
-#  VERSION: 1.32
-#  DATE:    Apr 25th 2017
+#  VERSION: 1.4
+#  DATE:    December 12th 2017
 #
 #===========================================================
 
@@ -55,12 +55,15 @@ def parse_poscar(filename):
     # Cartesian or fractional coordinates?
     ctype = gl.readline()[0].lower()
 
+    if ctype == 's':
+        ctype = gl.readline()[0].lower()
+
     if ctype == 'c':
         mult = np.linalg.inv(cell)
     elif ctype == 'd':
         mult = np.eye(3)
     else:
-        post_error('"{0}" is unknown POSCAR option'.format(plines[7].strip()))
+        post_error('"{0}" is unknown POSCAR option'.format(ctype))
 
     # Allocate storage for positions
     spos = np.zeros((len(symbols), 3))
@@ -76,7 +79,7 @@ def parse_poscar(filename):
     return cell, spos, symbols
 
 
-def parse_procar(filename):
+def parse_procar(filename, vasp_version):
     '''This function parses a PROCAR file. It returns a tuple
     consisting of following elements:
     
@@ -157,6 +160,7 @@ def get_absweights(i, j, s):
         for k in xrange(dim):
             # Fetch entire orbital weight block
             data = np.fromfile(gl, sep=" ", count=nions*(norbs+2))
+
             # Cast it into tabular shape
             data = data.reshape((nions, norbs+2))
 
@@ -170,25 +174,48 @@ def get_absweights(i, j, s):
     if '+ phase' in header_1:
         # Allocate storage for phases
         phases = np.zeros((npoints, nions*norbs, nbands, 2), complex)
+
+        if vasp_version < (5, 4, 4):
+            # Declare nested function that handles 
+            # parsing of complex weights
+            def get_weights(i, j, s):
+                # Read abs values of weights
+                get_absweights(i, j, s)
+
+                # Skip line with orbital names
+                gl.readline()
+
+                # Fetch entire phase block
+                data = np.fromfile(gl, sep=" ", count=2*nions*(norbs+1))
+                # Cast it into tabular shape
+                data = data.reshape((2*nions, norbs+1))
 
-        # Declare nested function that handles 
-        # parsing of complex weights
-        def get_weights(i, j, s):
-            # Read abs values of weights
-            get_absweights(i, j, s)
-
-            # Skip line with orbital names
-            gl.readline()
-
-            # Fetch entire phase block
-            data = np.fromfile(gl, sep=" ", count=2*nions*(norbs+1))
-            # Cast it into tabular shape
-            data = data.reshape((2*nions, norbs+1))
+                # Discard first column and store real and imaginary
+                # parts respectively
+                phases[i,:,j,s] =  data[::2,1:].flatten()
+                phases[i,:,j,s] += 1j*data[1::2,1:].flatten()
+        else:
+            # Declare nested function that handles 
+            # parsing of complex weights
+            def get_weights(i, j, s):
+                # Read abs values of weights
+                get_absweights(i, j, s)
+
+                # Skip line with orbital names
+                gl.readline()
+
+                # Fetch entire phase block
+                data = np.fromfile(gl, sep=" ", count=nions*(2*norbs+2))
+                # Cast it into tabular shape
+                data = data.reshape((nions, 2*norbs+2))
 
-            # Discard first column and store real and imaginary
-            # parts respectively
-            phases[i,:,j,s] =  data[::2,1:].flatten()
-            phases[i,:,j,s] += 1j*data[1::2,1:].flatten()
+                # Discard first column and store real and imaginary
+                # parts respectively
+                phases[i,:,j,s] =  data[:,1:-1:2].flatten()
+                phases[i,:,j,s] += 1j*data[:,2:-1:2].flatten()
+
+                # Skip line with charges
+                gl.readline()
     else:
         # Phases are None in this case
         phases = None
@@ -210,7 +237,7 @@ def get_weights(i, j, s):
         for j in xrange(nbands):
             # Parse band energy
             band_line = gl.readline().split()
-
+            
             bands[i, j, 0] = float(band_line[4])
             occupancies[i, j, 0] = float(band_line[-1])
 

src/unfolding/utils.py

@@ -4,21 +4,35 @@
 #  FILE:    utils.py
 #  AUTHOR:  Milan Tomic
 #  EMAIL:   [email protected]
-#  VERSION: 1.34
-#  DATE:    May 31st 2017
+#  VERSION: 1.4
+#  DATE:    December 12th 2017
 #
 #===========================================================
 
 import numpy as np
 import argparse
 import sys
 import fractions
+import traceback
 
 
-def post_error(error_info):
+def post_error(error_info, show_traceback=False):
     '''Write error message to sderr and exit program
     '''
-    sys.stderr.write('\nError: '+error_info+'\n\n')
+    message = '\nError: '+error_info+'\n\n'
+
+    if show_traceback:
+        # Show exception traceback
+        tb_msg = traceback.format_exc()
+
+        max_line_len = max(len(line) for line in tb_msg.split('\n'))
+
+        message += 'Operation failed due to the following exception:\n'
+        message += '='*max_line_len
+        message += '\n' + traceback.format_exc()
+        message += '='*max_line_len
+
+    sys.stderr.write(message)
 
     exit()
 
@@ -90,8 +104,18 @@ def translation(tstring):
     except:
         post_error('Unable to parse string: "{0}". Check help for valid '
                    'translation generator specification'.format(tstring))
-
 
+def version(vstring):
+    '''Parse string containing version information.
+    Valid form has dot separated digits. Returns tuple
+    of integers which are to be lexicographically compared.
+    '''
+    try:
+        return tuple(int(d) for d in vstring.split('.'))
+    except:
+        post_error('Unable to parse string: "{0}". The valid version '
+                   'is composed of dot separated digists'.format(vstring))
+
 def lcm(a, b):
     '''Return lowest common multiple.'''
     return a * b // fractions.gcd(a, b)