Ionizing · August 1, 2024 01:19 · Dec 28, 2021 · Apr 23, 2021 · Jan 25, 2021 · Jan 24, 2021
diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -9,8 +9,8 @@
     - 1:56, Jan 25th, 2021
         - Add E-fermi level correction
         - More logging info
-    - 17:25, Apr 23th, 2021
-        - Change plotimage y range and turn on minor ticks
+    - 10:05, Dec 28th, 2021
+        - Fix E-fermi extraction error, return the last one if there are many 'E-fermi' items.
 '''
 
 from argparse import ArgumentParser
@@ -38,7 +38,7 @@ def get_efermi(outcar="OUTCAR"):
             logger.warning("OUTCAR file not found. E-fermi set to 0.0eV")
             return None
         txt = open(outcar).read()
-        efermi = re.search(r'E-fermi :\s*([-+]?[0-9]+[.]?[0-9]*([eE][-+]?[0-9]+)?)', txt).groups()[0]
+        efermi = re.findall(r'E-fermi :\s*([-+]?[0-9]+[.]?[0-9]*([eE][-+]?[0-9]+)?)', txt)[-1][0]
         logger.info("Found E-fermi = {}".format(efermi))
         return float(efermi)
 
@@ -108,4 +108,4 @@ def parse_cml_arguments():
         plt.title(args.title)
 
     logger.info("Saving to {}".format(args.output))
-    plt.savefig(args.output, dpi=args.dpi)
+    plt.savefig(args.output, dpi=args.dpi)
diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -9,6 +9,8 @@
     - 1:56, Jan 25th, 2021
         - Add E-fermi level correction
         - More logging info
+    - 17:25, Apr 23th, 2021
+        - Change plotimage y range and turn on minor ticks
 '''
 
 from argparse import ArgumentParser
@@ -99,10 +101,11 @@ def parse_cml_arguments():
     plt.ylabel('Potential(eV)')
     plt.grid(color='gray', ls='-.')
     plt.xlim(0, np.max(x))
-    plt.ylim(np.min(y)-0.5, np.max(y)+0.5)
+    plt.ylim(np.max(y)-2, np.max(y)+0.5)
+    plt.minorticks_on()
 
     if args.title:
         plt.title(args.title)
 
     logger.info("Saving to {}".format(args.output))
-    plt.savefig(args.output, dpi=args.dpi)
+    plt.savefig(args.output, dpi=args.dpi)
diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -33,7 +33,7 @@ def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat', outcar="OUTCAR"
     '''
     def get_efermi(outcar="OUTCAR"):
         if not os.path.isfile(outcar):
-            logger.warn("OUTCAR file not found. E-fermi set to 0.0eV")
+            logger.warning("OUTCAR file not found. E-fermi set to 0.0eV")
             return None
         txt = open(outcar).read()
         efermi = re.search(r'E-fermi :\s*([-+]?[0-9]+[.]?[0-9]*([eE][-+]?[0-9]+)?)', txt).groups()[0]

diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -14,6 +14,7 @@
 from argparse import ArgumentParser
 import os
 import logging
+import re
 import numpy as np
 import matplotlib.pyplot as plt
 from ase.calculators.vasp import VaspChargeDensity

diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -1,27 +1,45 @@
 #!/usr/bin/env python3
 '''
 A script to plot work function using LOCPOT.
-
 Requirement: python3, numpy, matplotlib, ase
-
 Author: @Ionizing
 Date: 22:46, Jan 11th, 2021.
