import tabulate scripts for table files

tools/tabulate/README.md

@@ -0,0 +1,39 @@
+# Python scripts to generate tabulated potential files for LAMMPS
+
+This directory contains a Python module 'tabulate' that can be used to
+generate tabulated potential files for pair style table, bond style
+table, and angle style table
+
+To create tables, you need to define your energy and - optionally -
+force functions and then an instance of either the PairTabulate(),
+BondTabulate(), AngleTabulate(), or DihedralTabulate() class from the
+tabulate module and call its run() method to generate the table.  Most
+of the settings (number of points, minimum, maximum etc.)  are provided
+as command line flags.  The run() method may be called multiple times to
+generate multiple tables, for instance after changing parameters of the
+energy/force functions.  If the force function is not provided, the
+force will be determined through numerical differentiation.
+
+Please see the individual tabulation scripts in this folder for examples:
+
+| ------------------------------|----------------------------------------------------------------------------|
+| File                          | Description                                                                |
+| ------------------------------|----------------------------------------------------------------------------|
+| pair_lj_tabulate.py           | creates two Lennard-Jones pair potential tables with different parameters  |
+| bond_morse_tabulate.py        | creates a table for a Morse bond potential table                           |
+| angle_harmonic_tabulate.py    | creates a table for a harmonic angle potential table                       |
+| dihedral_harmonic_tabulate.py | creates a table for a harmonic dihedral potential table                    |
+| pair_hybrid_tabulate.py       | creates a Morse/Lennard-Jones hybrid potential table with smooth switching |
+| ------------------------------|----------------------------------------------------------------------------|
+
+Common command line flags:
+
+```
+options:
+  -h, --help                          show this help message and exit
+  --num-points NUM, -n NUM            Number of tabulated points (default: 1000)
+  --filename FILENAME, -f FILENAME    Name of output file (default: -)
+  --diff-num, -d                      Differentiate energy function numerically
+  --inner XMIN, -i XMIN               Inner cutoff of table (required for pair)
+  --outer XMAX, -o XMAX               Outer cutoff of table (required)
+```

tools/tabulate/angle_harmonic_tabulate.py

@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+
+from tabulate import AngleTabulate
+
+################################################################################
+import math
+
+def harmonic_energy(theta):
+    k = 100.0
+    thetazero = 120.0
+    # the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
+    deg2rad = math.pi / 180.0
+    t = (theta - thetazero) * deg2rad
+    f = k * t * t
+    return f
+
+################################################################################
+
+if __name__ == "__main__":
+    atable = AngleTabulate(harmonic_energy, units='real')
+    atable.run('HARM')

tools/tabulate/bond_morse_tabulate.py

@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+
+from tabulate import BondTabulate
+
+################################################################################
+import math
+
+def morse_energy(r):
+    depth = 1.0
+    width = 2.0
+    rzero = 1.2
+    ralpha = math.exp(-width*(r-rzero))
+    f = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha))
+    return f
+
+def morse_force(r):
+    depth = 1.0
+    width = 2.0
+    rzero = 1.2
+    ralpha = math.exp(-width*(r-rzero))
+    f = -2.0 * depth * width * (1.0 -ralpha) * ralpha
+    return f
+
+################################################################################
+
+if __name__ == "__main__":
+    btable = BondTabulate(morse_energy, morse_force, units='lj')
+    btable.run('MORSE')

tools/tabulate/dihedral_harmonic_tabulate.py

@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+
+from tabulate import DihedralTabulate
+
+################################################################################
+import math
+
+def harmonic_energy(theta):
+    k = 100.0
+    # the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
+    deg2rad = math.pi / 180.0
+    f = k * (1 - math.cos(2.0 * deg2rad * theta))
+    return f
+
+################################################################################
+
+if __name__ == "__main__":
+    dtable = DihedralTabulate(harmonic_energy, units='metal')
+    dtable.run('HARM')

