From ec87c71064ca25172d5b6385b7c4093169434b4d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 23 Feb 2023 00:28:58 -0500 Subject: [PATCH] add documentation for fix wall/lepton, fix wall/table and WallTabulate() --- doc/src/Commands_fix.rst | 2 + doc/src/Tools.rst | 3 +- doc/src/fix.rst | 2 + doc/src/fix_wall.rst | 257 +++++++++++++++----- doc/utils/sphinx-config/false_positives.txt | 1 + tools/tabulate/README.md | 18 +- tools/tabulate/wall_multi_tabulate.py | 10 +- 7 files changed, 224 insertions(+), 69 deletions(-) diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 63adc56ed3..5c164c8f0a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -263,6 +263,7 @@ OPT. * :doc:`wall/lj1043 ` * :doc:`wall/lj126 ` * :doc:`wall/lj93 (k) ` + * :doc:`wall/lepton ` * :doc:`wall/morse ` * :doc:`wall/piston ` * :doc:`wall/reflect (k) ` @@ -270,4 +271,5 @@ OPT. * :doc:`wall/region ` * :doc:`wall/region/ees ` * :doc:`wall/srd ` + * :doc:`wall/table ` * :doc:`widom ` diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index 20e3d2c094..2eae1aabc7 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -320,7 +320,8 @@ eam generate tool ----------------------------- The tools/eam_generate directory contains several one-file C programs -that convert an analytic formula into a tabulated :doc:`embedded atom method (EAM) ` setfl potential file. The potentials they +that convert an analytic formula into a tabulated :doc:`embedded atom +method (EAM) ` setfl potential file. The potentials they produce are in the potentials directory, and can be used with the :doc:`pair_style eam/alloy ` command. diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 0a926d570c..ec5922163c 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -415,6 +415,7 @@ accelerated styles exist. * :doc:`wall/lj1043 ` - Lennard-Jones 10--4--3 wall * :doc:`wall/lj126 ` - Lennard-Jones 12--6 wall * :doc:`wall/lj93 ` - Lennard-Jones 9--3 wall +* :doc:`wall/lepton ` - Custom Lepton expression wall * :doc:`wall/morse ` - Morse potential wall * :doc:`wall/piston ` - moving reflective piston wall * :doc:`wall/reflect ` - reflecting wall(s) @@ -422,6 +423,7 @@ accelerated styles exist. * :doc:`wall/region ` - use region surface as wall * :doc:`wall/region/ees ` - use region surface as wall for ellipsoidal particles * :doc:`wall/srd ` - slip/no-slip wall for SRD particles +* :doc:`wall/table ` - Tabulated potential wall wall * :doc:`widom ` - Widom insertions of atoms or molecules Restrictions diff --git a/doc/src/fix_wall.rst b/doc/src/fix_wall.rst index 9f20470bbc..7b38b6035c 100644 --- a/doc/src/fix_wall.rst +++ b/doc/src/fix_wall.rst @@ -4,7 +4,9 @@ .. index:: fix wall/lj1043 .. index:: fix wall/colloid .. index:: fix wall/harmonic +.. index:: fix wall/lepton .. index:: fix wall/morse +.. index:: fix wall/table fix wall/lj93 command ===================== @@ -23,20 +25,31 @@ fix wall/colloid command fix wall/harmonic command ========================= +fix wall/lepton command +========================= + fix wall/morse command ====================== +fix wall/table command +========================= + Syntax """""" .. parsed-literal:: - fix ID group-ID style face args ... keyword value ... + fix ID group-ID style [tabstyle] [N] face args ... keyword value ... * ID, group-ID are documented in :doc:`fix ` command -* style = *wall/lj93* or *wall/lj126* or *wall/lj1043* or *wall/colloid* or *wall/harmonic* or *wall/morse* +* style = *wall/lj93* or *wall/lj126* or *wall/lj1043* or *wall/colloid* or *wall/harmonic* or *wall/lepton* or *wall/morse* or *wall/table* +* tabstyle = *linear* or *spline* = method of table interpolation (only applies to *wall/table*) +* N = use N values in *linear* or *spline* interpolation (only applies to *wall/table*) * one or more face/arg pairs may be appended * face = *xlo* or *xhi* or *ylo* or *yhi* or *zlo* or *zhi* + +.. spacer + * args for styles *lj93* or *lj126* or *lj1043* or *colloid* or *harmonic* .. parsed-literal:: @@ -50,7 +63,19 @@ Syntax epsilon can be a variable (see below) sigma = size factor for wall-particle interaction (distance units) sigma can be a variable (see below) - cutoff = distance from wall at which wall-particle interaction is cut off (distance units) + cutoff = distance from wall at which wall-particle interactions are cut off (distance units) + +* args for style *lepton* + + .. parsed-literal:: + + args = coord expression cutoff + coord = position of wall = EDGE or constant or variable + EDGE = current lo or hi edge of simulation box + constant = number like 0.0 or -30.0 (distance units) + variable = :doc:`equal-style variable ` like v_x or v_wiggle + expression = Lepton expression for the potential (energy units) + cutoff = distance from wall at which wall-particle interactions are cut off (distance units) * args for style *morse* @@ -67,7 +92,20 @@ Syntax alpha can be a variable (see below) r_0 = distance of the potential minimum from the face of region (distance units) r_0 can be a variable (see below) - cutoff = distance from wall at which wall-particle interaction is cut off (distance units) + cutoff = distance from wall at which wall-particle interactions are cut off (distance units) + +* args for style *table* + + .. parsed-literal:: + + args = coord filename keyword cutoff + coord = position of wall = EDGE or constant or variable + EDGE = current lo or hi edge of simulation box + constant = number like 0.0 or -30.0 (distance units) + variable = :doc:`equal-style variable ` like v_x or v_wiggle + filename = file containing tabulated energy and force values + keyword = section identifier to select a specific table in table file + cutoff = distance from wall at which wall-particle interactions are cut off (distance units) * zero or more keyword/value pairs may be appended * keyword = *units* or *fld* or *pbc* @@ -91,9 +129,13 @@ Examples fix wallhi all wall/lj93 xlo -1.0 1.0 1.0 2.5 units box fix wallhi all wall/lj93 xhi EDGE 1.0 1.0 2.5 + fix wallhi all wall/harmonic xhi EDGE 100.0 0.0 4.0 units box fix wallhi all wall/morse xhi EDGE 1.0 1.0 1.0 2.5 units box fix wallhi all wall/lj126 v_wiggle 23.2 1.0 1.0 2.5 fix zwalls all wall/colloid zlo 0.0 1.0 1.0 0.858 zhi 40.0 1.0 1.0 0.858 + fix xwall mobile wall/table spline 200 EDGE -5.0 walltab.dat HARMONIC 4.0 + fix xwalls mobile wall/lepton xlo -5.0 "k*(r-rc)^2;k=100.0" 4.0 xhi 5.0 "k*(r-rc)^2;k=100.0" 4.0 + Description """"""""""" @@ -103,7 +145,7 @@ wall that interacts with the atoms in the group by generating a force on the atom in a direction perpendicular to the wall. The energy of wall-particle interactions depends on the style. -For style *wall/lj93*, the energy E is given by the 9/3 potential: +For style *wall/lj93*, the energy E is given by the 9-3 Lennard-Jones potential: .. math:: @@ -111,7 +153,7 @@ For style *wall/lj93*, the energy E is given by the 9/3 potential: \left(\frac{\sigma}{r}\right)^3 \right] \qquad r < r_c -For style *wall/lj126*, the energy E is given by the 12/6 potential: +For style *wall/lj126*, the energy E is given by the 12-6 Lennard-Jones potential: .. math:: @@ -119,7 +161,7 @@ For style *wall/lj126*, the energy E is given by the 12/6 potential: \left(\frac{\sigma}{r}\right)^6 \right] \qquad r < r_c -For style *wall/lj1043*, the energy E is given by the 10/4/3 potential: +For style *wall/lj1043*, the energy E is given by the 10-4-3 Lennard-Jones potential: .. math:: @@ -138,8 +180,8 @@ of the :doc:`pair_style colloid ` potential: & \left. - \frac{1}{6} \left(\frac{2R(D+R) + D(D+2R) \left[ \ln D - \ln (D+2R) \right]}{D(D+2R)} \right) \right] \qquad r < r_c -For style *wall/harmonic*, the energy E is given by a harmonic spring -potential: +For style *wall/harmonic*, the energy E is given by a repulsive-only harmonic +spring potential: .. math:: @@ -152,6 +194,56 @@ For style *wall/morse*, the energy E is given by a Morse potential: E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right] \qquad r < r_c +For style *wall/lepton*, the energy E is provided as an Lepton +expression string using "r" as the distance variable. The `Lepton +library `_, that the *wall/lepton* +style interfaces with, evaluates this expression string at run time to +compute the wall-particle energy. It also creates an analytical +representation of the first derivative of this expression with respect +to "r" and then uses that to compute the force between the wall and +atoms in the fix group. The Lepton expression must be either enclosed +in quotes or must not contain any whitespace so that LAMMPS recognizes +it as a single keyword. + +Optionally, the expression may use "rc" to refer to the cutoff distance +for the given wall. Further constants in the expression can be defined +in the same string as additional expressions separated by semi-colons. +The expression "k*(r-rc)^2;k=100.0" represents a repulsive-only harmonic +spring as in fix *wall/harmonic* with a force constant *K* (same as +:math:`\epsilon` above) of 100 energy units. More details on the Lepton +expression strings are given below. + +For style *wall/table*, the energy E and forces are determined from +interpolation tables listed in one or more files as a function of +distance. The interpolation tables are used to evaluate energy and +forces between particles and the wall similar to how analytic formulas +are used for the other wall styles. + +The interpolation tables are created as a pre-computation by fitting +cubic splines to the file values and interpolating energy and force +values at each of *N* distances. During a simulation, the tables are +used to interpolate energy and force values as needed for each wall and +particle separated by a distance *R*\ . The interpolation is done in +one of two styles: *linear* or *spline*\ . + +For the *linear* style, the distance *R* is used to find the 2 +surrounding table values from which an energy or force is computed by +linear interpolation. + +For the *spline* style, cubic spline coefficients are computed and +stored for each of the *N* values in the table, one set of splines for +energy, another for force. Note that these splines are different than +the ones used to pre-compute the *N* values. Those splines were fit +to the *Nfile* values in the tabulated file, where often *Nfile* < +*N*\ . The distance *R* is used to find the appropriate set of spline +coefficients which are used to evaluate a cubic polynomial which +computes the energy or force. + +For each wall a filename and a keyword must be provided as in the +examples above. The filename specifies a file containing tabulated +energy and force values. The keyword specifies a section of the file. +The format of this file is described below. + In all cases, *r* is the distance from the particle to the wall at position *coord*, and :math:`r_c` is the *cutoff* distance at which the particle and wall no longer interact. The energy of the wall @@ -180,11 +272,12 @@ box parameters and timestep and elapsed time. Thus it is easy to specify a time-dependent wall position. See examples below. For the *wall/lj93* and *wall/lj126* and *wall/lj1043* styles, -:math:`\epsilon` and :math:`\sigma` are the usual Lennard-Jones parameters, which -determine the strength and size of the particle as it interacts with -the wall. Epsilon has energy units. Note that this :math:`\epsilon` and -:math:`\sigma` may be different than any :math:`\epsilon` or :math:`\sigma` values defined -for a pair style that computes particle-particle interactions. +:math:`\epsilon` and :math:`\sigma` are the usual Lennard-Jones +parameters, which determine the strength and size of the particle as it +interacts with the wall. Epsilon has energy units. Note that this +:math:`\epsilon` and :math:`\sigma` may be different than any +:math:`\epsilon` or :math:`\sigma` values defined for a pair style that +computes particle-particle interactions. The *wall/lj93* interaction is derived by integrating over a 3d half-lattice of Lennard-Jones 12/6 particles. The *wall/lj126* @@ -207,11 +300,11 @@ are the number density of the constituent particles, in the wall and colloid respectively, in units of 1/volume. The *wall/colloid* interaction is derived by integrating over -constituent LJ particles of size :math:`\sigma` within the colloid particle -and a 3d half-lattice of Lennard-Jones 12/6 particles of size :math:`\sigma` -in the wall. As mentioned in the preceding paragraph, the density of -particles in the wall and colloid can be different, as specified by -the :math:`\epsilon` prefactor. +constituent LJ particles of size :math:`\sigma` within the colloid +particle and a 3d half-lattice of Lennard-Jones 12/6 particles of size +:math:`\sigma` in the wall. As mentioned in the preceding paragraph, +the density of particles in the wall and colloid can be different, as +specified by the :math:`\epsilon` prefactor. For the *wall/harmonic* style, :math:`\epsilon` is effectively the spring constant K, and has units (energy/distance\^2). The input parameter @@ -220,20 +313,21 @@ spring is at the *cutoff*\ . This is a repulsive-only spring since the interaction is truncated at the *cutoff* For the *wall/morse* style, the three parameters are in this order: -:math:`D_0` the depth of the potential, :math:`\alpha` the width parameter, and -:math:`r_0` the location of the minimum. :math:`D_0` has energy units, :math:`\alpha` -inverse distance units, and :math:`r_0` distance units. +:math:`D_0` the depth of the potential, :math:`\alpha` the width +parameter, and :math:`r_0` the location of the minimum. :math:`D_0` has +energy units, :math:`\alpha` inverse distance units, and :math:`r_0` +distance units. -For any wall, the :math:`\epsilon` and/or :math:`\sigma` and/or :math:`\alpha` parameter can -be specified -as an :doc:`equal-style variable `, in which case it should be +For any wall that supports them, the :math:`\epsilon` and/or +:math:`\sigma` and/or :math:`\alpha` parameter can be specified as an +:doc:`equal-style variable `, in which case it should be specified as v_name, where name is the variable name. As with a -variable wall position, the variable is evaluated each timestep and -the result becomes the current epsilon or sigma of the wall. -Equal-style variables can specify formulas with various mathematical -functions, and include :doc:`thermo_style ` command -keywords for the simulation box parameters and timestep and elapsed -time. Thus it is easy to specify a time-dependent wall interaction. +variable wall position, the variable is evaluated each timestep and the +result becomes the current epsilon or sigma of the wall. Equal-style +variables can specify formulas with various mathematical functions, and +include :doc:`thermo_style ` command keywords for the +simulation box parameters and timestep and elapsed time. Thus it is +easy to specify a time-dependent wall interaction. .. note:: @@ -266,20 +360,19 @@ define the lattice spacings. The *fld* keyword can be used with a *yes* setting to invoke the wall constraint before pairwise interactions are computed. This allows an implicit FLD model using :doc:`pair_style lubricateU ` -to include the wall force in its calculations. If the setting is -*no*, wall forces are imposed after pairwise interactions, in the -usual manner. +to include the wall force in its calculations. If the setting is *no*, +wall forces are imposed after pairwise interactions, in the usual +manner. -The *pbc* keyword can be used with a *yes* setting to allow walls to -be specified in a periodic dimension. See the -:doc:`boundary ` command for options on simulation box -boundaries. The default for *pbc* is *no*, which means the system -must be non-periodic when using a wall. But you may wish to use a -periodic box. E.g. to allow some particles to interact with the wall -via the fix group-ID, and others to pass through it and wrap around a -periodic box. In this case you should ensure that the wall if -sufficiently far enough away from the box boundary. If you do not, -then particles may interact with both the wall and with periodic +The *pbc* keyword can be used with a *yes* setting to allow walls to be +specified in a periodic dimension. See the :doc:`boundary ` +command for options on simulation box boundaries. The default for *pbc* +is *no*, which means the system must be non-periodic when using a wall. +But you may wish to use a periodic box. E.g. to allow some particles to +interact with the wall via the fix group-ID, and others to pass through +it and wrap around a periodic box. In this case you should ensure that +the wall if sufficiently far enough away from the box boundary. If you +do not, then particles may interact with both the wall and with periodic images on the other side of the box, which is probably not what you want. @@ -328,6 +421,57 @@ perturbation on the particles: ---------- +.. include:: lepton_expression.rst + +---------- + +Table file format +""""""""""""""""" + +Suitable tables for use with fix *wall/table* can be created by 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): + +.. parsed-literal:: + + # Tabulated wall potential UNITS: real + + HARMONIC (keyword is the first text on a line) + N 100 FP 200 200 + (blank line) + 1 0.04 1568.16 792.00 (index, distance to wall, energy, force) + 2 0.08 1536.64 784.00 + 3 0.12 1505.44 776.00 + ... + 99 3.96 0.16 8.00 + 100 4.00 0 0 + +A section begins with a non-blank line whose first character is not a +"#"; blank lines or lines starting with "#" can be used as comments +between sections. The first line begins with a keyword which identifies +the section. The line can contain additional text, but the initial text +must match the argument specified in the fix *wall/table* command. The +next line lists (in any order) one or more parameters for the table. +Each parameter is a keyword followed by one or more numeric values. + +The parameter "N" is required and its value is the number of table +entries that follow. Note that this may be different than the *N* +specified in the fix *wall/table* command. Let Ntable = *N* in the fix +command, and Nfile = "N" in the tabulated file. What LAMMPS does is a +preliminary interpolation by creating splines using the Nfile tabulated +values as nodal points. It uses these to interpolate as needed to +generate energy and force values at Ntable different points. The +resulting tables of length Ntable are then used as described above, when +computing energy and force for wall-particle interactions. This means that +if you want the interpolation tables of length Ntable to match exactly +what is in the tabulated file (with effectively no preliminary +interpolation), you should set Ntable = Nfile. + +---------- + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -354,16 +498,15 @@ fix. This allows to set at which level of the :doc:`r-RESPA ` integrator the fix is adding its forces. Default is the outermost level. -This fix computes a global scalar energy and a global vector of -forces, which can be accessed by various :doc:`output commands -`. Note that the scalar energy is the sum of -interactions with all defined walls. If you want the energy on a -per-wall basis, you need to use multiple fix wall commands. The -length of the vector is equal to the number of walls defined by the -fix. Each vector value is the normal force on a specific wall. Note -that an outward force on a wall will be a negative value for *lo* -walls and a positive value for *hi* walls. The scalar and vector -values calculated by this fix are "extensive". +This fix computes a global scalar energy and a global vector of forces, +which can be accessed by various :doc:`output commands `. +Note that the scalar energy is the sum of interactions with all defined +walls. If you want the energy on a per-wall basis, you need to use +multiple fix wall commands. The length of the vector is equal to the +number of walls defined by the fix. Each vector value is the normal +force on a specific wall. Note that an outward force on a wall will be +a negative value for *lo* walls and a positive value for *hi* walls. +The scalar and vector values calculated by this fix are "extensive". No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. @@ -386,7 +529,11 @@ invoked by the :doc:`minimize ` command. Restrictions """""""""""" - none + +Fix *wall/lepton* is part of the LEPTON package and only enabled if +LAMMPS was built with this package. See the :doc:`Build package +` page for more info. + Related commands """""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index fdacbf9c4c..d03300439f 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -3487,6 +3487,7 @@ sz Sz Tabbernor tabinner +tabstyle Tadmor Tafipolsky tagID diff --git a/tools/tabulate/README.md b/tools/tabulate/README.md index fb23e260a1..6f63d03a80 100644 --- a/tools/tabulate/README.md +++ b/tools/tabulate/README.md @@ -6,13 +6,14 @@ 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. +BondTabulate(), AngleTabulate(), DihedralTabulate(), or WallTabulate() +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 from the energy function through +numerical differentiation. Please see the individual tabulation scripts in this folder for examples: @@ -24,7 +25,8 @@ Please see the individual tabulation scripts in this folder for examples: | 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 | -| pair_zbladd_tabulate.py | creates a table for hybrid/overlay to use ZBL repulsion with an EAM potential | +| 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 | | ------------------------------|-------------------------------------------------------------------------------| Common command line flags: diff --git a/tools/tabulate/wall_multi_tabulate.py b/tools/tabulate/wall_multi_tabulate.py index 0458a929e2..b1b4aa6cbf 100755 --- a/tools/tabulate/wall_multi_tabulate.py +++ b/tools/tabulate/wall_multi_tabulate.py @@ -11,7 +11,7 @@ epsilon = 0.02 sigma = 2.0 depth = 20.0 width = 2.0 -rzero = 1.2 +r0 = 1.2 def harmonic_force(r): dr = r - rzero @@ -32,12 +32,12 @@ def lj126_energy(r): return f def morse_energy(r): - ralpha = math.exp(-width*(r-rzero)) + ralpha = math.exp(-width*(r-r0)) f = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha)) return f def morse_force(r): - ralpha = math.exp(-width*(r-rzero)) + ralpha = math.exp(-width*(r-r0)) f = -2.0 * depth * width * (1.0 -ralpha) * ralpha return f @@ -51,9 +51,9 @@ if __name__ == "__main__": sys.argv.append('--filename') sys.argv.append(fname) sys.argv.append('--num-points') - sys.argv.append('1000') + sys.argv.append('100') sys.argv.append('--inner') - sys.argv.append('0.02') + sys.argv.append('0.04') sys.argv.append('--outer') sys.argv.append('4.0') wtable = WallTabulate(harmonic_energy, harmonic_force, units='real')