+
+CHANGELOG:
+    - 1:56, Jan 25th, 2021
+        - Add E-fermi level correction
+        - More logging info
 '''
 
 from argparse import ArgumentParser
+import os
+import logging
 import numpy as np
 import matplotlib.pyplot as plt
-from ase.calculators.vasp.vasp import VaspChargeDensity
+from ase.calculators.vasp import VaspChargeDensity
 
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
 
-def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat'):
+
+def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat', outcar="OUTCAR"):
     '''
     Reads the LOCPOT file and calculate the average potential along `axis`.
      @in: See function argument.
     @out:
           - xvals: grid data along selected axis;
           - mean: averaged potential corresponding to `xvals`.
     '''
+    def get_efermi(outcar="OUTCAR"):
+        if not os.path.isfile(outcar):
+            logger.warn("OUTCAR file not found. E-fermi set to 0.0eV")
+            return None
+        txt = open(outcar).read()
+        efermi = re.search(r'E-fermi :\s*([-+]?[0-9]+[.]?[0-9]*([eE][-+]?[0-9]+)?)', txt).groups()[0]
+        logger.info("Found E-fermi = {}".format(efermi))
+        return float(efermi)
+
+    logger.info("Loading LOCPOT file {}".format(fname))
     locd = VaspChargeDensity(fname)
     cell = locd.atoms[0].cell
     latlens = np.linalg.norm(cell, axis=1)
@@ -33,17 +51,28 @@ def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat'):
     axes = tuple(axes)
 
     locpot = locd.chg[0]
-    mean = np.mean(locpot, axes) * vol  # must multiply with cell volume, similar to CHGCAR
+    # must multiply with cell volume, similar to CHGCAR
+    logger.info("Calculating workfunction along {} axis".format(axis))
+    mean = np.mean(locpot, axes) * vol
 
     xvals = np.linspace(0, latlens[iaxis], locpot.shape[iaxis])
 
     # save to 'locpot.dat'
-    np.savetxt(savefile, np.c_[xvals, mean], fmt='%13.5f', header='Distance(A) Potential(eV)')
+    efermi = get_efermi(outcar)
+    logger.info("Saving raw data to {}".format(savefile))
+    if efermi is None:
+        np.savetxt(savefile, np.c_[xvals, mean],
+                   fmt='%13.5f', header='Distance(A) Potential(eV) # E-fermi not corrected')
+    else:
+        mean -= efermi
+        np.savetxt(savefile, np.c_[xvals, mean],
+                   fmt='%13.5f', header='Distance(A) Potential(eV) # E-fermi shifted to 0.0eV')
     return (xvals, mean)
 
 
 def parse_cml_arguments():
-    parser = ArgumentParser(description='A tool to plot work function according to LOCPOT', add_help=True)
+    parser = ArgumentParser(
+        description='A tool to plot work function according to LOCPOT', add_help=True)
     parser.add_argument('-a', '--axis', type=str, action='store',
                         help='Which axis to be calculated: x, y or z. Default by z', default='z', choices=['x', 'y', 'z'])
     parser.add_argument('input', nargs='?', type=str,
@@ -57,13 +86,13 @@ def parse_cml_arguments():
     parser.add_argument('--title', type=str, action='store',
                         help='Title in output image. If none, no title is added, default is None', default=None)
     return parser.parse_args()
-    pass
 
 
 if '__main__' == __name__:
     args = parse_cml_arguments()
     x, y = locpot_mean(args.input, args.axis, args.write)
 
+    logger.info("Plotting to image")
     plt.plot(x, y, color='k')
     plt.xlabel('Distance(A)')
     plt.ylabel('Potential(eV)')
@@ -73,5 +102,6 @@ def parse_cml_arguments():
 
     if args.title:
         plt.title(args.title)
+
+    logger.info("Saving to {}".format(args.output))
     plt.savefig(args.output, dpi=args.dpi)
-    pass
diff --git a/plot-workfunc.py b/plot-workfunc.py
@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+'''
+A script to plot work function using LOCPOT.
+
+Requirement: python3, numpy, matplotlib, ase
+
+Author: @Ionizing
+Date: 22:46, Jan 11th, 2021.
+'''
+
+from argparse import ArgumentParser
+import numpy as np
+import matplotlib.pyplot as plt
+from ase.calculators.vasp.vasp import VaspChargeDensity
+
+
+def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat'):
+    '''
+    Reads the LOCPOT file and calculate the average potential along `axis`.
+     @in: See function argument.
+    @out:
+          - xvals: grid data along selected axis;
+          - mean: averaged potential corresponding to `xvals`.
+    '''
+    locd = VaspChargeDensity(fname)
+    cell = locd.atoms[0].cell
+    latlens = np.linalg.norm(cell, axis=1)
+    vol = np.linalg.det(cell)
+
+    iaxis = ['x', 'y', 'z'].index(axis.lower())
+    axes = [0, 1, 2]
+    axes.remove(iaxis)
+    axes = tuple(axes)
+
+    locpot = locd.chg[0]
+    mean = np.mean(locpot, axes) * vol  # must multiply with cell volume, similar to CHGCAR
+
+    xvals = np.linspace(0, latlens[iaxis], locpot.shape[iaxis])
+
+    # save to 'locpot.dat'
+    np.savetxt(savefile, np.c_[xvals, mean], fmt='%13.5f', header='Distance(A) Potential(eV)')
+    return (xvals, mean)
+
+
+def parse_cml_arguments():
+    parser = ArgumentParser(description='A tool to plot work function according to LOCPOT', add_help=True)
+    parser.add_argument('-a', '--axis', type=str, action='store',
+                        help='Which axis to be calculated: x, y or z. Default by z', default='z', choices=['x', 'y', 'z'])
+    parser.add_argument('input', nargs='?', type=str,
+                        help='The input file name, default by LOCPOT', default='LOCPOT')
+    parser.add_argument('-w', '--write', type=str, action='store',
+                        help='Save raw work function data to file, default by locpot.dat', default='locpot.dat')
+    parser.add_argument('-o', '--output', type=str, action='store',
+                        help='Output image file name, default by Workfunction.png', default='Workfunction.png')
+    parser.add_argument('--dpi', type=int, action='store',
+                        help='DPI of output image, default by 400', default=400)
+    parser.add_argument('--title', type=str, action='store',
+                        help='Title in output image. If none, no title is added, default is None', default=None)
+    return parser.parse_args()
+    pass
+
+
+if '__main__' == __name__:
+    args = parse_cml_arguments()
+    x, y = locpot_mean(args.input, args.axis, args.write)
+
+    plt.plot(x, y, color='k')
+    plt.xlabel('Distance(A)')
+    plt.ylabel('Potential(eV)')
+    plt.grid(color='gray', ls='-.')
+    plt.xlim(0, np.max(x))
+    plt.ylim(np.min(y)-0.5, np.max(y)+0.5)
+
+    if args.title:
+        plt.title(args.title)
+    plt.savefig(args.output, dpi=args.dpi)
+    pass