tools/tabulate/pair_hybrid_tabulate.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+
+from tabulate import PairTabulate
+
+################################################################################
+import math
+
+# hybrid potential using Morse for repulsive and LJ for attractive.
+# using tanh() as smooth switching function between the two.
+
+def hybrid_energy(r):
+    depth = 1.0
+    width = 5.0
+    epsilon = 1.0
+    sigma = 1.0
+    rzero = 1.2
+
+    f1 = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0))
+
+    ralpha = math.exp(-width*(r-rzero))
+    f2 = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha))
+
+    switch = 0.5*math.tanh(10.0*(r - rzero)) + 0.5
+    f = switch*f1 + (1.0 - switch)*f2
+    return f
+
+################################################################################
+
+if __name__ == "__main__":
+    ptable = PairTabulate(hybrid_energy, units='lj', comment='Morse repulsion + LJ attraction')
+    ptable.run('MORSE_LJ')

tools/tabulate/pair_lj_tabulate.py

@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+from tabulate import PairTabulate
+
+################################################################################
+import math
+epsilon = 1.0
+sigma = 1.0
+
+def lj_energy(r):
+    f = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0))
+    return f
+
+def lj_force(r):
+    epsilon = 1.0
+    sigma = 1.0
+    f = -4.0*epsilon*(-12.0*math.pow(sigma/r,12.0)/r + 6.0*math.pow(sigma/r,6.0)/r)
+    return f
+################################################################################
+
+if __name__ == "__main__":
+    ptable = PairTabulate(lj_energy, lj_force)
+    ptable.run('LJ_11')
+    epsilon = 1.0
+    sigma = 1.5
+    ptable.run('LJ_12')

tools/tabulate/tabulate.py

