Rewritten plot_forces for more generic case

lammps/plot_forces.py

@@ -25,6 +25,16 @@ def compute_energy(tp):
     e -= e[-1]
     return e
 
+units = {
+    "lj": {"Energy": "Reduced units", "Force": "Reduced units", "kb": 1},
+    "real": {"Energy": "kcal/mol", "Force": "kcal/mol/A", "kb": 0.001985875},
+    "metal": {"Energy": "eV", "Force": "eV/A", "kb": 8.6173332e-5},
+    "si": {"Energy": "J", "Force": "N", "kb": 1.380649e-23},
+    "cgs": {"Energy": "ergs", "Force": "dynes", "kb": 1.3806504e-16},
+    "electron": {"Energy": "Hartrees", "Force": "Hartree/Bohr", "kb": 3.16681534e-6},
+    "micro": {"Energy": "pg·um^2/us^2", "Force": "pg·um/us^2", "kb": 1.3806504e-8},
+    "nano": {"Energy": "ag·nm^2/ns^2", "Force": "ag·nm/ns^2", "kb": 0.013806504},
+    }
 
 def main():
 
@@ -32,12 +42,11 @@ def main():
         description="Plots Lammps tabulated forces."
     )
     parser.add_argument(
-        "-v",
-        "--verbose",
-        dest="verbosity",
-        default=0,
-        action="count",
-        help="display more information",
+        "-u",
+        "--units"
+        dest="units"
+        default="real"
+        help="Units of the file (LAMMPS units system)",
     )
     parser.add_argument(
         "-f",
@@ -50,19 +59,25 @@ def main():
         "-x",
         dest="xrange",
         default="",
-        help="xrange : separated (for negative values: -x=-3:10)",
+        help="xrange separated by : (for negative values: -x=-3:10)",
     )
     parser.add_argument(
         "-y",
         dest="yrange",
         default="",
-        help="yrange : separated",
+        help="yrange separated by :",
     )
     parser.add_argument(
         "-t",
         dest="temp",
-        default=-1, type=float,
-        help="temperature for kbT plot [default none]",
+        default=None, type=float,
+        help="temperature for KbT plot [default none]",
+    )
+    parser.add_argument(
+        "--recompute-energy",
+        dest="recompute",
+        action="store_true",
+        help="Recompute the energies from forces and distances through finite differences",
     )
     parser.add_argument(
         "-e",
@@ -75,15 +90,9 @@ def main():
     ##########
     # Manage arguments
     
-    # -v/--verbose
-
-    kb = 0.001985875  # kcal/K/mol
-
-    try:
-        utilsscript.init_logging(args.verbosity)
-    except NameError:
-        logging.basicConfig(level=logging.WARNING)
-        logging.warning("utilsscript lib not found, using default logging at warning level.")
+    uenergy = units(args.units)["Energy"]
+    uforce = units(args.units)["Force"]
+    kb = units(args.units)["kb"]
 
     infile = args.infile
     if not os.path.isfile(infile):
@@ -118,6 +127,7 @@ def main():
                 next(lines)
             except StopIteration:
                 break
+    if args.recompute:
         for tp in toplot:
             tp[1] = compute_energy(tp)
 
@@ -129,9 +139,13 @@ def main():
         hmin, hmax = axes[1].get_xlim()
         axes[1].hlines(0, hmin, hmax, color="black", linewidth=3, linestyles="dashdot")
 
+    if args.temp:
         if args.temp > 0:
             hmin, hmax = axes[0].get_xlim()
             axes[0].hlines(kb*args.temp, hmin, hmax, color="orange", label=r'$k_BT$', linewidth=3, linestyles="dashdot")
+            logging.info("KbT value= {:e}".format(kb*args.temp))
+        else:
+            logging.info("Invalid temperature value: {:e}".format(args.temp))
 
     if args.xrange:
         xmin, xmax = list(map(float, args.xrange.split(":")))
@@ -142,7 +156,7 @@ def main():
         axes[0].set_ylim(ymin, ymax)
         axes[1].set_ylim(ymin, ymax)
 
-    font = "Hack"
+    font = "DejaVu Sans"
     fontsize = 30
 
     # Setting axes 0