ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'lammps:develop' into fortran-fix-external

Browse Source
This commit is contained in:
hammondkd
2022-12-16 12:15:27 -06:00
committed by GitHub
parent dc0496ed48 a633915829
commit 14fd40acc5
13 changed files with 485 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
doc/src/Tools.rst
View File

@ -57,7 +57,7 @@ Pre-processing tools
* :ref:`msi2lmp <msi>`
* :ref:`polybond <polybond>`
* :ref:`stl_bin2txt <stlconvert>`
* :ref:`tabulate <tabulate>`
Post-processing tools
=====================
@ -1159,13 +1159,27 @@ For illustration purposes below is a part of the Tcl example script.
----------
.. _tabulate:
tabulate tool
--------------
.. versionadded:: TBD
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.
----------
.. _vim:
vim tool
------------------
The files in the tools/vim directory are add-ons to the VIM editor
that allow easier editing of LAMMPS input scripts. See the README.txt
The files in the ``tools/vim`` directory are add-ons to the VIM editor
that allow easier editing of LAMMPS input scripts. See the ``README.txt``
file for details.
These files were provided by Gerolf Ziegenhain (gerolf at

4
doc/src/angle_table.rst
View File

@ -59,6 +59,10 @@ format of this file is described below.
----------
Suitable tables for use with this angle style can be created using the
Python code in the ``tools/tabulate`` folder of the LAMMPS source code
distribution.
The format of a tabulated file is as follows (without the
parenthesized comments):

10
doc/src/bond_table.rst
View File

@ -59,6 +59,13 @@ this file is described below.
----------
Suitable tables for use with this bond style can be created by LAMMPS
itself from existing bond styles using the :doc:`bond_write
<bond_write>` command. This can be useful to have a template file for
testing the bond style settings and to build a compatible custom file.
Another option to generate tables is the Python code in the
``tools/tabulate`` folder of the LAMMPS source code distribution.
The format of a tabulated file is as follows (without the
parenthesized comments):
@ -149,7 +156,8 @@ info.
Related commands
""""""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`
:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`,
:doc:`bond_write <bond_write>`
Default
"""""""

4
doc/src/dihedral_table.rst
View File

@ -114,6 +114,10 @@ below.
----------
Suitable tables for use with this dihedral style can be created using
the Python code in the ``tools/tabulate`` folder of the LAMMPS source
code distribution.
The format of a tabulated file is as follows (without the
parenthesized comments). It can begin with one or more comment
or blank lines.

7
doc/src/pair_table.rst
View File

@ -120,6 +120,13 @@ best effect:
----------
Suitable tables in the correct format for use with these pair styles can
be created by LAMMPS itself using the :doc:`pair_write <pair_write>`
command. In combination with the :doc:`pair style python <pair_python>`
this can be a powerful mechanism to implement and test tables for use
with LAMMPS. Another option to generate tables is the Python code in
the ``tools/tabulate`` folder of the LAMMPS source code distribution.
The format of a tabulated file has an (optional) header followed by a
series of one or more sections, defined as follows (without the
parenthesized comments). The header must start with a `#` character

1
tools/README
View File

@ -48,6 +48,7 @@ smd convert Smooth Mach Dynamics triangles to VTK
spin perform a cubic polynomial interpolation of a GNEB MEP
stl_bin2txt convert binary STL files to ASCII
swig Interface file and demo scripts for SWIG wrappers for the LAMMPS C library interface
tabulate Python scripts to generate tabulated potential files for LAMMPS
tinker2lmp.py convert Tinker input files to LAMMPS input files
valgrind suppression files for use with valgrind's memcheck tool
vim add-ons to VIM editor for editing LAMMPS input scripts

39
tools/tabulate/README.md Normal file
View File

@ -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)
```

21
tools/tabulate/angle_harmonic_tabulate.py Executable file
View File

@ -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')

28
tools/tabulate/bond_morse_tabulate.py Executable file
View File

@ -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')

19
tools/tabulate/dihedral_harmonic_tabulate.py Executable file
View File

@ -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')

31
tools/tabulate/pair_hybrid_tabulate.py Executable file
View File

@ -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')

26
tools/tabulate/pair_lj_tabulate.py Executable file
View File

@ -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')

279
tools/tabulate/tabulate.py Executable file
View File

@ -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()