diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst
index e502480ce5..488de848bf 100644
--- a/doc/src/Tools.rst
+++ b/doc/src/Tools.rst
@@ -1250,10 +1250,10 @@ tabulate tool
.. versionadded:: 22Dec2022
-The ``tabulate`` folder contains Python scripts scripts to generate tabulated
-potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module
-in the ``tabulate.py`` file. Some example files demonstrating its use are
-included. See the README file for more information.
+The ``tabulate`` folder contains Python scripts scripts to generate and
+visualize tabulated potential files for LAMMPS. The bulk of the code is in the
+``tabulate`` module in the ``tabulate.py`` file. Some example files
+demonstrating its use are included. See the README file for more information.
----------
@@ -1276,7 +1276,7 @@ Those scripts were written by Steve Plimpton sjplimp at gmail.com
valgrind tool
-------------
-The ``valgrind`` folder contains additional suppressions fur LAMMPS when
+The ``valgrind`` folder contains additional suppressions for LAMMPS when
using `valgrind's `_ ` `memcheck tool
`_ to search for memory
access violation and memory leaks. These suppressions are automatically
diff --git a/tools/tabulate/README.md b/tools/tabulate/README.md
index d3b4df0738..995bb8abdd 100644
--- a/tools/tabulate/README.md
+++ b/tools/tabulate/README.md
@@ -27,6 +27,7 @@ Please see the individual tabulation scripts in this folder for examples:
| wall_harmonic_tabulate.py | creates a table for fix wall/table with a simple repulsive harmonic potential |
| wall_multi_tabulate.py | creates a table for fix wall/table with multiple tables |
| pair_bi_tabulate.py | creates a table from radial distribution file using Boltzmann Inversion |
+| plot_forces.py | plots and extracts tabulated forces from table files |
Common command line flags:
diff --git a/tools/tabulate/plot_forces.py b/tools/tabulate/plot_forces.py
new file mode 100755
index 0000000000..ee55a93be6
--- /dev/null
+++ b/tools/tabulate/plot_forces.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Author: Germain Clavier (Unicaen), germain.clavier at unicaen.fr
+
+"""
+Plot LAMMPS tabulated forces.
+"""
+
+import argparse
+import numpy as np
+import os
+import logging
+import sys
+from matplotlib import pyplot as plt
+
+logger = logging.getLogger(__name__)
+
+
+units = {
+ "lj": {
+ "Distance": "Reduced units",
+ "Energy": "Reduced units",
+ "Force": "Reduced units",
+ "kb": 1,
+ },
+ "real": {
+ "Distance": "[A]",
+ "Energy": "[kcal/mol]",
+ "Force": "[kcal/mol/A]",
+ "kb": 0.001985875,
+ },
+ "metal": {
+ "Distance": "[A]",
+ "Energy": "[eV]",
+ "Force": "[eV/A]",
+ "kb": 8.6173332e-5,
+ },
+ "si": {"Distance": "[m]", "Energy": "[J]", "Force": "[N]", "kb": 1.380649e-23},
+ "cgs": {
+ "Distance": "[cm]",
+ "Energy": "[ergs]",
+ "Force": "[dynes]",
+ "kb": 1.3806504e-16,
+ },
+ "electron": {
+ "Distance": "[Bohr]",
+ "Energy": "[Hartrees]",
+ "Force": "[Hartree/Bohr]",
+ "kb": 3.16681534e-6,
+ },
+ "micro": {
+ "Distance": "[um]",
+ "Energy": "[pg·um^2/us^2]",
+ "Force": "[pg·um/us^2]",
+ "kb": 1.3806504e-8,
+ },
+ "nano": {
+ "Distance": "[nm]",
+ "Energy": "[ag·nm^2/ns^2]",
+ "Force": "[ag·nm/ns^2]",
+ "kb": 0.013806504,
+ },
+}
+
+
+def compute_energy(tp):
+ r = tp[0]
+ fo = tp[2]
+ e = np.zeros(r.shape)
+ for i, (ri, fi) in enumerate(zip(r, fo)):
+ if i == 0:
+ continue
+ dr = ri - r[i - 1]
+ e[i] = e[i - 1] - dr * fo[i - 1]
+ e -= e[-1]
+ return e
+
+
+def main():
+
+ parser = argparse.ArgumentParser(
+ description="""
+ Plots LAMMPS tabulated forces. This script takes a table
+ file as an input and plots all the tabulated forces inside with their
+ corresponding energy. The forces label is the token used to name the
+ force in the file. It can be used to output all the forces in separate
+ files and/or recompute the energy from forces through finite difference
+ (assuming e(rc)=0). This script requires the matplotlib and numpy
+ Python libraries. Bitmap format is not supported.
+ """
+ )
+ parser.add_argument(
+ "-u",
+ "--units",
+ dest="units",
+ default="real",
+ help="Units of the file (LAMMPS units system)",
+ )
+ parser.add_argument(
+ "-f",
+ "--file",
+ dest="infile",
+ default="",
+ help="File to read",
+ )
+ parser.add_argument(
+ "-x",
+ dest="xrange",
+ default="",
+ help="xrange separated by : (for negative values use the '=' sign: -x=-3:10)",
+ )
+ parser.add_argument(
+ "-y",
+ dest="yrange",
+ default="",
+ help="yrange separated by :",
+ )
+ parser.add_argument(
+ "-t",
+ dest="temp",
+ default=None,
+ type=float,
+ help="temperature for KbT plot [default none]",
+ )
+ parser.add_argument(
+ "-d",
+ "--diff-num",
+ dest="recompute",
+ action="store_true",
+ help="Recompute the energies from forces and distances through finite differences",
+ )
+ parser.add_argument(
+ "-e",
+ dest="extract",
+ action="store_true",
+ help="Extract the forces in separate files",
+ )
+ args = parser.parse_args()
+ logging.basicConfig(level=logging.INFO)
+
+ ##########
+ # Manage arguments
+
+ udistance = units[args.units]["Distance"]
+ uenergy = units[args.units]["Energy"]
+ uforce = units[args.units]["Force"]
+ kb = units[args.units]["kb"]
+ rlabel = " ".join(["Rij", udistance])
+ elabel = " ".join(["E", uenergy])
+ flabel = " ".join(["F", uforce])
+ etitle = "Energy"
+ ftitle = "Force"
+ font = "DejaVu Sans"
+ fontsize = 30
+
+ infile = args.infile
+ if not os.path.isfile(infile):
+ logger.error("Input file not found")
+ sys.exit(1)
+
+ toplot = []
+ with open(infile, "r") as f:
+ lines = iter(f.readlines())
+ while True:
+ try:
+ r = []
+ force = []
+ ener = []
+ tok = []
+ while not tok:
+ tok = next(lines).partition("#")[0].rstrip()
+ logger.info("Found {} token".format(tok))
+ infos = next(lines).split()
+ npoints = int(infos[1])
+ next(lines)
+ if "bitmap" in infos:
+ logger.info("Unsupported bitmap format for token {:s}".format(tok))
+ for _ in range(npoints):
+ continue
+ else:
+ for i in range(npoints):
+ line = next(lines).split()
+ r.append(float(line[1]))
+ ener.append(float(line[2]))
+ force.append(float(line[3]))
+ r = np.array(r)
+ ener = np.array(ener)
+ force = np.array(force)
+ toplot.append([r, ener, force, tok])
+ tok = []
+ next(lines)
+ except StopIteration:
+ break
+ if args.recompute:
+ etitle = "Estimated energy"
+ for tp in toplot:
+ tp[1] = compute_energy(tp)
+
+ fig, axes = plt.subplots(1, 2)
+
+ for tp in toplot:
+ axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3)
+ axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3)
+ 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",
+ )
+ axes[0].text(hmax / 2.0, kb * args.temp, "KbT", fontsize=0.7 * fontsize)
+ logger.info("KbT value= {:e} {:s}".format(kb * args.temp, uenergy))
+ else:
+ logger.info("Invalid temperature value: {:e}".format(args.temp))
+
+ if args.xrange:
+ xmin, xmax = list(map(float, args.xrange.split(":")))
+ axes[0].set_xlim(xmin, xmax)
+ axes[1].set_xlim(xmin, xmax)
+ if args.yrange:
+ ymin, ymax = list(map(float, args.yrange.split(":")))
+ axes[0].set_ylim(ymin, ymax)
+ axes[1].set_ylim(ymin, ymax)
+
+ # Setting axes 0
+ axes[0].set_title(etitle, fontsize=fontsize)
+ axes[0].set_xlabel(
+ rlabel, fontname=font, fontsize=fontsize
+ ) # xlabel name, size 30pts
+ axes[0].set_ylabel(
+ elabel, fontname=font, fontsize=fontsize
+ ) # ylabel name, size 30pts
+ axes[0].tick_params(
+ axis="both", which="major", labelsize=fontsize
+ ) # Biggers ticks, bigger tick labels!
+
+ # Setting axes 1
+ axes[1].set_title(ftitle, fontsize=fontsize)
+ axes[1].legend(frameon=False, fontsize=fontsize) # Fat font, no frame
+ axes[1].set_xlabel(
+ rlabel, fontname=font, fontsize=fontsize
+ )
+ axes[1].set_ylabel(
+ flabel, fontname=font, fontsize=fontsize
+ )
+ axes[1].tick_params(
+ axis="both", which="major", labelsize=fontsize
+ )
+
+ figManager = plt.get_current_fig_manager()
+ figManager.window.showMaximized()
+ plt.show()
+
+ if args.extract:
+ for tp in toplot:
+ outfile = "".join([tp[3], ".plot"])
+ logger.info("Writing file {}".format(outfile))
+ with open(outfile, "w") as f:
+ f.write("# {} force extracted from {}\n".format(tp[3], infile))
+ f.write("# {:^20} {:^20} {:^20}\n".format("r", "energy", "force"))
+ for a, b, c in zip(tp[0], tp[1], tp[2]):
+ f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c))
+ return
+
+
+if __name__ == "__main__":
+ try:
+ main()
+ except KeyboardInterrupt:
+ raise SystemExit("User interruption.")