@@ -0,0 +1,279 @@
+# ----------------------------------------------------------------------
+#   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+#   https://www.lammps.org/ Sandia National Laboratories
+#   LAMMPS Development team: developers@lammps.org
+#
+#   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+#   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+#   certain rights in this software.  This software is distributed under
+#   the GNU General Public License.
+#
+#   See the README file in the top-level LAMMPS directory.
+# -------------------------------------------------------------------------
+"""
+Python classes and support functions to generate tabulated potential files for LAMMPS.
+"""
+
+# for python2/3 compatibility
+from __future__ import print_function
+
+import sys
+import argparse
+from os import path
+from datetime import datetime
+
+########################################################################
+
+def numdiff(x, func):
+    """ Get the value of the first derivative of a function 'func(x)' at 'x'
+        from numerical differentiation."""
+
+    # optimal delta x value for 2-point numerical differentiation of two floating point numbers
+    epsilon = 6.05504e-6
+    fval1 = func(x - epsilon)
+    fval2 = func(x + epsilon)
+    return 0.5 * (fval2-fval1) / epsilon
+
+########################################################################
+
+def mktable(tstyle, label, num, xmin, xmax, efunc, diff=False, ffunc=None):
+    """ Do the tabulation of the provided energy function. Compute force from
+    numerical differentiation if no force function is provided.  Also detect
+    minimum for use to determine potential shifting in bonded potentials."""
+
+    # must use numerical differentiation if no force function provided
+    if not ffunc:
+        diff = True
+
+    print("# Creating %s table %s with %d points from %g to %g" % (tstyle, label, num, xmin, xmax))
+
+    table = []
+    delx = (xmax - xmin) / (float(num) - 1.0)
+    emin = 999999999.0
+    xzero = 0.0
+
+    for i in range(0,num):
+        x = xmin + i*delx
+        energy = efunc(x)
+        if energy < emin:
+            emin = energy
+            xzero = x
+        if diff:
+            force = -numdiff(x, efunc)
+        else:
+            force = ffunc(x)
+        table.append([i+1, x, energy, force])
+
+    return table, xzero
+
+
+########################################################################
+# base class with shared functionality
+class Tabulate(object):
+    """Base tabulation class. Contains all shared functionality: common argument parsing,
+    output file handling, table writing"""
+
+    def __init__(self, style, efunc, ffunc=None, units=None, comment=None):
+        self.fp = sys.stdout
+        self.tstyle = style
+        self.energyfunc = efunc
+        self.forcefunc = ffunc
+        self.units = units
+        self.comment = comment
+        self.eshift = 0.0
+        self.args = None
+        self.diff = True
+
+        self.parser = argparse.ArgumentParser(description='Tool to generate tabulated '
+                                              + self.tstyle + ' potential files for LAMMPS')
+        self.parser.add_argument('--num-points', '-n', dest='num', default=1000, type=int,
+                                 help="Number of tabulated points")
+        self.parser.add_argument('--filename', '-f', dest='filename', default='-',
+                                 help="Name of output file")
+        self.parser.add_argument('--diff-num', '-d', dest='diff',default=False,
+                                 action='store_true',
+                                 help="Differentiate energy function numerically")
+        self.parser.add_argument('--inner', '-i', dest='xmin', required=True, type=float,
+                                 help="Inner cutoff of table")
+        self.parser.add_argument('--outer', '-o', dest='xmax', required=True, type=float,
+                                 help="Outer cutoff of table")
+
+    def openfile(self, label):
+        """Open table file, if needed and print label for new table entry"""
+        if self.args and self.args.filename != '-':
+            try:
+                if path.isfile(self.args.filename):
+                    self.fp = open(self.args.filename, 'a')
+                    print("# Appending table to file " + self.args.filename)
+                else:
+                    self.fp = open(self.args.filename, 'w')
+                    print("# Writing table to new file " + self.args.filename)
+                    self.fp.write('# DATE: ' + datetime.now().date().isoformat())
+                    if self.units:
+                        self.fp.write(' UNITS: ' + str(self.units))
+                    if self.comment:
+                        self.fp.write(' COMMENT: ' + str(self.comment) + '\n')
+            except IOError:
+                sys.exit("Cannot open file %s for writing table data" % self.args.filename)
+        self.fp.write('\n' + label + '\n')
+
+    def writetable(self, table, offset):
+        """ Formatted output tabulated data with 4 columns"""
+        for i,r,energy,force in table:
+            self.fp.write("%8d %- 22.15g %- 22.15g %- 22.15g\n" % (i, r, energy - offset, force))
+
+    def helpexit(self, text):
+        """ Convenience function to exit program with error and help message"""
+        sys.exit('\n' + text + '\n\n' + self.parser.format_help())
+
+################################################################################
+# create tabulation for pair styles
+class PairTabulate(Tabulate):
+    def __init__(self, efunc, ffunc=None, units=None, comment=None):
+        super(PairTabulate, self).__init__('pair', efunc, ffunc, units, comment)
+        self.parser.add_argument('--eshift', '-e', dest='eshift', default=False,
+                                 action='store_true',
+                                 help="Shift potential energy to be zero at outer cutoff")
+        try:
+            self.args = self.parser.parse_args()
+        except argparse.ArgumentError:
+            sys.exit()
+
+    def run(self, label):
+        # sanity checks
+        if self.args.num < 2:
+            self.helpexit('Expect 2 or more points in table for tabulation')
+        if self.args.xmin <= 0.0:
+            self.helpexit('Inner tabulation cutoff must be > 0 for pair style table')
+        if self.args.xmax <= self.args.xmin:
+            self.helpexit('Outer cutoff must be larger than inner cutoff')
+
+        self.diff = self.args.diff
+        if not self.forcefunc:
+            self.diff = True
+        offset = 0.0
+        if self.args.eshift:
+            offset=self.energyfunc(self.args.xmax)
+
+        table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
+                               self.energyfunc, self.args.diff, self.forcefunc)
+
+        # open table file
+        self.openfile(label)
+
+        # write pair style specific header
+        if self.forcefunc:
+            diffmin = -numdiff(self.args.xmin, self.forcefunc)
+            diffmax = -numdiff(self.args.xmax, self.forcefunc)
+            self.fp.write("N %d R %g %g FPRIME %- 22.15g %- 22.15g\n\n"
+                          % (self.args.num, self.args.xmin, self.args.xmax, diffmin, diffmax))
+        else:
+            self.fp.write("N %d R %g %g\n\n" % (self.args.num, self.args.xmin, self.args.xmax))
+
+        self.writetable(table, offset)
+        if self.args.filename != '-':
+            self.fp.close()
+
+
+################################################################################
+# shared functionality to create tabulation for bond or angle styles
+class BondAngleTabulate(Tabulate):
+    def __init__(self, style, efunc, ffunc=None, units=None, comment=None):
+        super(BondAngleTabulate, self).__init__(style, efunc, ffunc, units, comment)
+        self.parser.add_argument('--eshift', '-e', dest='eshift', default=False,
+                                 action='store_true',
+                                 help="Shift potential energy to be zero at minimum")
+        idx = [a.dest for a in self.parser._actions].index('xmin')
+        self.parser._actions[idx].required=False
+        self.parser._actions[idx].default=0.0
+        if style == 'angle':
+            idx = [a.dest for a in self.parser._actions].index('xmax')
+            self.parser._actions[idx].required=False
+            self.parser._actions[idx].default=180.0
+        try:
+            self.args = self.parser.parse_args()
+        except argparse.ArgumentError:
+            sys.exit()
+
+    def run(self, label):
+        # sanity checks
+        if self.args.num < 2:
+            self.helpexit('Expect 2 or more points in table for tabulation')
+        if self.args.xmin < 0.0:
+            self.helpexit('Inner cutoff must not be negative')
+        if self.tstyle == 'angle' and self.args.xmax > 180.0:
+            self.helpexit('Outer cutoff must not be larger than 180.0 degrees')
+
+        self.diff = self.args.diff
+        if not self.forcefunc:
+            self.diff = True
+
+        table, xzero = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
+                               self.energyfunc, self.args.diff, self.forcefunc)
+        print("# Minimum energy of tabulated potential is at %g" % xzero)
+        offset = 0.0
+        if self.args.eshift:
+            offset=self.energyfunc(xzero)
+
+        self.openfile(label)
+
+        if self.forcefunc:
+            diffmin = -numdiff(self.args.xmin, self.forcefunc)
+            diffmax = -numdiff(self.args.xmax, self.forcefunc)
+            self.fp.write("N %d FP %- 22.15g %- 22.15g EQ %g\n\n" %
+                          (self.args.num, diffmin, diffmax, xzero))
+        else:
+            self.fp.write("N %d EQ %g\n\n" % (self.args.num, xzero))
+
+        self.writetable(table, offset)
+        if self.args.filename != '-':
+            self.fp.close()
+
+################################################################################
+class BondTabulate(BondAngleTabulate):
+    def __init__(self, efunc, ffunc=None, units=None, comment=None):
+        super(BondTabulate, self).__init__('bond', efunc, ffunc, units, comment)
+
+################################################################################
+class AngleTabulate(BondAngleTabulate):
+    def __init__(self, efunc, ffunc=None, units=None, comment=None):
+        super(AngleTabulate, self).__init__('angle', efunc, ffunc, units, comment)
+
+################################################################################
+# create tabulation for dihdedral
+class DihedralTabulate(Tabulate):
+    def __init__(self, efunc, ffunc=None, units=None, comment=None):
+        super(DihedralTabulate, self).__init__('dihedral', efunc, ffunc, units, comment)
+        idx = [a.dest for a in self.parser._actions].index('xmin')
+        self.parser._actions[idx].required=False
+        self.parser._actions[idx].default=-180.0
+        idx = [a.dest for a in self.parser._actions].index('xmax')
+        self.parser._actions[idx].required=False
+        self.parser._actions[idx].default=179.999999
+        try:
+            self.args = self.parser.parse_args()
+        except argparse.ArgumentError:
+            sys.exit()
+
+    def run(self, label):
+        # sanity checks
+        if self.args.num < 2:
+            self.helpexit('Expect 2 or more points in table for tabulation')
+        if self.args.xmin < -180 or self.args.xmin > 0.0:
+            self.helpexit('Inner cutoff must be within -180.0 and 0.0 degrees')
+        if self.args.xmax < 0.0 or self.args.xmin > 360.0:
+            self.helpexit('Outer cutoff must be within 0.0 and 360.0 degrees')
+        if (self.args.xmax - self.args.xmin) >= 360.0:
+            self.helpexit('Inner and outer cutoff range must be less than 360.0 degrees')
+
+        self.diff = self.args.diff
+        if not self.forcefunc:
+            self.diff = True
+
+        table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
+                               self.energyfunc, self.args.diff, self.forcefunc)
+        self.openfile(label)
+        self.fp.write("N %d DEGREES \n\n" % (self.args.num))
+        self.writetable(table, 0.0)
+        if self.args.filename != '-':
+            self.fp.close()