diff --git a/tools/tabulate/tabulate.py b/tools/tabulate/tabulate.py index 3bba7c2ecf..c3cb3a256f 100755 --- a/tools/tabulate/tabulate.py +++ b/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") diff --git a/tools/tabulate/wall_harmonic_tabulate.py b/tools/tabulate/wall_harmonic_tabulate.py new file mode 100755 index 0000000000..9424856dbf --- /dev/null +++ b/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') diff --git a/tools/tabulate/wall_multi_tabulate.py b/tools/tabulate/wall_multi_tabulate.py new file mode 100755 index 0000000000..0458a929e2 --- /dev/null +++ b/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')