add support to create tables for fix wall/table

tools/tabulate/tabulate.py

@@ -279,6 +279,53 @@ class DihedralTabulate(Tabulate):
         if self.args.filename != '-':
             self.fp.close()
 
+################################################################################
+# create tabulation for fix wall/table
+class WallTabulate(Tabulate):
+    def __init__(self, efunc, ffunc=None, units=None, comment=None):
+        super(WallTabulate, self).__init__('wall', 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")
+        idx = [a.dest for a in self.parser._actions].index('xmin')
+        self.parser._actions[idx].required=False
+        self.parser._actions[idx].default=0.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')
+
+        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(self.args.xmax)
+
+        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\n\n" % (self.args.num, diffmin, diffmax))
+        else:
+            self.fp.write("N %d\n\n" % (self.args.num))
+
+        self.writetable(table, offset)
+        if self.args.filename != '-':
+            self.fp.close()
+
 ################################################################################
 if __name__ == "__main__":
     sys.exit("The tabulate module is not meant to be executed directly")

tools/tabulate/wall_harmonic_tabulate.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+from tabulate import WallTabulate
+
+################################################################################
+import math
+
+k = 100.0
+rzero = 4.0
+
+def harmonic_force(r):
+    dr = r - rzero
+    f = -2.0 * k * dr
+    return f
+
+def harmonic_energy(r):
+    dr = r - rzero
+    f = k * dr * dr
+    return f
+
+################################################################################
+
+if __name__ == "__main__":
+    wtable = WallTabulate(harmonic_energy, harmonic_force, units='real')
+    wtable.run('HARMONIC')

tools/tabulate/wall_multi_tabulate.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+
+from tabulate import WallTabulate
+import os, sys
+################################################################################
+import math
+
+k = 100.0
+rzero = 4.0
+epsilon = 0.02
+sigma = 2.0
+depth = 20.0
+width = 2.0
+rzero = 1.2
+
+def harmonic_force(r):
+    dr = r - rzero
+    f = -2.0 * k * dr
+    return f
+
+def harmonic_energy(r):
+    dr = r - rzero
+    f = k * dr * dr
+    return f
+
+def lj126_force(r):
+    f = -4.0*epsilon*(-12.0*math.pow(sigma/r,12.0)/r + 6.0*math.pow(sigma/r,6.0)/r)
+    return f
+
+def lj126_energy(r):
+    f = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0))
+    return f
+
+def morse_energy(r):
+    ralpha = math.exp(-width*(r-rzero))
+    f = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha))
+    return f
+
+def morse_force(r):
+    ralpha = math.exp(-width*(r-rzero))
+    f = -2.0 * depth * width * (1.0 -ralpha) * ralpha
+    return f
+
+
+################################################################################
+
+if __name__ == "__main__":
+    fname = 'walltab.dat'
+    if os.path.exists(fname):
+        os.remove(fname)
+    sys.argv.append('--filename')
+    sys.argv.append(fname)
+    sys.argv.append('--num-points')
+    sys.argv.append('1000')
+    sys.argv.append('--inner')
+    sys.argv.append('0.02')
+    sys.argv.append('--outer')
+    sys.argv.append('4.0')
+    wtable = WallTabulate(harmonic_energy, harmonic_force, units='real')
+    wtable.run('HARMONIC')
+    wtable = WallTabulate(lj126_energy, lj126_force, units='real')
+    wtable.run('LJ126')
+    wtable = WallTabulate(morse_energy, morse_force, units='real')
+    wtable.run('MORSE')