Merge pull request #3279 from schererc/develop
Addition of two new MANYBODY pair styles (sw/angle/table and threebody/table)
This commit is contained in:
Axel Kohlmeyer
@ -269,6 +269,7 @@ OPT.
|
||||
* :doc:`spin/neel <pair_spin_neel>`
|
||||
* :doc:`srp <pair_srp>`
|
||||
* :doc:`sw (giko) <pair_sw>`
|
||||
* :doc:`sw/angle/table <pair_sw_angle_table>`
|
||||
* :doc:`sw/mod (o) <pair_sw>`
|
||||
* :doc:`table (gko) <pair_table>`
|
||||
* :doc:`table/rx (k) <pair_table_rx>`
|
||||
@ -279,6 +280,7 @@ OPT.
|
||||
* :doc:`tersoff/table (o) <pair_tersoff>`
|
||||
* :doc:`tersoff/zbl (gko) <pair_tersoff_zbl>`
|
||||
* :doc:`thole <pair_thole>`
|
||||
* :doc:`threebody/table <pair_threebody_table>`
|
||||
* :doc:`tip4p/cut (o) <pair_coul>`
|
||||
* :doc:`tip4p/long (o) <pair_coul>`
|
||||
* :doc:`tip4p/long/soft (o) <pair_fep_soft>`
|
||||
|
||||
@ -348,6 +348,7 @@ accelerated styles exist.
|
||||
* :doc:`spin/neel <pair_spin_neel>` -
|
||||
* :doc:`srp <pair_srp>` -
|
||||
* :doc:`sw <pair_sw>` - Stillinger-Weber 3-body potential
|
||||
* :doc:`sw/angle/table <pair_sw_angle_table>` - Stillinger-Weber potential with tabulated angular term
|
||||
* :doc:`sw/mod <pair_sw>` - modified Stillinger-Weber 3-body potential
|
||||
* :doc:`table <pair_table>` - tabulated pair potential
|
||||
* :doc:`table/rx <pair_table_rx>` -
|
||||
@ -358,6 +359,7 @@ accelerated styles exist.
|
||||
* :doc:`tersoff/table <pair_tersoff>` -
|
||||
* :doc:`tersoff/zbl <pair_tersoff_zbl>` - Tersoff/ZBL 3-body potential
|
||||
* :doc:`thole <pair_thole>` - Coulomb interactions with thole damping
|
||||
* :doc:`threebody/table <pair_threebody_table>` - generic tabulated three-body potential
|
||||
* :doc:`tip4p/cut <pair_coul>` - Coulomb for TIP4P water w/out LJ
|
||||
* :doc:`tip4p/long <pair_coul>` - long-range Coulomb for TIP4P water w/out LJ
|
||||
* :doc:`tip4p/long/soft <pair_fep_soft>` -
|
||||
|
||||
311
doc/src/pair_sw_angle_table.rst
Normal file
311
doc/src/pair_sw_angle_table.rst
Normal file
@ -0,0 +1,311 @@
|
||||
.. index:: pair_style sw/angle/table
|
||||
|
||||
pair_style sw/angle/table command
|
||||
=================================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style style
|
||||
|
||||
* style = *sw/angle/table*
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style sw/angle/table
|
||||
pair_coeff * * spce.sw type
|
||||
|
||||
Used in example input script:
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
examples/PACKAGES/manybody_table/in.spce_sw
|
||||
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The *sw/angle/table* style is a modification of the original
|
||||
:doc:`pair_style sw <pair_sw>`. It has been developed for coarse-grained
|
||||
simulations (of water) (:ref:`Scherer1 <Scherer1>`), but can be employed
|
||||
for all kinds of systems. It computes a modified 3-body
|
||||
:ref:`Stillinger-Weber <Stillinger3>` potential for the energy E of a
|
||||
system of atoms as
|
||||
|
||||
.. math::
|
||||
|
||||
E & = \sum_i \sum_{j > i} \phi_2 (r_{ij}) +
|
||||
\sum_i \sum_{j \neq i} \sum_{k > j}
|
||||
\phi_3 (r_{ij}, r_{ik}, \theta_{ijk}) \\
|
||||
\phi_2(r_{ij}) & = A_{ij} \epsilon_{ij} \left[ B_{ij} (\frac{\sigma_{ij}}{r_{ij}})^{p_{ij}} -
|
||||
(\frac{\sigma_{ij}}{r_{ij}})^{q_{ij}} \right]
|
||||
\exp \left( \frac{\sigma_{ij}}{r_{ij} - a_{ij} \sigma_{ij}} \right) \\
|
||||
\phi_3(r_{ij},r_{ik},\theta_{ijk}) & = f^{\textrm{3b}}\left(\theta_{ijk}\right)
|
||||
\exp \left( \frac{\gamma_{ij} \sigma_{ij}}{r_{ij} - a_{ij} \sigma_{ij}} \right)
|
||||
\exp \left( \frac{\gamma_{ik} \sigma_{ik}}{r_{ik} - a_{ik} \sigma_{ik}} \right)
|
||||
|
||||
where :math:`\phi_2` is a two-body term and :math:`\phi_3` is a
|
||||
three-body term. The summations in the formula are over all neighbors J
|
||||
and K of atom I within a cutoff distance :math:`a \sigma`. In contrast
|
||||
to the original *sw* style, *sw/angle/table* allows for a flexible
|
||||
three-body term :math:`f^{\textrm{3b}}\left(\theta_{ijk}\right)` which
|
||||
is read in as a tabulated interaction. It can be parameterized with the
|
||||
csg_fmatch app of VOTCA as available at:
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca.
|
||||
|
||||
Only a single pair_coeff command is used with the *sw/angle/table* style
|
||||
which specifies a modified Stillinger-Weber potential file with
|
||||
parameters for all needed elements. These are mapped to LAMMPS atom
|
||||
types by specifying N_el additional arguments after the ".sw" filename
|
||||
in the pair_coeff command, where N_el is the number of LAMMPS atom
|
||||
types:
|
||||
|
||||
* ".sw" filename
|
||||
* N_el element names = mapping of SW elements to atom types
|
||||
|
||||
See the :doc:`pair_coeff <pair_coeff>` page for alternate ways to
|
||||
specify the path for the potential file.
|
||||
|
||||
As an example, imagine a file SiC.sw has Stillinger-Weber values for Si
|
||||
and C. If your LAMMPS simulation has 4 atoms types and you want the
|
||||
first 3 to be Si, and the fourth to be C, you would use the following
|
||||
pair_coeff command:
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_coeff * * SiC.sw Si Si Si C
|
||||
|
||||
The first 2 arguments must be \* \* so as to span all LAMMPS atom types.
|
||||
The first three Si arguments map LAMMPS atom types 1,2,3 to the Si
|
||||
element in the SW file. The final C argument maps LAMMPS atom type 4 to
|
||||
the C element in the SW file. If a mapping value is specified as NULL,
|
||||
the mapping is not performed. This can be used when a *sw/angle/table*
|
||||
potential is used as part of the *hybrid* pair style. The NULL values
|
||||
are placeholders for atom types that will be used with other potentials.
|
||||
|
||||
The (modified) Stillinger-Weber files have a ".sw" suffix. Lines that
|
||||
are not blank or comments (starting with #) define parameters for a
|
||||
triplet of elements. The parameters in a single entry correspond to the
|
||||
two-body and three-body coefficients in the formula above. Here, also
|
||||
the suffix ".sw" is used though the original Stillinger-Weber file
|
||||
format is supplemented with four additional lines per parameter block to
|
||||
specify the tabulated three-body interaction. A single entry then
|
||||
contains:
|
||||
|
||||
* element 1 (the center atom in a 3-body interaction)
|
||||
* element 2
|
||||
* element 3
|
||||
* :math:`\epsilon` (energy units)
|
||||
* :math:`\sigma` (distance units)
|
||||
* a
|
||||
* :math:`\lambda`
|
||||
* :math:`\gamma`
|
||||
* :math:`\cos\theta_0`
|
||||
* A
|
||||
* B
|
||||
* p
|
||||
* q
|
||||
* tol
|
||||
* filename
|
||||
* keyword
|
||||
* style
|
||||
* N
|
||||
|
||||
The A, B, p, and q parameters are used only for two-body interactions.
|
||||
The :math:`\lambda` and :math:`\cos\theta_0` parameters, only used for
|
||||
three-body interactions in the original Stillinger-Weber style, are read
|
||||
in but ignored in this modified pair style. The :math:`\epsilon`
|
||||
parameter is only used for two-body interactions in this modified pair
|
||||
style and not for the three-body terms. The :math:`\sigma` and *a*
|
||||
parameters are used for both two-body and three-body
|
||||
interactions. :math:`\gamma` is used only in the three-body
|
||||
interactions, but is defined for pairs of atoms. The non-annotated
|
||||
parameters are unitless.
|
||||
|
||||
LAMMPS introduces an additional performance-optimization parameter tol
|
||||
that is used for both two-body and three-body interactions. In the
|
||||
Stillinger-Weber potential, the interaction energies become negligibly
|
||||
small at atomic separations substantially less than the theoretical
|
||||
cutoff distances. LAMMPS therefore defines a virtual cutoff distance
|
||||
based on a user defined tolerance tol. The use of the virtual cutoff
|
||||
distance in constructing atom neighbor lists can significantly reduce
|
||||
the neighbor list sizes and therefore the computational cost. LAMMPS
|
||||
provides a *tol* value for each of the three-body entries so that they
|
||||
can be separately controlled. If tol = 0.0, then the standard
|
||||
Stillinger-Weber cutoff is used.
|
||||
|
||||
The additional parameters *filename*, *keyword*, *style*, and *N* refer
|
||||
to the tabulated angular potential
|
||||
:math:`f^{\textrm{3b}}\left(\theta_{ijk}\right)`. The tabulated angular
|
||||
potential has to be of the format as used in the :doc:`angle_style table
|
||||
<angle_table>` command:
|
||||
|
||||
An interpolation tables of length *N* is created. The interpolation is
|
||||
done in one of 2 *styles*: *linear* or *spline*. For the *linear*
|
||||
style, the angle is used to find 2 surrounding table values from which
|
||||
an energy or its derivative is computed by linear interpolation. For the
|
||||
*spline* style, a cubic spline coefficients are computed and stored at
|
||||
each of the *N* values in the table. The angle is used to find the
|
||||
appropriate set of coefficients which are used to evaluate a cubic
|
||||
polynomial which computes the energy or derivative.
|
||||
|
||||
The *filename* specifies the file containing the tabulated energy and
|
||||
derivative values of :math:`f^{\textrm{3b}}\left(\theta_{ijk}\right)`.
|
||||
The *keyword* then specifies a section of the file. The format of this
|
||||
file is as follows (without the parenthesized comments):
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
# Angle potential for harmonic (one or more comment or blank lines)
|
||||
|
||||
HAM (keyword is the first text on line)
|
||||
N 181 FP 0 0 EQ 90.0 (N, FP, EQ parameters)
|
||||
(blank line)
|
||||
1 0.0 200.5 2.5 (index, angle, energy, derivative)
|
||||
2 1.0 198.0 2.5
|
||||
...
|
||||
181 180.0 0.0 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 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 Stillinger-Weber potential file. Let Nsw = *N* in the
|
||||
".sw" file, and Nfile = "N" in the tabulated angular 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 derivative values at Ntable different
|
||||
points. The resulting tables of length Nsw are then used as described
|
||||
above, when computing energy and force for individual angles and their
|
||||
atoms. This means that if you want the interpolation tables of length
|
||||
Nsw to match exactly what is in the tabulated file (with effectively no
|
||||
preliminary interpolation), you should set Nsw = Nfile.
|
||||
|
||||
The "FP" parameter is optional. If used, it is followed by two values
|
||||
fplo and fphi, which are the second derivatives at the innermost and
|
||||
outermost angle settings. These values are needed by the spline
|
||||
construction routines. If not specified by the "FP" parameter, they are
|
||||
estimated (less accurately) by the first two and last two derivative
|
||||
values in the table.
|
||||
|
||||
The "EQ" parameter is also optional. If used, it is followed by a the
|
||||
equilibrium angle value, which is used, for example, by the :doc:`fix
|
||||
shake <fix_shake>` command. If not used, the equilibrium angle is set to
|
||||
180.0.
|
||||
|
||||
Following a blank line, the next N lines of the angular table file list
|
||||
the tabulated values. On each line, the first value is the index from 1
|
||||
to N, the second value is the angle value (in degrees), the third value
|
||||
is the energy (in energy units), and the fourth is -dE/d(theta) (also in
|
||||
energy units). The third term is the energy of the 3-atom configuration
|
||||
for the specified angle. The last term is the derivative of the energy
|
||||
with respect to the angle (in degrees, not radians). Thus the units of
|
||||
the last term are still energy, not force. The angle values must
|
||||
increase from one line to the next. The angle values must also begin
|
||||
with 0.0 and end with 180.0, i.e. span the full range of possible
|
||||
angles.
|
||||
|
||||
Note that one angular potential file can contain many sections, each
|
||||
with a tabulated potential. LAMMPS reads the file section by section
|
||||
until it finds one that matches the specified *keyword* of appropriate
|
||||
section of the ".sw" file.
|
||||
|
||||
The Stillinger-Weber potential file must contain entries for all the
|
||||
elements listed in the pair_coeff command. It can also contain entries
|
||||
for additional elements not being used in a particular simulation;
|
||||
LAMMPS ignores those entries.
|
||||
|
||||
For a single-element simulation, only a single entry is required
|
||||
(e.g. SiSiSi). For a two-element simulation, the file must contain 8
|
||||
entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that
|
||||
specify SW parameters for all permutations of the two elements
|
||||
interacting in three-body configurations. Thus for 3 elements, 27
|
||||
entries would be required, etc.
|
||||
|
||||
As annotated above, the first element in the entry is the center atom in
|
||||
a three-body interaction. Thus an entry for SiCC means a Si atom with 2
|
||||
C atoms as neighbors. The parameter values used for the two-body
|
||||
interaction come from the entry where the second and third elements are
|
||||
the same. Thus the two-body parameters for Si interacting with C, comes
|
||||
from the SiCC entry. The three-body angular potential
|
||||
:math:`f^{\textrm{3b}}\left(\theta_{ijk}\right)` can in principle be
|
||||
specific to the three elements of the configuration. However, the user
|
||||
must ensure that it makes physically sense. Note also that the function
|
||||
:math:`\phi_3` contains two exponential screening factors with parameter
|
||||
values from the ij pair and ik pairs. So :math:`\phi_3` for a C atom
|
||||
bonded to a Si atom and a second C atom will depend on the three-body
|
||||
parameters for the CSiC entry, and also on the two-body parameters for
|
||||
the CCC and CSiSi entries. Since the order of the two neighbors is
|
||||
arbitrary, the three-body parameters and the tabulated angular potential
|
||||
for entries CSiC and CCSi should be the same. Similarly, the two-body
|
||||
parameters for entries SiCC and CSiSi should also be the same. The
|
||||
parameters used only for two-body interactions (A, B, p, and q) in
|
||||
entries whose second and third element are different (e.g. SiCSi) are
|
||||
not used for anything and can be set to 0.0 if desired. This is also
|
||||
true for the parameters in :math:`\phi_3` that are taken from the ij and
|
||||
ik pairs (:math:`\sigma`, *a*, :math:`\gamma`)
|
||||
|
||||
Additional input files and reference data can be found at:
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-tutorials/spce/3body_sw
|
||||
|
||||
----------
|
||||
|
||||
Mixing, shift, table, tail correction, restart, rRESPA info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
For atom type pairs I,J and I != J, where types I and J correspond to
|
||||
two different element types, mixing is performed by LAMMPS as described
|
||||
above from values in the potential file, but not for the tabulated
|
||||
angular potential file.
|
||||
|
||||
This pair style does not support the :doc:`pair_modify <pair_modify>`
|
||||
shift, table, and tail options.
|
||||
|
||||
This pair style does not write its information to :doc:`binary restart
|
||||
files <restart>`, since it is stored in potential files. Thus, you need
|
||||
to re-specify the pair_style and pair_coeff commands in an input script
|
||||
that reads a restart file.
|
||||
|
||||
This pair style can only be used via the *pair* keyword of the
|
||||
:doc:`run_style respa <run_style>` command. It does not support the
|
||||
*inner*, *middle*, *outer* keywords.
|
||||
|
||||
----------
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This pair style is part of the MANYBODY package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` page for more info.
|
||||
|
||||
This pair style requires the :doc:`newton <newton>` setting to be "on"
|
||||
for pair interactions.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style sw <pair_sw>`, :doc:`pair_style threebody/table <pair_threebody_table>`
|
||||
|
||||
|
||||
----------
|
||||
|
||||
.. _Stillinger3:
|
||||
|
||||
**(Stillinger)** Stillinger and Weber, Phys Rev B, 31, 5262 (1985).
|
||||
|
||||
.. _Scherer1:
|
||||
|
||||
**(Scherer1)** C. Scherer and D. Andrienko, Phys. Chem. Chem. Phys. 20, 22387-22394 (2018).
|
||||
|
||||
281
doc/src/pair_threebody_table.rst
Normal file
281
doc/src/pair_threebody_table.rst
Normal file
@ -0,0 +1,281 @@
|
||||
.. index:: pair_style threebody/table
|
||||
|
||||
pair_style threebody/table command
|
||||
==================================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style style
|
||||
|
||||
* style = *threebody/table*
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style threebody/table
|
||||
pair_coeff * * spce2.3b type1 type2
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce.3b type
|
||||
|
||||
Used in example input scripts:
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
examples/PACKAGES/manybody_table/in.spce
|
||||
examples/PACKAGES/manybody_table/in.spce2
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The *threebody/table* style is a pair style for generic tabulated
|
||||
three-body interactions. It has been developed for (coarse-grained)
|
||||
simulations (of water) with Kernel-based machine learning (ML)
|
||||
potentials (:ref:`Scherer2 <Scherer2>`). As for many other MANYBODY
|
||||
package pair styles the energy of a system is computed as a sum over
|
||||
three-body terms:
|
||||
|
||||
.. math::
|
||||
|
||||
E = \sum_i \sum_{j \neq i} \sum_{k > j} \phi_3 (r_{ij}, r_{ik}, \theta_{ijk})
|
||||
|
||||
The summations in the formula are over all neighbors J and K of atom I
|
||||
within a cutoff distance :math:`cut`. In contrast to the
|
||||
Stillinger-Weber potential, all forces are not calculated analytically,
|
||||
but read in from a three-body force/energy table which can be generated
|
||||
with the csg_ml app of VOTCA as available at:
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca.
|
||||
|
||||
Only a single pair_coeff command is used with the *threebody/table*
|
||||
style which specifies a threebody potential (".3b") file with parameters
|
||||
for all needed elements. These are then mapped to LAMMPS atom types by
|
||||
specifying N_el additional arguments after the ".3b" filename in the
|
||||
pair_coeff command, where N_el is the number of LAMMPS atom types:
|
||||
|
||||
* ".3b" filename
|
||||
* N_el element names = mapping of threebody elements to atom types
|
||||
|
||||
See the :doc:`pair_coeff <pair_coeff>` page for alternate ways to
|
||||
specify the path for the potential file.
|
||||
|
||||
As an example, imagine a file SiC.3b has three-body values for Si and C.
|
||||
If your LAMMPS simulation has 4 atoms types and you want the first 3 to
|
||||
be Si, and the fourth to be C, you would use the following pair_coeff
|
||||
command:
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_coeff * * SiC.3b Si Si Si C
|
||||
|
||||
The first 2 arguments must be \* \* so as to span all LAMMPS atom types.
|
||||
The first three Si arguments map LAMMPS atom types 1,2,3 to the Si
|
||||
element in the ".3b" file. The final C argument maps LAMMPS atom type 4
|
||||
to the C element in the threebody file. If a mapping value is specified
|
||||
as NULL, the mapping is not performed. This can be used when a
|
||||
*threebody/table* potential is used as part of the *hybrid* pair style.
|
||||
The NULL values are placeholders for atom types that will be used with
|
||||
other potentials.
|
||||
|
||||
The three-body files have a ".3b" suffix. Lines that are not blank or
|
||||
comments (starting with #) define parameters for a triplet of
|
||||
elements. The parameters in a single entry specify to the (three-body)
|
||||
cutoff distance and the tabulated three-body interaction. A single entry
|
||||
then contains:
|
||||
|
||||
* element 1 (the center atom in a 3-body interaction)
|
||||
* element 2
|
||||
* element 3
|
||||
* cut (distance units)
|
||||
* filename
|
||||
* keyword
|
||||
* style
|
||||
* N
|
||||
|
||||
The parameter :math:`cut` is the (three-body) cutoff distance. When set
|
||||
to 0, no interaction is calculated for this element triplet. The
|
||||
parameters *filename*, *keyword*, *style*, and *N* refer to the
|
||||
tabulated three-body potential.
|
||||
|
||||
The tabulation is done on a three-dimensional grid of the two distances
|
||||
:math:`r_{ij}` and :math:`r_{ik}` as well as the angle
|
||||
:math:`\theta_{ijk}` which is constructed in the following way. There
|
||||
are two different cases. If element 2 and element 3 are of the same
|
||||
type (e.g. SiCC), the distance :math:`r_{ij}` is varied in "N" steps
|
||||
from rmin to rmax and the distance :math:`r_{ik}` is varied from
|
||||
:math:`r_{ij}` to rmax. This can be done, due to the symmetry of the
|
||||
triplet. If element 2 and element 3 are not of the same type
|
||||
(e.g. SiCSi), there is no additional symmetry and the distance
|
||||
:math:`r_{ik}` is also varied from rmin to rmax in "N" steps. The angle
|
||||
:math:`\theta_{ijk}` is always varied in "2N" steps from (0.0 +
|
||||
180.0/(4N)) to (180.0 - 180.0/(4N)). Therefore, the total number of
|
||||
table entries is "M = N * N * (N+1)" for the symmetric (element 2 and
|
||||
element 3 are of the same type) and "M = 2 * N * N * N" for the general
|
||||
case (element 2 and element 3 are not of the same type).
|
||||
|
||||
The forces on all three particles I, J, and K of a triplet of this type
|
||||
of three-body interaction potential (:math:`\phi_3 (r_{ij}, r_{ik},
|
||||
\theta_{ijk})`) lie within the plane defined by the three inter-particle
|
||||
distance vectors :math:`{\mathbf r}_{ij}`, :math:`{\mathbf r}_{ik}`, and
|
||||
:math:`{\mathbf r}_{jk}`. This property is used to project the forces
|
||||
onto the inter-particle distance vectors as follows
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{pmatrix}
|
||||
{\mathbf f}_{i} \\
|
||||
{\mathbf f}_{j} \\
|
||||
{\mathbf f}_{k} \\
|
||||
\end{pmatrix} =
|
||||
\begin{pmatrix}
|
||||
f_{i1} & f_{i2} & 0 \\
|
||||
f_{j1} & 0 & f_{j2} \\
|
||||
0 & f_{k1} & f_{k2} \\
|
||||
\end{pmatrix}
|
||||
\begin{pmatrix}
|
||||
{\mathbf r}_{ij} \\
|
||||
{\mathbf r}_{ik} \\
|
||||
{\mathbf r}_{jk} \\
|
||||
\end{pmatrix}
|
||||
|
||||
and then tabulate the 6 force constants :math:`f_{i1}`, :math:`f_{i2}`,
|
||||
:math:`f_{j1}`, :math:`f_{j2}`, :math:`f_{k1}`, and :math:`f_{k2}`, as
|
||||
well as the energy of a triplet e. Due to symmetry reasons, the
|
||||
following relations hold: :math:`f_{i1}=-f_{j1}`,
|
||||
:math:`f_{i2}=-f_{k1}`, and :math:`f_{j2}=-f_{k2}`. As in this pair
|
||||
style the forces are read in directly, a correct MD simulation is also
|
||||
performed in the case that the triplet energies are set to e=0.
|
||||
|
||||
The *filename* specifies the file containing the tabulated energy and
|
||||
derivative values of :math:`\phi_3 (r_{ij}, r_{ik}, \theta_{ijk})`. The
|
||||
*keyword* then specifies a section of the file. The format of this file
|
||||
is as follows (without the parenthesized comments):
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
# Tabulated three-body potential for spce water (one or more comment or blank lines)
|
||||
|
||||
ENTRY1 (keyword is the first text on line)
|
||||
N 12 rmin 2.55 rmax 3.65 (N, rmin, rmax parameters)
|
||||
(blank line)
|
||||
1 2.55 2.55 3.75 -867.212 -611.273 867.212 21386.8 611.273 -21386.8 0.0 (index, r_ij, r_ik, theta, f_i1, f_i2, f_j1, f_j2, f_k1, f_k2, e)
|
||||
2 2.55 2.55 11.25 -621.539 -411.189 621.539 5035.95 411.189 -5035.95 0.0
|
||||
...
|
||||
1872 3.65 3.65 176.25 -0.00215132 -0.00412886 0.00215137 0.00111754 0.00412895 -0.00111757 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 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. It should be the same than the parameter
|
||||
"N" of the ".3b" file, otherwise its value is overwritten. "N"
|
||||
determines the number of table entries "M" that follow: "M = N * N *
|
||||
(N+1)" (symmetric triplet, e.g. SiCC) or "M = 2 * N * N * N" (asymmetric
|
||||
triplet, e.g. SiCSi). Therefore "M = 12 * 12 * 13 = 1872" in the above
|
||||
symmetric example. The parameters "rmin" and "rmax" are also required
|
||||
and determine the minimum and maximum of the inter-particle distances
|
||||
:math:`r_{ij}` and :math:`r_{ik}`.
|
||||
|
||||
Following a blank line, the next M lines of the angular table file list
|
||||
the tabulated values. On each line, the first value is the index from 1
|
||||
to M, the second value is the distance :math:`r_{ij}`, the third value
|
||||
is the distance :math:`r_{ik}`, the fourth value is the angle
|
||||
:math:`\theta_{ijk})`, the next six values are the force constants
|
||||
:math:`f_{i1}`, :math:`f_{i2}`, :math:`f_{j1}`, :math:`f_{j2}`,
|
||||
:math:`f_{k1}`, and :math:`f_{k2}`, and the last value is the energy e.
|
||||
|
||||
Note that one three-body potential file can contain many sections, each
|
||||
with a tabulated potential. LAMMPS reads the file section by section
|
||||
until it finds one that matches the specified *keyword* of appropriate
|
||||
section of the ".3b" file.
|
||||
|
||||
At the moment, only the *style* *linear* is allowed and
|
||||
implemented. After reading in the force table, it is internally stored
|
||||
in LAMMPS as a lookup table. For each triplet configuration occurring in
|
||||
the simulation within the cutoff distance, the next nearest tabulated
|
||||
triplet configuration is looked up. No interpolation is done. This
|
||||
allows for a very efficient force calculation with the stored force
|
||||
constants and energies. Due to the know table structure, the lookup can
|
||||
be done efficiently. It has been tested (:ref:`Scherer2 <Scherer2>`)
|
||||
that with a reasonably small bin size, the accuracy and speed is
|
||||
comparable to that of a Stillinger-Weber potential with tabulated
|
||||
three-body interactions (:doc:`pair_style sw/angle/table
|
||||
<pair_sw_angle_table>`) while the table format of this pair style allows
|
||||
for more flexible three-body interactions.
|
||||
|
||||
As for the Stillinger-Weber potential, the three-body potential file
|
||||
must contain entries for all the elements listed in the pair_coeff
|
||||
command. It can also contain entries for additional elements not being
|
||||
used in a particular simulation; LAMMPS ignores those entries.
|
||||
|
||||
For a single-element simulation, only a single entry is required
|
||||
(e.g. SiSiSi). For a two-element simulation, the file must contain 8
|
||||
entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that
|
||||
specify threebody parameters for all permutations of the two elements
|
||||
interacting in three-body configurations. Thus for 3 elements, 27
|
||||
entries would be required, etc.
|
||||
|
||||
As annotated above, the first element in the entry is the center atom in
|
||||
a three-body interaction. Thus an entry for SiCC means a Si atom with 2
|
||||
C atoms as neighbors. The tabulated three-body forces can in principle
|
||||
be specific to the three elements of the configuration. However, the
|
||||
user must ensure that it makes physically sense. E.g., the tabulated
|
||||
three-body forces for the entries CSiC and CCSi should be the same
|
||||
exchanging :math:`r_{ij}` with r_{ik}, :math:`f_{j1}` with
|
||||
:math:`f_{k1}`, and :math:`f_{j2}` with :math:`f_{k2}`.
|
||||
|
||||
Additional input files and reference data can be found at:
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-tutorials/ml
|
||||
|
||||
----------
|
||||
|
||||
Mixing, shift, table, tail correction, restart, rRESPA info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
As all interactions are tabulated, no mixing is performed.
|
||||
|
||||
This pair style does not support the :doc:`pair_modify <pair_modify>`
|
||||
shift, table, and tail options.
|
||||
|
||||
This pair style does not write its information to :doc:`binary restart
|
||||
files <restart>`, since it is stored in potential files. Thus, you need
|
||||
to re-specify the pair_style and pair_coeff commands in an input script
|
||||
that reads a restart file.
|
||||
|
||||
This pair style can only be used via the *pair* keyword of the
|
||||
:doc:`run_style respa <run_style>` command. It does not support the
|
||||
*inner*, *middle*, *outer* keywords.
|
||||
|
||||
----------
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This pair style is part of the MANYBODY package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` page for more info.
|
||||
|
||||
This pair style requires the :doc:`newton <newton>` setting to be "on"
|
||||
for pair interactions.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`, :doc:`pair sw/angle/table <pair_sw_angle_table>`
|
||||
|
||||
|
||||
----------
|
||||
|
||||
.. _Scherer2:
|
||||
|
||||
**(Scherer2)** C. Scherer, R. Scheid, D. Andrienko, and T. Bereau, J. Chem. Theor. Comp. 16, 3194-3204 (2020).
|
||||
|
||||
@ -94,6 +94,7 @@ amu
|
||||
Amzallag
|
||||
analytical
|
||||
Anders
|
||||
Andrienko
|
||||
Andzelm
|
||||
Ang
|
||||
anglegrad
|
||||
@ -239,6 +240,7 @@ benchmarking
|
||||
Bennet
|
||||
Berardi
|
||||
Beraun
|
||||
Bereau
|
||||
berendsen
|
||||
Berendsen
|
||||
berger
|
||||
@ -560,6 +562,7 @@ Crozier
|
||||
Cryst
|
||||
Crystallogr
|
||||
Csanyi
|
||||
csg
|
||||
csh
|
||||
cshrc
|
||||
CSiC
|
||||
@ -873,6 +876,7 @@ Eindhoven
|
||||
Eisenforschung
|
||||
Ejtehadi
|
||||
El
|
||||
el
|
||||
elaplong
|
||||
elastance
|
||||
Electroneg
|
||||
@ -1081,6 +1085,7 @@ fm
|
||||
fmackay
|
||||
fmag
|
||||
fmass
|
||||
fmatch
|
||||
fmm
|
||||
fmt
|
||||
fmtlib
|
||||
@ -2406,6 +2411,7 @@ Nstep
|
||||
Nsteplast
|
||||
Nstop
|
||||
nsub
|
||||
Nsw
|
||||
Nswap
|
||||
nt
|
||||
Nt
|
||||
@ -3036,6 +3042,7 @@ scalexz
|
||||
scaleyz
|
||||
Scalfi
|
||||
Schaik
|
||||
Scheid
|
||||
Schilfgarde
|
||||
Schimansky
|
||||
Schiotz
|
||||
@ -3213,6 +3220,7 @@ statvolt
|
||||
stdin
|
||||
stdio
|
||||
stdlib
|
||||
stdout
|
||||
steelblue
|
||||
Stegailov
|
||||
Steinbach
|
||||
|
||||
1875
examples/PACKAGES/manybody_table/1-1-1.table
Normal file
1875
examples/PACKAGES/manybody_table/1-1-1.table
Normal file
File diff suppressed because it is too large
Load Diff
3459
examples/PACKAGES/manybody_table/1-1-2.table
Normal file
3459
examples/PACKAGES/manybody_table/1-1-2.table
Normal file
File diff suppressed because it is too large
Load Diff
56
examples/PACKAGES/manybody_table/README
Normal file
56
examples/PACKAGES/manybody_table/README
Normal file
@ -0,0 +1,56 @@
|
||||
Three examples inputs for pair styles threebody/table (in.spce and
|
||||
in.spce2) and sw/angle/table (in.spce_sw). All inputs are for the
|
||||
simulation of coarse-grained SPC/E water.
|
||||
|
||||
A water molecule is represented by one coarse-grained (CG) bead.
|
||||
|
||||
For the two threebody/table examples the three-body (force) tables are
|
||||
in the files 1-1-1.table and 1-1-2.table. These have been parametrized
|
||||
with the kernel-based machine learning (ML) with the VOTCA package
|
||||
(https://gitlab.mpcdf.mpg.de/votca/votca). For a example on the
|
||||
parametrization, have a look at
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-examples/guide
|
||||
and
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-examples/ml.
|
||||
|
||||
In both cases, the parametrization is done according to a sample system,
|
||||
using the three-body forces of a Stillinger-Weber potential with
|
||||
tabulated angular forces (sw/angle/table) (see in.spce_sw). These then
|
||||
are learned with the covariant meshing technique with the settings files
|
||||
used in
|
||||
https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-examples/ml/3body/with_binning.
|
||||
|
||||
The first example, in.spce uses the LAMMPS data file (data.spce) with
|
||||
the starting configuration of 1000 CG water molecules, and a threebody
|
||||
file (spce.3b) which loads the 1-1-1.table file.
|
||||
|
||||
A hybrid/overlay pair style is used to sum a tabulated pairwise
|
||||
interaction (table_CG_CG.txt) with the tabulated threebody interactions.
|
||||
The tabulated pair interaction is the effectively the same as in the
|
||||
what is used by the in.spce_sw input using sw/angle/table pair style.
|
||||
|
||||
IMPORTANT NOTE: The threebody tables are parameterized without storing
|
||||
energies (the last column of the threebody table files is zero). This
|
||||
is due to a current limitation of the paramerization procedure.
|
||||
|
||||
This in.spce2 example uses the data.spce2 file and the threebody input
|
||||
spce2.3b. This is essentially the same simulation as in in.spce only
|
||||
the atom type of the first 100 CG water molecules has been changed from
|
||||
1 to 2. This is done to demonstrate how to run a simulation with
|
||||
different atom types.
|
||||
|
||||
For this (artificial) two-element simulation, the threebody file now
|
||||
contain 8 entries for: type1 type1 type1, type1 type1 type2, type1 type2
|
||||
type1, type1 type2 type2, type2 type1 type1, type2 type1 type2, type2
|
||||
type2 type1, type2 type2 type2. Each entry has the same structure as
|
||||
above. However, entries where the second and the third element are
|
||||
different require a different force table (1-1-2.table) instead of
|
||||
(1-1-1.table). 1-1-2.table contains exactly the force constants as
|
||||
1-1-1.table. However it has to have the asymmetric structure where both
|
||||
interparticle distances (r_ij and r_ik) are varied from rmin to rmax and
|
||||
therefore contains "M = 2 * N * N * N" (2 * 12 * 12 * 12 = 3456)
|
||||
entries.
|
||||
|
||||
The third example, in.spce_sw, contains the analytical twobody interactions
|
||||
and replaces the threebody term of the Stillinger-Weber potential with an
|
||||
angle/table style potential which is stored in the table_CG_CG_CG.txt file.
|
||||
1015
examples/PACKAGES/manybody_table/data.spce
Normal file
1015
examples/PACKAGES/manybody_table/data.spce
Normal file
File diff suppressed because it is too large
Load Diff
1016
examples/PACKAGES/manybody_table/data.spce2
Normal file
1016
examples/PACKAGES/manybody_table/data.spce2
Normal file
File diff suppressed because it is too large
Load Diff
21
examples/PACKAGES/manybody_table/in.spce
Normal file
21
examples/PACKAGES/manybody_table/in.spce
Normal file
@ -0,0 +1,21 @@
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with one atom type
|
||||
read_data data.spce
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce.3b type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
22
examples/PACKAGES/manybody_table/in.spce2
Normal file
22
examples/PACKAGES/manybody_table/in.spce2
Normal file
@ -0,0 +1,22 @@
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with two atom types
|
||||
read_data data.spce2
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff * * table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce2.3b type1 type2
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce2 id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
|
||||
20
examples/PACKAGES/manybody_table/in.spce_sw
Normal file
20
examples/PACKAGES/manybody_table/in.spce_sw
Normal file
@ -0,0 +1,20 @@
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
read_data data.spce
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 sw/angle/table
|
||||
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * sw/angle/table spce.sw type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform mom yes
|
||||
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 10 spce_sw.dump id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
88
examples/PACKAGES/manybody_table/log.1Jun22.spce.g++.1
Normal file
88
examples/PACKAGES/manybody_table/log.1Jun22.spce.g++.1
Normal file
@ -0,0 +1,88 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with one atom type
|
||||
read_data data.spce
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce.3b type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair threebody/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.487 | 5.487 | 5.487 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -5377.8719 0 -4484.5232 -320.10184
|
||||
100 296.01121 -5418.505 0 -4537.0342 -223.39972
|
||||
200 295.27654 -5430.3033 0 -4551.0202 794.29126
|
||||
300 302.16526 -5445.8048 0 -4546.0083 -11.568299
|
||||
400 308.59003 -5434.7181 0 -4515.7896 1.7337645
|
||||
500 295.346 -5436.0896 0 -4556.5996 778.73307
|
||||
600 293.14671 -5422.6082 0 -4549.6673 -148.64256
|
||||
700 307.63238 -5465.187 0 -4549.1103 285.18556
|
||||
800 313.16537 -5466.4124 0 -4533.8594 489.99301
|
||||
900 303.42954 -5506.3208 0 -4602.7595 360.05608
|
||||
1000 299.50926 -5446.8981 0 -4555.0107 993.95615
|
||||
Loop time of 5.15787 on 1 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 33.502 ns/day, 0.716 hours/ns, 193.879 timesteps/s
|
||||
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 4.8428 | 4.8428 | 4.8428 | 0.0 | 93.89
|
||||
Neigh | 0.27527 | 0.27527 | 0.27527 | 0.0 | 5.34
|
||||
Comm | 0.020461 | 0.020461 | 0.020461 | 0.0 | 0.40
|
||||
Output | 0.00020949 | 0.00020949 | 0.00020949 | 0.0 | 0.00
|
||||
Modify | 0.01198 | 0.01198 | 0.01198 | 0.0 | 0.23
|
||||
Other | | 0.007163 | | | 0.14
|
||||
|
||||
Nlocal: 1000 ave 1000 max 1000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 5900 ave 5900 max 5900 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 191126 ave 191126 max 191126 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 382252 ave 382252 max 382252 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 382252
|
||||
Ave neighs/atom = 382.252
|
||||
Neighbor list builds = 27
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:05
|
||||
88
examples/PACKAGES/manybody_table/log.1Jun22.spce.g++.4
Normal file
88
examples/PACKAGES/manybody_table/log.1Jun22.spce.g++.4
Normal file
@ -0,0 +1,88 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with one atom type
|
||||
read_data data.spce
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce.3b type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair threebody/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.87 | 3.87 | 3.87 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -5377.8719 0 -4484.5232 -320.10184
|
||||
100 296.01121 -5418.505 0 -4537.0342 -223.39972
|
||||
200 295.27654 -5430.3033 0 -4551.0202 794.29126
|
||||
300 302.16526 -5445.8048 0 -4546.0083 -11.568299
|
||||
400 308.59003 -5434.7181 0 -4515.7896 1.7337642
|
||||
500 295.346 -5436.0896 0 -4556.5996 778.73304
|
||||
600 293.14648 -5422.6082 0 -4549.6681 -148.67071
|
||||
700 307.6975 -5465.3018 0 -4549.0312 287.70203
|
||||
800 314.09436 -5467.6073 0 -4532.2879 522.73489
|
||||
900 300.85843 -5503.7551 0 -4607.85 491.78041
|
||||
1000 302.84638 -5468.3331 0 -4566.5083 338.05123
|
||||
Loop time of 1.60997 on 4 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 107.331 ns/day, 0.224 hours/ns, 621.131 timesteps/s
|
||||
94.2% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.28 | 1.2942 | 1.3018 | 0.8 | 80.38
|
||||
Neigh | 0.069763 | 0.070222 | 0.070788 | 0.2 | 4.36
|
||||
Comm | 0.19379 | 0.20979 | 0.23226 | 3.6 | 13.03
|
||||
Output | 0.0001711 | 0.00046947 | 0.0013621 | 0.0 | 0.03
|
||||
Modify | 0.021613 | 0.030016 | 0.038718 | 4.8 | 1.86
|
||||
Other | | 0.005309 | | | 0.33
|
||||
|
||||
Nlocal: 250 ave 257 max 240 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
Nghost: 3488.75 ave 3504 max 3478 min
|
||||
Histogram: 1 1 0 0 0 1 0 0 0 1
|
||||
Neighs: 47828 ave 49169 max 45782 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
FullNghs: 95656 ave 98253 max 91425 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
|
||||
Total # of neighbors = 382624
|
||||
Ave neighs/atom = 382.624
|
||||
Neighbor list builds = 27
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:01
|
||||
89
examples/PACKAGES/manybody_table/log.1Jun22.spce2.g++.1
Normal file
89
examples/PACKAGES/manybody_table/log.1Jun22.spce2.g++.1
Normal file
@ -0,0 +1,89 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with two atom types
|
||||
read_data data.spce2
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff * * table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce2.3b type1 type2
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce2 id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair threebody/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.487 | 5.487 | 5.487 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -5377.8719 0 -4484.5232 -320.10184
|
||||
100 296.01121 -5418.505 0 -4537.0342 -223.39972
|
||||
200 295.27654 -5430.3033 0 -4551.0202 794.29126
|
||||
300 302.16526 -5445.8048 0 -4546.0083 -11.568299
|
||||
400 308.59003 -5434.7181 0 -4515.7896 1.7337645
|
||||
500 295.346 -5436.0896 0 -4556.5996 778.73307
|
||||
600 293.14671 -5422.6082 0 -4549.6673 -148.64256
|
||||
700 307.63238 -5465.187 0 -4549.1103 285.18556
|
||||
800 313.16537 -5466.4124 0 -4533.8594 489.99301
|
||||
900 303.42954 -5506.3208 0 -4602.7595 360.05608
|
||||
1000 299.50926 -5446.8981 0 -4555.0107 993.95615
|
||||
Loop time of 5.16365 on 1 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 33.465 ns/day, 0.717 hours/ns, 193.661 timesteps/s
|
||||
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 4.8443 | 4.8443 | 4.8443 | 0.0 | 93.81
|
||||
Neigh | 0.27931 | 0.27931 | 0.27931 | 0.0 | 5.41
|
||||
Comm | 0.020302 | 0.020302 | 0.020302 | 0.0 | 0.39
|
||||
Output | 0.00022712 | 0.00022712 | 0.00022712 | 0.0 | 0.00
|
||||
Modify | 0.011944 | 0.011944 | 0.011944 | 0.0 | 0.23
|
||||
Other | | 0.007616 | | | 0.15
|
||||
|
||||
Nlocal: 1000 ave 1000 max 1000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 5900 ave 5900 max 5900 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 191126 ave 191126 max 191126 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 382252 ave 382252 max 382252 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 382252
|
||||
Ave neighs/atom = 382.252
|
||||
Neighbor list builds = 27
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:05
|
||||
89
examples/PACKAGES/manybody_table/log.1Jun22.spce2.g++.4
Normal file
89
examples/PACKAGES/manybody_table/log.1Jun22.spce2.g++.4
Normal file
@ -0,0 +1,89 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
# data file with two atom types
|
||||
read_data data.spce2
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.004 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 threebody/table
|
||||
|
||||
#pair coefficients
|
||||
pair_coeff * * table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * threebody/table spce2.3b type1 type2
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform loop geom mom yes
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 100 dump.spce2 id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair threebody/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.87 | 3.87 | 3.87 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -5377.8719 0 -4484.5232 -320.10184
|
||||
100 296.01121 -5418.505 0 -4537.0342 -223.39972
|
||||
200 295.27654 -5430.3033 0 -4551.0202 794.29126
|
||||
300 302.16526 -5445.8048 0 -4546.0083 -11.568299
|
||||
400 308.59003 -5434.7181 0 -4515.7896 1.7337642
|
||||
500 295.346 -5436.0896 0 -4556.5996 778.73304
|
||||
600 293.14648 -5422.6082 0 -4549.6681 -148.67071
|
||||
700 307.6975 -5465.3018 0 -4549.0312 287.70203
|
||||
800 314.09436 -5467.6073 0 -4532.2879 522.73489
|
||||
900 300.85843 -5503.7551 0 -4607.85 491.78041
|
||||
1000 302.84638 -5468.3331 0 -4566.5083 338.05123
|
||||
Loop time of 1.54853 on 4 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 111.590 ns/day, 0.215 hours/ns, 645.773 timesteps/s
|
||||
96.0% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.2599 | 1.2908 | 1.3259 | 2.1 | 83.36
|
||||
Neigh | 0.069097 | 0.071294 | 0.075502 | 0.9 | 4.60
|
||||
Comm | 0.12731 | 0.15884 | 0.19196 | 5.7 | 10.26
|
||||
Output | 0.00017674 | 0.0026016 | 0.0098653 | 8.2 | 0.17
|
||||
Modify | 0.0093453 | 0.011999 | 0.014575 | 2.3 | 0.77
|
||||
Other | | 0.01295 | | | 0.84
|
||||
|
||||
Nlocal: 250 ave 257 max 240 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
Nghost: 3488.75 ave 3504 max 3478 min
|
||||
Histogram: 1 1 0 0 0 1 0 0 0 1
|
||||
Neighs: 47828 ave 49169 max 45782 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
FullNghs: 95656 ave 98253 max 91425 min
|
||||
Histogram: 1 0 0 0 0 1 0 0 1 1
|
||||
|
||||
Total # of neighbors = 382624
|
||||
Ave neighs/atom = 382.624
|
||||
Neighbor list builds = 27
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:01
|
||||
87
examples/PACKAGES/manybody_table/log.1Jun22.spce_sw.g++.1
Normal file
87
examples/PACKAGES/manybody_table/log.1Jun22.spce_sw.g++.1
Normal file
@ -0,0 +1,87 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
read_data data.spce
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 sw/angle/table
|
||||
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * sw/angle/table spce.sw type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform mom yes
|
||||
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 10 spce_sw.dump id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair sw/angle/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.487 | 5.487 | 5.487 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -4572.9581 0 -3679.6093 -402.23914
|
||||
100 286.5642 -4508.6912 0 -3655.352 -610.63256
|
||||
200 291.59063 -4465.6368 0 -3597.3297 -218.54913
|
||||
300 298.40301 -4460.64 0 -3572.0468 302.96636
|
||||
400 305.99618 -4460.1128 0 -3548.9084 -68.022415
|
||||
500 301.94233 -4440.337 0 -3541.2043 179.36975
|
||||
600 308.95709 -4485.8412 0 -3565.8197 -95.917517
|
||||
700 291.69015 -4489.4465 0 -3620.843 -56.044939
|
||||
800 294.95653 -4496.904 0 -3618.5738 563.3456
|
||||
900 295.50533 -4478.1134 0 -3598.149 89.234288
|
||||
1000 308.63559 -4471.1612 0 -3552.0971 906.33706
|
||||
Loop time of 5.39753 on 1 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 32.015 ns/day, 0.750 hours/ns, 185.270 timesteps/s
|
||||
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 5.0954 | 5.0954 | 5.0954 | 0.0 | 94.40
|
||||
Neigh | 0.26201 | 0.26201 | 0.26201 | 0.0 | 4.85
|
||||
Comm | 0.020497 | 0.020497 | 0.020497 | 0.0 | 0.38
|
||||
Output | 0.00021869 | 0.00021869 | 0.00021869 | 0.0 | 0.00
|
||||
Modify | 0.011849 | 0.011849 | 0.011849 | 0.0 | 0.22
|
||||
Other | | 0.007577 | | | 0.14
|
||||
|
||||
Nlocal: 1000 ave 1000 max 1000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 5862 ave 5862 max 5862 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 191262 ave 191262 max 191262 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 382524 ave 382524 max 382524 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 382524
|
||||
Ave neighs/atom = 382.524
|
||||
Neighbor list builds = 26
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:05
|
||||
87
examples/PACKAGES/manybody_table/log.1Jun22.spce_sw.g++.4
Normal file
87
examples/PACKAGES/manybody_table/log.1Jun22.spce_sw.g++.4
Normal file
@ -0,0 +1,87 @@
|
||||
LAMMPS (4 May 2022)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units real
|
||||
atom_style atomic
|
||||
|
||||
read_data data.spce
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (31.0648 31.0648 31.0648)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
1000 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
|
||||
pair_style hybrid/overlay table linear 1200 sw/angle/table
|
||||
|
||||
pair_coeff 1 1 table table_CG_CG.txt VOTCA
|
||||
pair_coeff * * sw/angle/table spce.sw type
|
||||
|
||||
fix 1 all nvt temp 300.0 300.0 200.0
|
||||
|
||||
velocity all create 300 432567 dist uniform mom yes
|
||||
|
||||
timestep 2.0
|
||||
|
||||
thermo 100
|
||||
#dump 2 all custom 10 spce_sw.dump id type x y z fx fy fz
|
||||
|
||||
run 1000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 5 5 5
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair table, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair sw/angle/table, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.87 | 3.87 | 3.87 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 300 -4572.9581 0 -3679.6093 -402.23914
|
||||
100 286.5642 -4508.6912 0 -3655.352 -610.63256
|
||||
200 291.59063 -4465.6368 0 -3597.3297 -218.54913
|
||||
300 298.40301 -4460.64 0 -3572.0468 302.96636
|
||||
400 305.99618 -4460.1128 0 -3548.9084 -68.022415
|
||||
500 301.94233 -4440.337 0 -3541.2043 179.36975
|
||||
600 308.95709 -4485.8412 0 -3565.8197 -95.917517
|
||||
700 291.69015 -4489.4465 0 -3620.843 -56.044939
|
||||
800 294.95653 -4496.904 0 -3618.5738 563.3456
|
||||
900 295.50533 -4478.1134 0 -3598.149 89.234292
|
||||
1000 308.63559 -4471.1612 0 -3552.0971 906.33708
|
||||
Loop time of 1.57073 on 4 procs for 1000 steps with 1000 atoms
|
||||
|
||||
Performance: 110.012 ns/day, 0.218 hours/ns, 636.646 timesteps/s
|
||||
97.1% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.3064 | 1.3347 | 1.3706 | 2.0 | 84.97
|
||||
Neigh | 0.065344 | 0.070299 | 0.07737 | 1.7 | 4.48
|
||||
Comm | 0.117 | 0.15297 | 0.18685 | 6.5 | 9.74
|
||||
Output | 0.00016937 | 0.00055648 | 0.0017097 | 0.0 | 0.04
|
||||
Modify | 0.0073414 | 0.0079027 | 0.0085343 | 0.6 | 0.50
|
||||
Other | | 0.0043 | | | 0.27
|
||||
|
||||
Nlocal: 250 ave 254 max 247 min
|
||||
Histogram: 1 1 0 0 0 1 0 0 0 1
|
||||
Nghost: 3473.25 ave 3490 max 3450 min
|
||||
Histogram: 1 0 0 0 0 1 0 1 0 1
|
||||
Neighs: 47815.5 ave 48520 max 47134 min
|
||||
Histogram: 1 0 1 0 0 0 1 0 0 1
|
||||
FullNghs: 95631 ave 97203 max 94083 min
|
||||
Histogram: 1 0 1 0 0 0 1 0 0 1
|
||||
|
||||
Total # of neighbors = 382524
|
||||
Ave neighs/atom = 382.524
|
||||
Neighbor list builds = 26
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:01
|
||||
8
examples/PACKAGES/manybody_table/spce.3b
Normal file
8
examples/PACKAGES/manybody_table/spce.3b
Normal file
@ -0,0 +1,8 @@
|
||||
type
|
||||
type
|
||||
type
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
18
examples/PACKAGES/manybody_table/spce.sw
Normal file
18
examples/PACKAGES/manybody_table/spce.sw
Normal file
@ -0,0 +1,18 @@
|
||||
type
|
||||
type
|
||||
type
|
||||
1 #epsilon in kcal/mol
|
||||
1 #sigma in dimensionless
|
||||
3.7 # a in Ang
|
||||
1.0 #lambda dimensionless
|
||||
0.8 #gamma in Ang
|
||||
0.0 #costheta0 dimensionless
|
||||
0 #two body part A=0
|
||||
0 #two body part B=0
|
||||
0 #two body part p=0
|
||||
0 #two body part q=0
|
||||
0.0 # use the standard Stillinger-Weber cutoff
|
||||
table_CG_CG_CG.txt
|
||||
VOTCA
|
||||
linear
|
||||
1001
|
||||
71
examples/PACKAGES/manybody_table/spce2.3b
Normal file
71
examples/PACKAGES/manybody_table/spce2.3b
Normal file
@ -0,0 +1,71 @@
|
||||
type1
|
||||
type1
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type1
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type2
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type2
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type1
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type1
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type2
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type2
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
1203
examples/PACKAGES/manybody_table/table_CG_CG.txt
Normal file
1203
examples/PACKAGES/manybody_table/table_CG_CG.txt
Normal file
File diff suppressed because it is too large
Load Diff
1004
examples/PACKAGES/manybody_table/table_CG_CG_CG.txt
Normal file
1004
examples/PACKAGES/manybody_table/table_CG_CG_CG.txt
Normal file
File diff suppressed because it is too large
Load Diff
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -997,6 +997,8 @@
|
||||
/neb.h
|
||||
/netcdf_units.cpp
|
||||
/netcdf_units.h
|
||||
/pair_3b_table.cpp
|
||||
/pair_3b_table.h
|
||||
/pair_adp.cpp
|
||||
/pair_adp.h
|
||||
/pair_agni.cpp
|
||||
@ -1291,6 +1293,8 @@
|
||||
/pair_sph_taitwater_morris.h
|
||||
/pair_sw.cpp
|
||||
/pair_sw.h
|
||||
/pair_sw_3b_table.cpp
|
||||
/pair_sw_3b_table.h
|
||||
/pair_sw_mod.cpp
|
||||
/pair_sw_mod.h
|
||||
/pair_tersoff.cpp
|
||||
|
||||
@ -55,7 +55,7 @@ class PairSW : public Pair {
|
||||
|
||||
void settings(int, char **) override;
|
||||
virtual void allocate();
|
||||
void read_file(char *);
|
||||
virtual void read_file(char *);
|
||||
virtual void setup_params();
|
||||
void twobody(Param *, double, double &, int, double &);
|
||||
virtual void threebody(Param *, Param *, Param *, double, double, double *, double *, double *,
|
||||
|
||||
753
src/MANYBODY/pair_sw_angle_table.cpp
Normal file
753
src/MANYBODY/pair_sw_angle_table.cpp
Normal file
@ -0,0 +1,753 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Christoph Scherer (MPIP Mainz)
|
||||
scherer@mpip-mainz.mpg.de
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_sw_angle_table.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "table_file_reader.h"
|
||||
#include "potential_file_reader.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
using MathConst::DEG2RAD;
|
||||
using MathConst::MY_PI;
|
||||
using MathConst::RAD2DEG;
|
||||
|
||||
#define DELTA 4
|
||||
|
||||
enum { LINEAR, SPLINE };
|
||||
|
||||
static constexpr double TINY = 1.0e-10;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSWAngleTable::PairSWAngleTable(LAMMPS *lmp) : PairSW(lmp), table_params(nullptr)
|
||||
{
|
||||
unit_convert_flag = utils::NOCONVERT;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
check if allocated, since class can be destructed when incomplete
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
PairSWAngleTable::~PairSWAngleTable()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
for (int m = 0; m < nparams; m++) free_param(&table_params[m]); // free_param will call free_table
|
||||
memory->destroy(params);
|
||||
memory->destroy(table_params);
|
||||
memory->destroy(elem3param);
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(neighshort);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
|
||||
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
|
||||
tagint itag,jtag;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
|
||||
double rsq,rsq1,rsq2;
|
||||
double delr1[3],delr2[3],fj[3],fk[3];
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
tagint *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
double fxtmp,fytmp,fztmp;
|
||||
|
||||
// loop over full neighbor list of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itag = tag[i];
|
||||
itype = map[type[i]];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
// two-body interactions, skip half of them
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
int numshort = 0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
|
||||
jtype = map[type[j]];
|
||||
ijparam = elem3param[itype][jtype][jtype];
|
||||
if (rsq >= params[ijparam].cutsq) {
|
||||
continue;
|
||||
} else {
|
||||
neighshort[numshort++] = j;
|
||||
if (numshort >= maxshort) {
|
||||
maxshort += maxshort/2;
|
||||
memory->grow(neighshort,maxshort,"pair:neighshort");
|
||||
}
|
||||
}
|
||||
|
||||
jtag = tag[j];
|
||||
if (itag > jtag) {
|
||||
if ((itag+jtag) % 2 == 0) continue;
|
||||
} else if (itag < jtag) {
|
||||
if ((itag+jtag) % 2 == 1) continue;
|
||||
} else {
|
||||
if (x[j][2] < ztmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
|
||||
}
|
||||
|
||||
twobody(¶ms[ijparam],rsq,fpair,eflag,evdwl);
|
||||
|
||||
fxtmp += delx*fpair;
|
||||
fytmp += dely*fpair;
|
||||
fztmp += delz*fpair;
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,0.0,fpair,delx,dely,delz);
|
||||
}
|
||||
|
||||
jnumm1 = numshort - 1;
|
||||
|
||||
for (jj = 0; jj < jnumm1; jj++) {
|
||||
j = neighshort[jj];
|
||||
jtype = map[type[j]];
|
||||
ijparam = elem3param[itype][jtype][jtype];
|
||||
delr1[0] = x[j][0] - xtmp;
|
||||
delr1[1] = x[j][1] - ytmp;
|
||||
delr1[2] = x[j][2] - ztmp;
|
||||
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
|
||||
|
||||
double fjxtmp,fjytmp,fjztmp;
|
||||
fjxtmp = fjytmp = fjztmp = 0.0;
|
||||
|
||||
for (kk = jj+1; kk < numshort; kk++) {
|
||||
k = neighshort[kk];
|
||||
ktype = map[type[k]];
|
||||
ikparam = elem3param[itype][ktype][ktype];
|
||||
ijkparam = elem3param[itype][jtype][ktype];
|
||||
|
||||
delr2[0] = x[k][0] - xtmp;
|
||||
delr2[1] = x[k][1] - ytmp;
|
||||
delr2[2] = x[k][2] - ztmp;
|
||||
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
|
||||
|
||||
threebody_table(¶ms[ijparam],¶ms[ikparam],&table_params[ijkparam],
|
||||
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
|
||||
|
||||
fxtmp -= fj[0] + fk[0];
|
||||
fytmp -= fj[1] + fk[1];
|
||||
fztmp -= fj[2] + fk[2];
|
||||
fjxtmp += fj[0];
|
||||
fjytmp += fj[1];
|
||||
fjztmp += fj[2];
|
||||
f[k][0] += fk[0];
|
||||
f[k][1] += fk[1];
|
||||
f[k][2] += fk[2];
|
||||
|
||||
if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
|
||||
}
|
||||
f[j][0] += fjxtmp;
|
||||
f[j][1] += fjytmp;
|
||||
f[j][2] += fjztmp;
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::read_file(char *file)
|
||||
{
|
||||
if (params) {
|
||||
for (int m = 0; m < nparams; m++) free_param(&table_params[m]); // free_param will call free_table
|
||||
memory->destroy(params);
|
||||
memory->destroy(table_params);
|
||||
memory->destroy(elem3param);
|
||||
}
|
||||
|
||||
nparams = maxparam = 0;
|
||||
|
||||
// open file on proc 0
|
||||
|
||||
if (comm->me == 0) {
|
||||
PotentialFileReader reader(lmp, file, "sw", unit_convert_flag);
|
||||
char *line;
|
||||
|
||||
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
|
||||
std::string iname = values.next_string();
|
||||
std::string jname = values.next_string();
|
||||
std::string kname = values.next_string();
|
||||
|
||||
// ielement,jelement,kelement = 1st args
|
||||
// if all 3 args are in element list, then parse this line
|
||||
// else skip to next entry in file
|
||||
int ielement, jelement, kelement;
|
||||
|
||||
for (ielement = 0; ielement < nelements; ielement++)
|
||||
if (iname == elements[ielement]) break;
|
||||
if (ielement == nelements) continue;
|
||||
for (jelement = 0; jelement < nelements; jelement++)
|
||||
if (jname == elements[jelement]) break;
|
||||
if (jelement == nelements) continue;
|
||||
for (kelement = 0; kelement < nelements; kelement++)
|
||||
if (kname == elements[kelement]) break;
|
||||
if (kelement == nelements) continue;
|
||||
|
||||
// load up parameter settings and error check their values
|
||||
|
||||
if (nparams == maxparam) {
|
||||
maxparam += DELTA;
|
||||
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
|
||||
"pair:params");
|
||||
table_params = (ParamTable *) memory->srealloc(table_params,maxparam*sizeof(ParamTable),
|
||||
"pair:table_params");
|
||||
|
||||
// make certain all addional allocated storage is initialized
|
||||
// to avoid false positives when checking with valgrind
|
||||
|
||||
memset(params + nparams, 0, DELTA*sizeof(Param));
|
||||
}
|
||||
|
||||
params[nparams].ielement = ielement;
|
||||
params[nparams].jelement = jelement;
|
||||
params[nparams].kelement = kelement;
|
||||
params[nparams].epsilon = values.next_double();
|
||||
params[nparams].sigma = values.next_double();
|
||||
params[nparams].littlea = values.next_double();
|
||||
params[nparams].lambda = values.next_double();
|
||||
params[nparams].gamma = values.next_double();
|
||||
params[nparams].costheta = values.next_double();
|
||||
params[nparams].biga = values.next_double();
|
||||
params[nparams].bigb = values.next_double();
|
||||
params[nparams].powerp = values.next_double();
|
||||
params[nparams].powerq = values.next_double();
|
||||
params[nparams].tol = values.next_double();
|
||||
|
||||
// read parameters of angle table
|
||||
std::string name = values.next_string();
|
||||
table_params[nparams].tablenamelength = name.length()+1;
|
||||
table_params[nparams].tablename = utils::strdup(name);
|
||||
|
||||
name = values.next_string();
|
||||
table_params[nparams].keywordlength = name.length()+1;
|
||||
table_params[nparams].keyword = utils::strdup(name);
|
||||
|
||||
name = values.next_string();
|
||||
if (name == "linear") table_params[nparams].tabstyle = LINEAR;
|
||||
else if (name == "spline") table_params[nparams].tabstyle = SPLINE;
|
||||
else error->all(FLERR,"Unknown table style {} of angle table file", name);
|
||||
table_params[nparams].tablength = values.next_int();
|
||||
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
if (params[nparams].epsilon < 0.0 || params[nparams].sigma < 0.0 ||
|
||||
params[nparams].littlea < 0.0 || params[nparams].lambda < 0.0 ||
|
||||
params[nparams].gamma < 0.0 || params[nparams].biga < 0.0 ||
|
||||
params[nparams].bigb < 0.0 || params[nparams].powerp < 0.0 ||
|
||||
params[nparams].powerq < 0.0 || params[nparams].tol < 0.0 ||
|
||||
table_params[nparams].tabstyle < 0.0 || table_params[nparams].tablength < 0.0)
|
||||
error->one(FLERR,"Illegal Stillinger-Weber parameter");
|
||||
|
||||
nparams++;
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
|
||||
|
||||
if (comm->me != 0) {
|
||||
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
|
||||
table_params = (ParamTable *) memory->srealloc(table_params,maxparam*sizeof(ParamTable), "pair:table_params");
|
||||
}
|
||||
|
||||
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
|
||||
MPI_Bcast(table_params, maxparam*sizeof(ParamTable), MPI_BYTE, 0, world);
|
||||
|
||||
// for each set of parameters, broadcast table name and keyword and read angle table
|
||||
for (int m = 0; m < nparams; ++m){
|
||||
if (comm->me != 0) {
|
||||
table_params[m].tablename = new char[table_params[m].tablenamelength];
|
||||
table_params[m].keyword = new char[table_params[m].keywordlength];
|
||||
}
|
||||
MPI_Bcast(table_params[m].tablename, table_params[m].tablenamelength, MPI_CHAR, 0, world);
|
||||
MPI_Bcast(table_params[m].keyword, table_params[m].keywordlength, MPI_CHAR, 0, world);
|
||||
|
||||
// initialize angtable
|
||||
memory->create(table_params[m].angtable,1,"table_params:angtable");
|
||||
null_table(table_params[m].angtable);
|
||||
|
||||
// call read_table to read corresponding tabulated angle file (only called by process 0)
|
||||
if (comm->me == 0) read_table(table_params[m].angtable,table_params[m].tablename,table_params[m].keyword);
|
||||
|
||||
// broadcast read in angtable to all processes
|
||||
bcast_table(table_params[m].angtable);
|
||||
|
||||
// the following table manipulations are done in all processes
|
||||
// error check on table parameters
|
||||
if (table_params[m].angtable->ninput <= 1) error->one(FLERR,"Invalid angle table length");
|
||||
|
||||
// error check on parameter range of angle table
|
||||
double alo,ahi;
|
||||
alo = table_params[m].angtable->afile[0];
|
||||
ahi = table_params[m].angtable->afile[table_params[m].angtable->ninput-1];
|
||||
if (fabs(alo-0.0) > TINY || fabs(ahi-180.0) > TINY)
|
||||
error->all(FLERR,"Angle table must range from 0 to 180 degrees");
|
||||
|
||||
// convert theta from degrees to radians
|
||||
for (int i = 0; i < table_params[m].angtable->ninput; ++i){
|
||||
table_params[m].angtable->afile[i] *= MY_PI/180.0;
|
||||
table_params[m].angtable->ffile[i] *= 180.0/MY_PI;
|
||||
}
|
||||
|
||||
// spline read-in table and compute a,e,f vectors within table
|
||||
spline_table(table_params[m].angtable);
|
||||
// compute_table needs parameter params[m].length for this specific interaction as
|
||||
// read in value length from .sw file can be different from value in angle table file
|
||||
compute_table(table_params[m].angtable,table_params[m].tablength);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::threebody_table(Param *paramij, Param *paramik, ParamTable *table_paramijk,
|
||||
double rsq1, double rsq2, double *delr1, double *delr2,
|
||||
double *fj, double *fk, int eflag, double &eng)
|
||||
{
|
||||
double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
|
||||
double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
|
||||
double rinv12,cs,facexp;
|
||||
double ftheta,facradtable,frad1table,frad2table,var;
|
||||
double acosprime,gradj1,gradj2,gradk1,gradk2,fprimetheta;
|
||||
|
||||
r1 = sqrt(rsq1);
|
||||
rinvsq1 = 1.0/rsq1;
|
||||
rainv1 = 1.0/(r1 - paramij->cut);
|
||||
gsrainv1 = paramij->sigma_gamma * rainv1;
|
||||
gsrainvsq1 = gsrainv1*rainv1/r1;
|
||||
expgsrainv1 = exp(gsrainv1);
|
||||
|
||||
r2 = sqrt(rsq2);
|
||||
rinvsq2 = 1.0/rsq2;
|
||||
rainv2 = 1.0/(r2 - paramik->cut);
|
||||
gsrainv2 = paramik->sigma_gamma * rainv2;
|
||||
gsrainvsq2 = gsrainv2*rainv2/r2;
|
||||
expgsrainv2 = exp(gsrainv2);
|
||||
|
||||
facexp = expgsrainv1*expgsrainv2;
|
||||
|
||||
rinv12 = 1.0/(r1*r2);
|
||||
cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
|
||||
|
||||
var = acos(cs);
|
||||
|
||||
// look up energy (f(theta), ftheta) and force (df(theta)/dtheta, fprimetheta) at
|
||||
// angle theta (var) in angle table belonging to parameter set paramijk
|
||||
uf_lookup(table_paramijk, var, ftheta, fprimetheta);
|
||||
|
||||
acosprime = 1.0 / (sqrt(1 - cs*cs ) );
|
||||
|
||||
facradtable = facexp*ftheta;
|
||||
frad1table = facradtable*gsrainvsq1;
|
||||
frad2table = facradtable*gsrainvsq2;
|
||||
gradj1 = acosprime * cs * rinvsq1 * facexp * fprimetheta;
|
||||
gradj2 = acosprime * rinv12 * facexp * fprimetheta;
|
||||
gradk1 = acosprime * cs * rinvsq2 * facexp * fprimetheta;
|
||||
gradk2 = acosprime * rinv12 * facexp * fprimetheta;
|
||||
|
||||
fj[0] = delr1[0]*(frad1table+gradj1)-delr2[0]*gradj2;
|
||||
fj[1] = delr1[1]*(frad1table+gradj1)-delr2[1]*gradj2;
|
||||
fj[2] = delr1[2]*(frad1table+gradj1)-delr2[2]*gradj2;
|
||||
|
||||
fk[0] = delr2[0]*(frad2table+gradk1)-delr1[0]*gradk2;
|
||||
fk[1] = delr2[1]*(frad2table+gradk1)-delr1[1]*gradk2;
|
||||
fk[2] = delr2[2]*(frad2table+gradk1)-delr1[2]*gradk2;
|
||||
|
||||
if (eflag) eng = facradtable;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
read table file, only called by proc 0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::read_table(Table *tb, char *file, char *keyword)
|
||||
{
|
||||
TableFileReader reader(lmp, file, "angletable");
|
||||
|
||||
char *line = reader.find_section_start(keyword);
|
||||
|
||||
if (!line) { error->one(FLERR, "Did not find keyword in table file"); }
|
||||
|
||||
// read args on 2nd line of section
|
||||
// allocate table arrays for file values
|
||||
|
||||
line = reader.next_line();
|
||||
param_extract(tb, line);
|
||||
memory->create(tb->afile, tb->ninput, "angle:afile");
|
||||
memory->create(tb->efile, tb->ninput, "angle:efile");
|
||||
memory->create(tb->ffile, tb->ninput, "angle:ffile");
|
||||
|
||||
// read a,e,f table values from file
|
||||
|
||||
int cerror = 0;
|
||||
reader.skip_line();
|
||||
for (int i = 0; i < tb->ninput; i++) {
|
||||
line = reader.next_line(4);
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
values.next_int();
|
||||
tb->afile[i] = values.next_double();
|
||||
tb->efile[i] = values.next_double();
|
||||
tb->ffile[i] = values.next_double();
|
||||
} catch (TokenizerException &) {
|
||||
++cerror;
|
||||
}
|
||||
}
|
||||
|
||||
// warn if data was read incompletely, e.g. columns were missing
|
||||
|
||||
if (cerror)
|
||||
error->warning(FLERR, "{} of {} lines in table incomplete or could not be parsed", cerror,
|
||||
tb->ninput);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
build spline representation of e,f over entire range of read-in table
|
||||
this function sets these values in e2file,f2file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::spline_table(Table *tb)
|
||||
{
|
||||
memory->create(tb->e2file, tb->ninput, "angle:e2file");
|
||||
memory->create(tb->f2file, tb->ninput, "angle:f2file");
|
||||
|
||||
double ep0 = -tb->ffile[0];
|
||||
double epn = -tb->ffile[tb->ninput - 1];
|
||||
spline(tb->afile, tb->efile, tb->ninput, ep0, epn, tb->e2file);
|
||||
|
||||
if (tb->fpflag == 0) {
|
||||
tb->fplo = (tb->ffile[1] - tb->ffile[0]) / (tb->afile[1] - tb->afile[0]);
|
||||
tb->fphi = (tb->ffile[tb->ninput - 1] - tb->ffile[tb->ninput - 2]) /
|
||||
(tb->afile[tb->ninput - 1] - tb->afile[tb->ninput - 2]);
|
||||
}
|
||||
|
||||
double fp0 = tb->fplo;
|
||||
double fpn = tb->fphi;
|
||||
spline(tb->afile, tb->ffile, tb->ninput, fp0, fpn, tb->f2file);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute a,e,f vectors from splined values
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::compute_table(Table *tb,int length)
|
||||
{
|
||||
// delta = table spacing in angle for N-1 bins
|
||||
|
||||
int tlm1 = length - 1;
|
||||
tb->delta = MY_PI / tlm1;
|
||||
tb->invdelta = 1.0 / tb->delta;
|
||||
tb->deltasq6 = tb->delta * tb->delta / 6.0;
|
||||
|
||||
// N-1 evenly spaced bins in angle from 0 to PI
|
||||
// ang,e,f = value at lower edge of bin
|
||||
// de,df values = delta values of e,f
|
||||
// ang,e,f are N in length so de,df arrays can compute difference
|
||||
|
||||
memory->create(tb->ang, length, "angle:ang");
|
||||
memory->create(tb->e, length, "angle:e");
|
||||
memory->create(tb->de, length, "angle:de");
|
||||
memory->create(tb->f, length, "angle:f");
|
||||
memory->create(tb->df, length, "angle:df");
|
||||
memory->create(tb->e2, length, "angle:e2");
|
||||
memory->create(tb->f2, length, "angle:f2");
|
||||
|
||||
double a;
|
||||
for (int i = 0; i < length; i++) {
|
||||
a = i * tb->delta;
|
||||
tb->ang[i] = a;
|
||||
tb->e[i] = splint(tb->afile, tb->efile, tb->e2file, tb->ninput, a);
|
||||
tb->f[i] = splint(tb->afile, tb->ffile, tb->f2file, tb->ninput, a);
|
||||
}
|
||||
|
||||
for (int i = 0; i < tlm1; i++) {
|
||||
tb->de[i] = tb->e[i + 1] - tb->e[i];
|
||||
tb->df[i] = tb->f[i + 1] - tb->f[i];
|
||||
}
|
||||
// get final elements from linear extrapolation
|
||||
tb->de[tlm1] = 2.0 * tb->de[tlm1 - 1] - tb->de[tlm1 - 2];
|
||||
tb->df[tlm1] = 2.0 * tb->df[tlm1 - 1] - tb->df[tlm1 - 2];
|
||||
|
||||
double ep0 = -tb->f[0];
|
||||
double epn = -tb->f[tlm1];
|
||||
spline(tb->ang, tb->e, length, ep0, epn, tb->e2);
|
||||
spline(tb->ang, tb->f, length, tb->fplo, tb->fphi, tb->f2);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
broadcast read-in table info from proc 0 to other procs
|
||||
this function communicates these values in Table:
|
||||
ninput,afile,efile,ffile,fpflag,fplo,fphi,theta0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::bcast_table(Table *tb)
|
||||
{
|
||||
MPI_Bcast(&tb->ninput, 1, MPI_INT, 0, world);
|
||||
|
||||
int me;
|
||||
MPI_Comm_rank(world, &me);
|
||||
if (me > 0) {
|
||||
memory->create(tb->afile, tb->ninput, "angle:afile");
|
||||
memory->create(tb->efile, tb->ninput, "angle:efile");
|
||||
memory->create(tb->ffile, tb->ninput, "angle:ffile");
|
||||
}
|
||||
|
||||
MPI_Bcast(tb->afile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->efile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->ffile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
|
||||
MPI_Bcast(&tb->fpflag, 1, MPI_INT, 0, world);
|
||||
if (tb->fpflag) {
|
||||
MPI_Bcast(&tb->fplo, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&tb->fphi, 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
MPI_Bcast(&tb->theta0, 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::null_table(Table *tb)
|
||||
{
|
||||
tb->afile = tb->efile = tb->ffile = nullptr;
|
||||
tb->e2file = tb->f2file = nullptr;
|
||||
tb->ang = tb->e = tb->de = nullptr;
|
||||
tb->f = tb->df = tb->e2 = tb->f2 = nullptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::free_table(Table *tb)
|
||||
{
|
||||
memory->destroy(tb->afile);
|
||||
memory->destroy(tb->efile);
|
||||
memory->destroy(tb->ffile);
|
||||
memory->destroy(tb->e2file);
|
||||
memory->destroy(tb->f2file);
|
||||
|
||||
memory->destroy(tb->ang);
|
||||
memory->destroy(tb->e);
|
||||
memory->destroy(tb->de);
|
||||
memory->destroy(tb->f);
|
||||
memory->destroy(tb->df);
|
||||
memory->destroy(tb->e2);
|
||||
memory->destroy(tb->f2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::free_param(ParamTable *pm)
|
||||
{
|
||||
// call free_table to destroy associated angle table
|
||||
free_table(pm->angtable);
|
||||
// then destroy associated angle table
|
||||
delete[] pm->keyword;
|
||||
delete[] pm->tablename;
|
||||
memory->sfree(pm->angtable);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract attributes from parameter line in table section
|
||||
format of line: N value FP fplo fphi EQ theta0
|
||||
N is required, other params are optional
|
||||
|
||||
only called by read_table, only called by proc 0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::param_extract(Table *tb, char *line)
|
||||
{
|
||||
tb->ninput = 0;
|
||||
tb->fpflag = 0;
|
||||
tb->theta0 = MY_PI;
|
||||
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
|
||||
while (values.has_next()) {
|
||||
std::string word = values.next_string();
|
||||
|
||||
if (word == "N") {
|
||||
tb->ninput = values.next_int();
|
||||
} else if (word == "FP") {
|
||||
tb->fpflag = 1;
|
||||
tb->fplo = values.next_double();
|
||||
tb->fphi = values.next_double();
|
||||
tb->fplo *= RAD2DEG * RAD2DEG;
|
||||
tb->fphi *= RAD2DEG * RAD2DEG;
|
||||
} else if (word == "EQ") {
|
||||
tb->theta0 = DEG2RAD * values.next_double();
|
||||
} else {
|
||||
error->one(FLERR, "Invalid keyword in angle table parameters");
|
||||
}
|
||||
}
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
if (tb->ninput == 0) error->one(FLERR, "Angle table parameters did not set N");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
spline and splint routines modified from Numerical Recipes
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
|
||||
{
|
||||
int i, k;
|
||||
double p, qn, sig, un;
|
||||
double *u = new double[n];
|
||||
|
||||
if (yp1 > 0.99e300)
|
||||
y2[0] = u[0] = 0.0;
|
||||
else {
|
||||
y2[0] = -0.5;
|
||||
u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
|
||||
}
|
||||
for (i = 1; i < n - 1; i++) {
|
||||
sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
|
||||
p = sig * y2[i - 1] + 2.0;
|
||||
y2[i] = (sig - 1.0) / p;
|
||||
u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
|
||||
u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
|
||||
}
|
||||
if (ypn > 0.99e300)
|
||||
qn = un = 0.0;
|
||||
else {
|
||||
qn = 0.5;
|
||||
un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
|
||||
}
|
||||
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
|
||||
for (k = n - 2; k >= 0; k--) y2[k] = y2[k] * y2[k + 1] + u[k];
|
||||
|
||||
delete[] u;
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairSWAngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
|
||||
{
|
||||
int klo, khi, k;
|
||||
double h, b, a, y;
|
||||
|
||||
klo = 0;
|
||||
khi = n - 1;
|
||||
while (khi - klo > 1) {
|
||||
k = (khi + klo) >> 1;
|
||||
if (xa[k] > x)
|
||||
khi = k;
|
||||
else
|
||||
klo = k;
|
||||
}
|
||||
h = xa[khi] - xa[klo];
|
||||
a = (xa[khi] - x) / h;
|
||||
b = (x - xa[klo]) / h;
|
||||
y = a * ya[klo] + b * ya[khi] +
|
||||
((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0;
|
||||
return y;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
calculate potential u and force f at angle x
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSWAngleTable::uf_lookup(ParamTable *pm, double x, double &u, double &f)
|
||||
{
|
||||
if (!std::isfinite(x)) { error->one(FLERR, "Illegal angle in angle style table"); }
|
||||
|
||||
double fraction,a,b;
|
||||
|
||||
// invdelta is based on tablength-1
|
||||
int itable = static_cast<int>(x * pm->angtable->invdelta);
|
||||
if (itable < 0) itable = 0;
|
||||
if (itable >= pm->tablength) itable = pm->tablength - 1;
|
||||
|
||||
if (pm->tabstyle == LINEAR) {
|
||||
fraction = (x - pm->angtable->ang[itable]) * pm->angtable->invdelta;
|
||||
u = pm->angtable->e[itable] + fraction*pm->angtable->de[itable];
|
||||
f = pm->angtable->f[itable] + fraction*pm->angtable->df[itable];
|
||||
} else if (pm->tabstyle == SPLINE) {
|
||||
fraction = (x - pm->angtable->ang[itable]) * pm->angtable->invdelta;
|
||||
|
||||
b = (x - pm->angtable->ang[itable]) * pm->angtable->invdelta;
|
||||
a = 1.0 - b;
|
||||
u = a * pm->angtable->e[itable] + b * pm->angtable->e[itable+1] +
|
||||
((a * a * a - a) * pm->angtable->e2[itable] + (b * b * b - b) * pm->angtable->e2[itable+1]) *
|
||||
pm->angtable->deltasq6;
|
||||
f = a * pm->angtable->f[itable] + b * pm->angtable->f[itable+1] +
|
||||
((a * a * a - a) * pm->angtable->f2[itable] + (b * b * b - b) * pm->angtable->f2[itable+1]) *
|
||||
pm->angtable->deltasq6;
|
||||
}
|
||||
}
|
||||
77
src/MANYBODY/pair_sw_angle_table.h
Normal file
77
src/MANYBODY/pair_sw_angle_table.h
Normal file
@ -0,0 +1,77 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
// clang-format off
|
||||
PairStyle(sw/angle/table,PairSWAngleTable);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SW_ANGLE_TABLE_H
|
||||
#define LMP_PAIR_SW_ANGLE_TABLE_H
|
||||
|
||||
#include "pair_sw.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSWAngleTable : public PairSW {
|
||||
public:
|
||||
PairSWAngleTable(class LAMMPS *);
|
||||
~PairSWAngleTable() override;
|
||||
void compute(int, int) override;
|
||||
|
||||
static constexpr int NPARAMS_PER_LINE = 18;
|
||||
|
||||
// use struct Table as in class AngleTable
|
||||
struct Table {
|
||||
int ninput, fpflag;
|
||||
double fplo, fphi, theta0;
|
||||
double *afile, *efile, *ffile;
|
||||
double *e2file, *f2file;
|
||||
double delta, invdelta, deltasq6;
|
||||
double *ang, *e, *de, *f, *df, *e2, *f2;
|
||||
};
|
||||
|
||||
struct ParamTable {
|
||||
int tablenamelength; // length of table name
|
||||
char *tablename; // name of associated angular table
|
||||
int keywordlength; // length of key in table
|
||||
char *keyword; // key in table
|
||||
int tabstyle, tablength; // length of interpolation table (not ninput) and style
|
||||
Table *angtable; // angle table
|
||||
};
|
||||
|
||||
protected:
|
||||
ParamTable *table_params; // tabulated parameter set for an I-J-K interaction
|
||||
|
||||
void read_file(char *) override;
|
||||
void threebody_table(Param *, Param *, ParamTable *, double, double, double *, double *,
|
||||
double *, double *, int, double &);
|
||||
|
||||
void read_table(Table *, char *, char *);
|
||||
void spline_table(Table *);
|
||||
void compute_table(Table *, int length);
|
||||
void bcast_table(Table *);
|
||||
void null_table(Table *);
|
||||
void free_table(Table *);
|
||||
void free_param(ParamTable *);
|
||||
void param_extract(Table *, char *);
|
||||
void spline(double *, double *, int, double, double, double *);
|
||||
double splint(double *, double *, double *, int, double);
|
||||
void uf_lookup(ParamTable *, double, double &, double &);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
809
src/MANYBODY/pair_threebody_table.cpp
Normal file
809
src/MANYBODY/pair_threebody_table.cpp
Normal file
@ -0,0 +1,809 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Christoph Scherer (MPIP Mainz)
|
||||
scherer@mpip-mainz.mpg.de
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_threebody_table.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "potential_file_reader.h"
|
||||
#include "table_file_reader.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using MathConst::MY_PI;
|
||||
|
||||
#define DELTA 4
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairThreebodyTable::PairThreebodyTable(LAMMPS *lmp) :
|
||||
Pair(lmp), params(nullptr), neighshort(nullptr)
|
||||
{
|
||||
single_enable = 0;
|
||||
restartinfo = 0;
|
||||
one_coeff = 1;
|
||||
manybody_flag = 1;
|
||||
centroidstressflag = CENTROID_NOTAVAIL;
|
||||
|
||||
maxshort = 10;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
check if allocated, since class can be destructed when incomplete
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
PairThreebodyTable::~PairThreebodyTable()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
for (int m = 0; m < nparams; m++) free_param(¶ms[m]); // free_param will call free_table
|
||||
memory->sfree(params);
|
||||
memory->destroy(elem3param);
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(neighshort);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::compute(int eflag, int vflag)
|
||||
{
|
||||
int i, j, k, ii, jj, kk, inum, jnum, jnumm1;
|
||||
int itype, jtype, ktype, ijparam, ijkparam;
|
||||
tagint itag, jtag;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl;
|
||||
double rsq, rsq1, rsq2;
|
||||
double delr1[3], delr2[3], fi[3], fj[3], fk[3];
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
tagint *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
// loop over full neighbor list of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itag = tag[i];
|
||||
itype = map[type[i]];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
// two-body interactions, skip half of them
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
int numshort = 0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
jtype = map[type[j]];
|
||||
ijparam = elem3param[itype][jtype][jtype];
|
||||
if (rsq >= params[ijparam].cutsq) {
|
||||
continue;
|
||||
} else {
|
||||
neighshort[numshort++] = j;
|
||||
if (numshort >= maxshort) {
|
||||
maxshort += maxshort / 2;
|
||||
memory->grow(neighshort, maxshort, "pair:neighshort");
|
||||
}
|
||||
}
|
||||
|
||||
jtag = tag[j];
|
||||
if (itag > jtag) {
|
||||
if ((itag + jtag) % 2 == 0) continue;
|
||||
} else if (itag < jtag) {
|
||||
if ((itag + jtag) % 2 == 1) continue;
|
||||
} else {
|
||||
if (x[j][2] < ztmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
|
||||
}
|
||||
|
||||
//two-body interactions are not computed
|
||||
}
|
||||
|
||||
jnumm1 = numshort - 1;
|
||||
|
||||
for (jj = 0; jj < jnumm1; jj++) {
|
||||
j = neighshort[jj];
|
||||
jtype = map[type[j]];
|
||||
ijparam = elem3param[itype][jtype][jtype];
|
||||
delr1[0] = x[j][0] - xtmp;
|
||||
delr1[1] = x[j][1] - ytmp;
|
||||
delr1[2] = x[j][2] - ztmp;
|
||||
rsq1 = delr1[0] * delr1[0] + delr1[1] * delr1[1] + delr1[2] * delr1[2];
|
||||
|
||||
double fjxtmp, fjytmp, fjztmp;
|
||||
fjxtmp = fjytmp = fjztmp = 0.0;
|
||||
|
||||
for (kk = jj + 1; kk < numshort; kk++) {
|
||||
k = neighshort[kk];
|
||||
ktype = map[type[k]];
|
||||
ijkparam = elem3param[itype][jtype][ktype];
|
||||
|
||||
delr2[0] = x[k][0] - xtmp;
|
||||
delr2[1] = x[k][1] - ytmp;
|
||||
delr2[2] = x[k][2] - ztmp;
|
||||
rsq2 = delr2[0] * delr2[0] + delr2[1] * delr2[1] + delr2[2] * delr2[2];
|
||||
|
||||
threebody(¶ms[ijkparam], rsq1, rsq2, delr1, delr2, fi, fj, fk, eflag, evdwl);
|
||||
|
||||
fxtmp += fi[0];
|
||||
fytmp += fi[1];
|
||||
fztmp += fi[2];
|
||||
fjxtmp += fj[0];
|
||||
fjytmp += fj[1];
|
||||
fjztmp += fj[2];
|
||||
f[k][0] += fk[0];
|
||||
f[k][1] += fk[1];
|
||||
f[k][2] += fk[2];
|
||||
|
||||
if (evflag) ev_tally3(i, j, k, evdwl, 0.0, fj, fk, delr1, delr2);
|
||||
}
|
||||
f[j][0] += fjxtmp;
|
||||
f[j][1] += fjytmp;
|
||||
f[j][2] += fjztmp;
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int np1 = atom->ntypes + 1;
|
||||
|
||||
memory->create(setflag, np1, np1, "pair:setflag");
|
||||
memory->create(cutsq, np1, np1, "pair:cutsq");
|
||||
memory->create(neighshort, maxshort, "pair:neighshort");
|
||||
map = new int[np1];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::settings(int narg, char ** /*arg*/)
|
||||
{
|
||||
if (narg != 0) error->all(FLERR, "Illegal pair_style command");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::coeff(int narg, char **arg)
|
||||
{
|
||||
if (!allocated) allocate();
|
||||
|
||||
map_element2type(narg - 3, arg + 3);
|
||||
|
||||
// read potential file and initialize potential parameters
|
||||
|
||||
if (params) {
|
||||
for (int m = 0; m < nparams; m++) free_param(¶ms[m]); // free_param will call free_table
|
||||
memory->sfree(params);
|
||||
params = nullptr;
|
||||
}
|
||||
read_file(arg[2]);
|
||||
setup_params();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::init_style()
|
||||
{
|
||||
if (atom->tag_enable == 0) error->all(FLERR, "Pair style threebody/table requires atom IDs");
|
||||
if (force->newton_pair == 0)
|
||||
error->all(FLERR, "Pair style threebody/table requires newton pair on");
|
||||
|
||||
// need a full neighbor list
|
||||
|
||||
neighbor->add_request(this, NeighConst::REQ_FULL);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairThreebodyTable::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
|
||||
|
||||
return cutmax;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::read_file(char *file)
|
||||
{
|
||||
params = nullptr;
|
||||
nparams = maxparam = 0;
|
||||
|
||||
// open file on proc 0
|
||||
|
||||
if (comm->me == 0) {
|
||||
PotentialFileReader reader(lmp, file, "threebody", unit_convert_flag);
|
||||
char *line;
|
||||
|
||||
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
|
||||
std::string iname = values.next_string();
|
||||
std::string jname = values.next_string();
|
||||
std::string kname = values.next_string();
|
||||
|
||||
// ielement,jelement,kelement = 1st args
|
||||
// if all 3 args are in element list, then parse this line
|
||||
// else skip to next entry in file
|
||||
int ielement, jelement, kelement;
|
||||
|
||||
for (ielement = 0; ielement < nelements; ielement++)
|
||||
if (iname == elements[ielement]) break;
|
||||
if (ielement == nelements) continue;
|
||||
for (jelement = 0; jelement < nelements; jelement++)
|
||||
if (jname == elements[jelement]) break;
|
||||
if (jelement == nelements) continue;
|
||||
for (kelement = 0; kelement < nelements; kelement++)
|
||||
if (kname == elements[kelement]) break;
|
||||
if (kelement == nelements) continue;
|
||||
|
||||
// load up parameter settings and error check their values
|
||||
|
||||
if (nparams == maxparam) {
|
||||
maxparam += DELTA;
|
||||
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
|
||||
|
||||
// make certain all addional allocated storage is initialized
|
||||
// to avoid false positives when checking with valgrind
|
||||
|
||||
memset(params + nparams, 0, DELTA * sizeof(Param));
|
||||
}
|
||||
|
||||
params[nparams].ielement = ielement;
|
||||
params[nparams].jelement = jelement;
|
||||
params[nparams].kelement = kelement;
|
||||
// if jelement = kelement, symmetric is true, if not then it is false
|
||||
params[nparams].symmetric = false;
|
||||
if (params[nparams].jelement == params[nparams].kelement) params[nparams].symmetric = true;
|
||||
// read cut off
|
||||
params[nparams].cut = values.next_double();
|
||||
|
||||
// read parameters of angle table
|
||||
std::string name = values.next_string();
|
||||
params[nparams].tablenamelength = name.length() + 1;
|
||||
params[nparams].tablename = utils::strdup(name);
|
||||
|
||||
name = values.next_string();
|
||||
params[nparams].keywordlength = name.length() + 1;
|
||||
params[nparams].keyword = utils::strdup(name);
|
||||
|
||||
name = values.next_string();
|
||||
if (name != "linear") error->all(FLERR, "Unknown table style {} in threebody table", name);
|
||||
params[nparams].tablength = values.next_int();
|
||||
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
if (params[nparams].cut < 0.0 || params[nparams].tablength < 0.0)
|
||||
error->one(FLERR, "Illegal threebody/table parameters");
|
||||
|
||||
nparams++;
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
|
||||
|
||||
if (comm->me != 0)
|
||||
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
|
||||
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
|
||||
|
||||
// for each set of parameters, broadcast table name and keyword and read threebody table
|
||||
for (int m = 0; m < nparams; ++m) {
|
||||
if (comm->me != 0) {
|
||||
params[m].tablename = new char[params[m].tablenamelength];
|
||||
params[m].keyword = new char[params[m].keywordlength];
|
||||
}
|
||||
|
||||
MPI_Bcast(params[m].tablename, params[m].tablenamelength, MPI_CHAR, 0, world);
|
||||
MPI_Bcast(params[m].keyword, params[m].keywordlength, MPI_CHAR, 0, world);
|
||||
|
||||
// initialize threebodytable
|
||||
memory->create(params[m].mltable, 1, "param:threebodytable");
|
||||
null_table(params[m].mltable);
|
||||
|
||||
//call read_table to read corresponding tabulated threebody file (only called by process 0)
|
||||
if (comm->me == 0) {
|
||||
read_table(params[m].mltable, params[m].tablename, params[m].keyword, params[m].symmetric);
|
||||
}
|
||||
|
||||
// broadcast read in threebodytable to all processes
|
||||
bcast_table(params[m].mltable, params[m].symmetric);
|
||||
|
||||
// error check on table parameters
|
||||
if (params[m].mltable->ninput <= 1) error->one(FLERR, "Invalid threebody table length");
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::setup_params()
|
||||
{
|
||||
int i, j, k, m, n;
|
||||
double rtmp;
|
||||
|
||||
// set elem3param for all triplet combinations
|
||||
// must be a single exact match to lines read from file
|
||||
// do not allow for ACB in place of ABC
|
||||
|
||||
memory->destroy(elem3param);
|
||||
memory->create(elem3param, nelements, nelements, nelements, "pair:elem3param");
|
||||
|
||||
for (i = 0; i < nelements; i++)
|
||||
for (j = 0; j < nelements; j++)
|
||||
for (k = 0; k < nelements; k++) {
|
||||
n = -1;
|
||||
for (m = 0; m < nparams; m++) {
|
||||
if (i == params[m].ielement && j == params[m].jelement && k == params[m].kelement) {
|
||||
if (n >= 0) error->all(FLERR, "Potential file has duplicate entry");
|
||||
n = m;
|
||||
}
|
||||
}
|
||||
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
|
||||
elem3param[i][j][k] = n;
|
||||
}
|
||||
|
||||
// compute parameter values derived from inputs
|
||||
|
||||
// set cutsq using shortcut to reduce neighbor list for accelerated
|
||||
// calculations. cut must remain unchanged as it is a potential parameter
|
||||
// (cut = a)
|
||||
|
||||
for (m = 0; m < nparams; m++) {
|
||||
rtmp = params[m].cut;
|
||||
params[m].cutsq = rtmp * rtmp;
|
||||
}
|
||||
|
||||
// set cutmax to max of all params
|
||||
|
||||
cutmax = 0.0;
|
||||
for (m = 0; m < nparams; m++) {
|
||||
rtmp = sqrt(params[m].cutsq);
|
||||
if (rtmp > cutmax) cutmax = rtmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
read table file, only called by proc 0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::read_table(Table *tb, char *file, char *keyword, bool symmetric)
|
||||
{
|
||||
TableFileReader reader(lmp, file, "threebodytable");
|
||||
|
||||
char *line = reader.find_section_start(keyword);
|
||||
|
||||
if (!line) { error->one(FLERR, "Did not find keyword in table file"); }
|
||||
|
||||
// read args on 2nd line of section
|
||||
// allocate table arrays for file values
|
||||
|
||||
line = reader.next_line();
|
||||
param_extract(tb, line);
|
||||
|
||||
// if it is a symmetric threebody interaction, less table entries are required
|
||||
if (symmetric == true) {
|
||||
memory->create(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r12file");
|
||||
memory->create(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r13file");
|
||||
memory->create(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:thetafile");
|
||||
memory->create(tb->f11file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f11file");
|
||||
memory->create(tb->f12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f12file");
|
||||
memory->create(tb->f21file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f21file");
|
||||
memory->create(tb->f22file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f22file");
|
||||
memory->create(tb->f31file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f31file");
|
||||
memory->create(tb->f32file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f32file");
|
||||
memory->create(tb->efile, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:efile");
|
||||
}
|
||||
// else, more (full) table entries are required
|
||||
else {
|
||||
memory->create(tb->r12file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:r12file");
|
||||
memory->create(tb->r13file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:r13file");
|
||||
memory->create(tb->thetafile, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:thetafile");
|
||||
memory->create(tb->f11file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f11file");
|
||||
memory->create(tb->f12file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f12file");
|
||||
memory->create(tb->f21file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f21file");
|
||||
memory->create(tb->f22file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f22file");
|
||||
memory->create(tb->f31file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f31file");
|
||||
memory->create(tb->f32file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f32file");
|
||||
memory->create(tb->efile, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:efile");
|
||||
}
|
||||
|
||||
// read threebody table values from file
|
||||
|
||||
int cerror = 0;
|
||||
reader.skip_line();
|
||||
// if it is a symmetric threebody interaction, less table entries are required
|
||||
if (symmetric == true) {
|
||||
for (int i = 0; i < tb->ninput * tb->ninput * (tb->ninput + 1); i++) {
|
||||
line = reader.next_line(11);
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
values.next_int();
|
||||
tb->r12file[i] = values.next_double();
|
||||
tb->r13file[i] = values.next_double();
|
||||
tb->thetafile[i] = values.next_double();
|
||||
tb->f11file[i] = values.next_double();
|
||||
tb->f12file[i] = values.next_double();
|
||||
tb->f21file[i] = values.next_double();
|
||||
tb->f22file[i] = values.next_double();
|
||||
tb->f31file[i] = values.next_double();
|
||||
tb->f32file[i] = values.next_double();
|
||||
tb->efile[i] = values.next_double();
|
||||
} catch (TokenizerException &) {
|
||||
++cerror;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < 2 * tb->ninput * tb->ninput * tb->ninput; i++) {
|
||||
line = reader.next_line(11);
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
values.next_int();
|
||||
tb->r12file[i] = values.next_double();
|
||||
tb->r13file[i] = values.next_double();
|
||||
tb->thetafile[i] = values.next_double();
|
||||
tb->f11file[i] = values.next_double();
|
||||
tb->f12file[i] = values.next_double();
|
||||
tb->f21file[i] = values.next_double();
|
||||
tb->f22file[i] = values.next_double();
|
||||
tb->f31file[i] = values.next_double();
|
||||
tb->f32file[i] = values.next_double();
|
||||
tb->efile[i] = values.next_double();
|
||||
} catch (TokenizerException &) {
|
||||
++cerror;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// warn if data was read incompletely, e.g. columns were missing
|
||||
|
||||
if (cerror)
|
||||
error->warning(FLERR, "{} of {} lines in table incomplete or could not be parsed", cerror,
|
||||
tb->ninput);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract attributes from parameter line in table section
|
||||
format of line: N value FP fplo fphi EQ theta0
|
||||
N is required, other params are optional
|
||||
|
||||
only called by read_table, only called by proc 0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::param_extract(Table *tb, char *line)
|
||||
{
|
||||
tb->ninput = 0;
|
||||
tb->rmin = 0.0;
|
||||
tb->rmax = 0.0;
|
||||
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
|
||||
while (values.has_next()) {
|
||||
std::string word = values.next_string();
|
||||
|
||||
if (word == "N") {
|
||||
tb->ninput = values.next_int();
|
||||
} else if (word == "rmin") {
|
||||
tb->rmin = values.next_double();
|
||||
} else if (word == "rmax") {
|
||||
tb->rmax = values.next_double();
|
||||
} else {
|
||||
error->one(FLERR, "Invalid keyword {} in angle table parameters", word);
|
||||
}
|
||||
}
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
if (tb->ninput == 0) error->one(FLERR, "threebodytable parameters did not set N");
|
||||
if (tb->rmin == 0.0) error->one(FLERR, "threebodytable parameters did not set rmin");
|
||||
if (tb->rmax == 0.0) error->one(FLERR, "threebodytable parameters did not set rmax");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
broadcast read-in table info from proc 0 to other procs
|
||||
this function communicates these values in Table:
|
||||
ninput,afile,efile,ffile,fpflag,fplo,fphi,theta0
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::bcast_table(Table *tb, bool symmetric)
|
||||
{
|
||||
MPI_Bcast(&tb->ninput, 1, MPI_INT, 0, world);
|
||||
|
||||
int me;
|
||||
MPI_Comm_rank(world, &me);
|
||||
if (me > 0) {
|
||||
// if it is a symmetric threebody interaction, less table entries are required
|
||||
if (symmetric == true) {
|
||||
memory->create(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r12file");
|
||||
memory->create(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r13file");
|
||||
memory->create(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1),
|
||||
"mltable:thetafile");
|
||||
memory->create(tb->f11file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f11file");
|
||||
memory->create(tb->f12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f12file");
|
||||
memory->create(tb->f21file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f21file");
|
||||
memory->create(tb->f22file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f22file");
|
||||
memory->create(tb->f31file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f31file");
|
||||
memory->create(tb->f32file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:f32file");
|
||||
memory->create(tb->efile, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:efile");
|
||||
}
|
||||
// else, more (full) table entries are required
|
||||
else {
|
||||
memory->create(tb->r12file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:r12file");
|
||||
memory->create(tb->r13file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:r13file");
|
||||
memory->create(tb->thetafile, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:thetafile");
|
||||
memory->create(tb->f11file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f11file");
|
||||
memory->create(tb->f12file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f12file");
|
||||
memory->create(tb->f21file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f21file");
|
||||
memory->create(tb->f22file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f22file");
|
||||
memory->create(tb->f31file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f31file");
|
||||
memory->create(tb->f32file, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:f32file");
|
||||
memory->create(tb->efile, 2 * tb->ninput * tb->ninput * tb->ninput, "mltable:efile");
|
||||
}
|
||||
}
|
||||
|
||||
// if it is a symmetric threebody interaction, less table entries are required
|
||||
if (symmetric == true) {
|
||||
MPI_Bcast(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f11file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f12file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f21file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f22file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f31file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f32file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->efile, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world);
|
||||
}
|
||||
// else, more (full) table entries are required
|
||||
else {
|
||||
MPI_Bcast(tb->r12file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->r13file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->thetafile, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f11file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f12file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f21file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f22file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f31file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->f32file, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->efile, 2 * tb->ninput * tb->ninput * tb->ninput, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
|
||||
MPI_Bcast(&tb->rmin, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&tb->rmax, 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::free_param(Param *pm)
|
||||
{
|
||||
// call free_table to destroy associated threebodytables
|
||||
free_table(pm->mltable);
|
||||
// then destroy associated threebodytable
|
||||
delete[] pm->tablename;
|
||||
delete[] pm->keyword;
|
||||
memory->sfree(pm->mltable);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::free_table(Table *tb)
|
||||
{
|
||||
memory->destroy(tb->r12file);
|
||||
memory->destroy(tb->r13file);
|
||||
memory->destroy(tb->thetafile);
|
||||
memory->destroy(tb->f11file);
|
||||
memory->destroy(tb->f12file);
|
||||
memory->destroy(tb->f21file);
|
||||
memory->destroy(tb->f22file);
|
||||
memory->destroy(tb->f31file);
|
||||
memory->destroy(tb->f32file);
|
||||
memory->destroy(tb->efile);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::null_table(Table *tb)
|
||||
{
|
||||
tb->r12file = tb->r13file = tb->thetafile = nullptr;
|
||||
tb->f11file = tb->f12file = nullptr;
|
||||
tb->f21file = tb->f22file = nullptr;
|
||||
tb->f31file = tb->f32file = nullptr;
|
||||
tb->efile = nullptr;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
calculate potential u and force f at angle x
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double theta, double &f11,
|
||||
double &f12, double &f21, double &f22, double &f31, double &f32,
|
||||
double &u)
|
||||
{
|
||||
int i, itable, nr12, nr13, ntheta;
|
||||
double dr, dtheta;
|
||||
dr = (pm->mltable->rmax - pm->mltable->rmin) / (pm->mltable->ninput - 1);
|
||||
dtheta = (180.0 - 0.0) / (pm->mltable->ninput * 2);
|
||||
|
||||
//lookup scheme
|
||||
|
||||
// if it is a symmetric threebody interaction, less table entries are required
|
||||
if (pm->symmetric == true) {
|
||||
nr12 = (r12 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr;
|
||||
if (r12 == (pm->mltable->rmin - 0.5 * dr)) { nr12 = 0; }
|
||||
nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr;
|
||||
if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; }
|
||||
nr13 -= nr12;
|
||||
ntheta = (theta - 0.00000001) / dtheta;
|
||||
if (theta == 180.0) { ntheta = 79; }
|
||||
itable = 0;
|
||||
for (i = 0; i < nr12; i++) { itable += (pm->mltable->ninput - i); }
|
||||
itable += nr13;
|
||||
itable *= (pm->mltable->ninput * 2);
|
||||
itable += ntheta;
|
||||
} else {
|
||||
// else, more (full) table entries are required
|
||||
nr12 = (r12 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr;
|
||||
if (r12 == (pm->mltable->rmin - 0.5 * dr)) { nr12 = 0; }
|
||||
nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr;
|
||||
if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; }
|
||||
ntheta = (theta - 0.00000001) / dtheta;
|
||||
if (theta == 180.0) { ntheta = 79; }
|
||||
itable = nr12 * (pm->mltable->ninput);
|
||||
itable += nr13;
|
||||
itable *= (pm->mltable->ninput * 2);
|
||||
itable += ntheta;
|
||||
}
|
||||
|
||||
f11 = pm->mltable->f11file[itable];
|
||||
f12 = pm->mltable->f12file[itable];
|
||||
f21 = pm->mltable->f21file[itable];
|
||||
f22 = pm->mltable->f22file[itable];
|
||||
f31 = pm->mltable->f31file[itable];
|
||||
f32 = pm->mltable->f32file[itable];
|
||||
u = pm->mltable->efile[itable];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairThreebodyTable::threebody(Param *paramijk, double rsq1, double rsq2, double *delr1,
|
||||
double *delr2, double *fi, double *fj, double *fk, int eflag,
|
||||
double &eng)
|
||||
{
|
||||
double r12, r13, theta, rinv, cs;
|
||||
|
||||
double f11, f12, f21, f22, f31, f32, u, temp;
|
||||
bool swapped;
|
||||
double dr;
|
||||
dr = (paramijk->mltable->rmax - paramijk->mltable->rmin) / (paramijk->mltable->ninput - 1);
|
||||
|
||||
//if swap indices or not
|
||||
swapped = false;
|
||||
|
||||
r12 = sqrt(rsq1);
|
||||
r13 = sqrt(rsq2);
|
||||
rinv = 1.0 / (r12 * r13);
|
||||
cs = (delr1[0] * delr2[0] + delr1[1] * delr2[1] + delr1[2] * delr2[2]) * rinv;
|
||||
//compute angle between r12 and r13 in degrees
|
||||
theta = acos(cs) * 180.0 / MY_PI;
|
||||
|
||||
//if r12 > r13 swap them, as in lookup table always r13 > r12 do to symmetry reasons
|
||||
if (r12 > r13) {
|
||||
temp = r12;
|
||||
r12 = r13;
|
||||
r13 = temp;
|
||||
swapped = true;
|
||||
}
|
||||
|
||||
//look up forces and energy in table belonging to parameter set paramijk
|
||||
|
||||
//only do lookup and add three-body interactions if r12 and r13 are both between rmin and rmax
|
||||
|
||||
if ((r12 >= (paramijk->mltable->rmin - 0.5 * dr)) &&
|
||||
(r13 <= (paramijk->mltable->rmax + 0.5 * dr)) &&
|
||||
(r13 >= (paramijk->mltable->rmin - 0.5 * dr)) &&
|
||||
(r13 <= (paramijk->mltable->rmax + 0.5 * dr))) {
|
||||
uf_lookup(paramijk, r12, r13, theta, f11, f12, f21, f22, f31, f32, u);
|
||||
} else {
|
||||
f11 = f12 = f21 = f22 = f31 = f32 = u = 0.0;
|
||||
}
|
||||
|
||||
// if the indices have been swapped, swap them back
|
||||
if (swapped == true) {
|
||||
temp = r12;
|
||||
r12 = r13;
|
||||
r13 = temp;
|
||||
temp = f11;
|
||||
f11 = f12;
|
||||
f12 = temp;
|
||||
temp = f21;
|
||||
f21 = f31;
|
||||
f31 = temp;
|
||||
temp = f22;
|
||||
f22 = -f32;
|
||||
f32 = -temp;
|
||||
}
|
||||
|
||||
fi[0] = delr1[0] * f11 + delr2[0] * f12;
|
||||
fi[1] = delr1[1] * f11 + delr2[1] * f12;
|
||||
fi[2] = delr1[2] * f11 + delr2[2] * f12;
|
||||
|
||||
fj[0] = delr1[0] * f21 + (delr2[0] - delr1[0]) * f22;
|
||||
fj[1] = delr1[1] * f21 + (delr2[1] - delr1[1]) * f22;
|
||||
fj[2] = delr1[2] * f21 + (delr2[2] - delr1[2]) * f22;
|
||||
|
||||
fk[0] = delr2[0] * f31 + (delr2[0] - delr1[0]) * f32;
|
||||
fk[1] = delr2[1] * f31 + (delr2[1] - delr1[1]) * f32;
|
||||
fk[2] = delr2[2] * f31 + (delr2[2] - delr1[2]) * f32;
|
||||
|
||||
if (eflag) eng = u;
|
||||
}
|
||||
89
src/MANYBODY/pair_threebody_table.h
Normal file
89
src/MANYBODY/pair_threebody_table.h
Normal file
@ -0,0 +1,89 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
// clang-format off
|
||||
PairStyle(threebody/table,PairThreebodyTable);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_THREEBODY_TABLE_H
|
||||
#define LMP_PAIR_THREEBODY_TABLE_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairThreebodyTable : public Pair {
|
||||
public:
|
||||
PairThreebodyTable(class LAMMPS *);
|
||||
~PairThreebodyTable() override;
|
||||
void compute(int, int) override;
|
||||
void coeff(int, char **) override;
|
||||
double init_one(int, int) override;
|
||||
void init_style() override;
|
||||
|
||||
static constexpr int NPARAMS_PER_LINE = 8;
|
||||
|
||||
// no write or read from binary restart file
|
||||
|
||||
// struct for threebody/table
|
||||
struct Table {
|
||||
int ninput;
|
||||
double rmin, rmax;
|
||||
double *r12file, *r13file, *thetafile, *f11file, *f12file, *f21file, *f22file, *f31file,
|
||||
*f32file, *efile;
|
||||
};
|
||||
|
||||
struct Param {
|
||||
double cut, cutsq;
|
||||
int ielement, jelement, kelement;
|
||||
bool symmetric; // whether it is a symmetric table or not
|
||||
int tablenamelength; // length of table name
|
||||
char *tablename; // name of associated angular table
|
||||
int keywordlength; // length of key in table
|
||||
char *keyword; // key in table
|
||||
int tabstyle, tablength; // length of interpolation table (not ninput) and style
|
||||
Table *mltable; // threebody table
|
||||
};
|
||||
|
||||
protected:
|
||||
double cutmax; // max cutoff for all elements
|
||||
Param *params; // parameter set for an I-J-K interaction
|
||||
int maxshort; // size of short neighbor list array
|
||||
int *neighshort; // short neighbor list array
|
||||
|
||||
void settings(int, char **) override;
|
||||
virtual void allocate();
|
||||
void read_file(char *);
|
||||
virtual void setup_params();
|
||||
void threebody(Param *, double, double, double *, double *, double *, double *, double *, int,
|
||||
double &);
|
||||
|
||||
void read_table(Table *, char *, char *, bool);
|
||||
void bcast_table(Table *, bool);
|
||||
void null_table(Table *);
|
||||
|
||||
void free_table(Table *);
|
||||
void free_param(Param *);
|
||||
|
||||
void param_extract(Table *, char *);
|
||||
|
||||
void uf_lookup(Param *, double, double, double, double &, double &, double &, double &, double &,
|
||||
double &, double &);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
1875
unittest/force-styles/tests/1-1-1.table
Normal file
1875
unittest/force-styles/tests/1-1-1.table
Normal file
File diff suppressed because it is too large
Load Diff
3459
unittest/force-styles/tests/1-1-2.table
Normal file
3459
unittest/force-styles/tests/1-1-2.table
Normal file
File diff suppressed because it is too large
Load Diff
97
unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml
Normal file
97
unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml
Normal file
@ -0,0 +1,97 @@
|
||||
---
|
||||
lammps_version: 4 May 2022
|
||||
tags: generated
|
||||
date_generated: Wed Jun 1 15:17:22 2022
|
||||
epsilon: 1e-12
|
||||
skip_tests: single
|
||||
prerequisites: ! |
|
||||
pair sw/angle/table
|
||||
pair table
|
||||
pre_commands: ! |
|
||||
variable units index real
|
||||
variable newton_pair delete
|
||||
variable newton_pair index on
|
||||
shell cp ${input_dir}/table_CG_CG_CG.txt .
|
||||
post_commands: ! ""
|
||||
input_file: in.metal
|
||||
pair_style: hybrid/overlay table linear 1200 sw/angle/table
|
||||
pair_coeff: ! |
|
||||
* * table ${input_dir}/table_CG_CG.txt VOTCA
|
||||
* * sw/angle/table ${input_dir}/spce.sw type type
|
||||
extract: ! ""
|
||||
natoms: 32
|
||||
init_vdwl: 2428.9633428241673
|
||||
init_coul: 0
|
||||
init_stress: ! |2-
|
||||
2.1799808252981104e+04 2.2098847758334741e+04 1.5227317495796946e+04 5.5915487284825185e+03 4.2330141376182507e+02 -2.1635648030093221e+03
|
||||
init_forces: ! |2
|
||||
1 -3.5272456427265922e+02 4.7896703188444229e+02 7.6160632749768567e+01
|
||||
2 -2.2734507500176406e+01 -1.5821645847684510e+02 1.1028261572305695e+02
|
||||
3 -3.6932223674029137e+02 1.7315628979419187e+03 -5.2288079251214720e+02
|
||||
4 3.1116044194470057e+01 -3.1009878149268530e+01 -6.8414290700926671e+01
|
||||
5 1.2157567668233901e+03 1.3361425827696451e+03 -2.4428700622833938e+02
|
||||
6 2.1498603159595226e+02 1.2379258234822432e+01 -1.8192892150594932e+02
|
||||
7 -5.8065242995364076e+02 6.3615913954340272e+02 -5.8940661342540871e+01
|
||||
8 -2.9330102622037936e+02 -1.4478456371145094e+02 3.4992669050834974e+02
|
||||
9 5.4581529315399578e+01 -1.0085658890730177e+02 -4.3539166606697755e+01
|
||||
10 -7.5328757518773557e+02 -1.9208550331031577e+03 -5.9929086772884966e+02
|
||||
11 -6.2073979508185595e+01 -9.3172505877146349e+01 8.8201736909256510e+01
|
||||
12 5.2022495622775352e+02 -6.4668600468680108e+02 -1.8086255931588799e+01
|
||||
13 -1.4637585277113917e+02 -2.4193749312797078e+01 -1.3497675472534843e+02
|
||||
14 2.2785633726795228e+02 -1.4050021950202930e+02 4.2957377860079254e+02
|
||||
15 -4.1589593903912913e+01 5.6849936807240290e+01 -5.3315771137404397e+01
|
||||
16 5.0207265346701280e+02 -4.3553084670415353e+02 2.2270110539464073e+02
|
||||
17 2.7243217976852867e+02 6.7842110608020960e+02 -1.8488293016613730e+02
|
||||
18 -1.6339467540544510e+03 -8.6208840396403559e+02 -5.2809809085219297e+02
|
||||
19 -1.8146991394127588e+03 -1.4248970633821093e+03 1.6246778777497133e+02
|
||||
20 5.0143947854312678e+01 3.0349353798587607e+01 -7.6753179337391444e+01
|
||||
21 1.1359392702527382e+03 6.7780617382057903e+02 -2.1777379118829096e+01
|
||||
22 2.6318213617558456e+01 -1.1442799194941128e+02 -4.0723882345600529e+01
|
||||
23 1.0173532367943421e+03 1.4870722398544501e+03 -3.3061556638580618e+02
|
||||
24 1.8951324945224176e+03 1.3655558041004167e+03 6.3746947970957035e+02
|
||||
25 2.1139286860441129e+02 -1.4343085616543428e+02 -2.4472193090284622e+02
|
||||
26 -6.6054117554868481e+02 -1.7214679588856484e+03 1.1872792057456782e+03
|
||||
27 3.8554823693482177e+02 -2.4263768110018356e+02 -1.4505783275426307e+01
|
||||
28 -1.6156920382667545e+02 3.2681073686927527e+02 4.0195534333261003e+02
|
||||
29 1.0269877810330740e+03 1.0972018261937728e+03 -4.9239365569732279e+01
|
||||
30 -7.7183246664884393e+01 -1.1163723935859770e+02 -5.6015149765282524e+02
|
||||
31 2.7330076741933460e+01 -6.2134053241130312e+02 3.7926314422192496e+02
|
||||
32 -1.8451713394504984e+03 -9.7754451225108528e+02 -6.8151426644039077e+01
|
||||
run_vdwl: 2428.764401023566
|
||||
run_coul: 0
|
||||
run_stress: ! |2-
|
||||
2.1807179009069081e+04 2.2096249577665836e+04 1.5217251424717178e+04 5.5876293741471409e+03 4.2481794037948595e+02 -2.1641073132805273e+03
|
||||
run_forces: ! |2
|
||||
1 -3.5657232205619187e+02 4.8051759317525114e+02 7.4589821279043520e+01
|
||||
2 -2.2890330563537752e+01 -1.5847929947391452e+02 1.1057024581491029e+02
|
||||
3 -3.7048240284136381e+02 1.7274339155425446e+03 -5.2721183867080663e+02
|
||||
4 3.1213113517876284e+01 -3.0972108629752253e+01 -6.8362160774369471e+01
|
||||
5 1.2131778379678060e+03 1.3324719919494626e+03 -2.4196312117000102e+02
|
||||
6 2.1588387878389867e+02 1.2017400433555963e+01 -1.8323099068613300e+02
|
||||
7 -5.8298397256607473e+02 6.3858638821865122e+02 -5.9272595884065503e+01
|
||||
8 -2.9450691866540427e+02 -1.4488601302098704e+02 3.5084177622838757e+02
|
||||
9 5.5071730152343321e+01 -1.0126346692429934e+02 -4.3948685147789718e+01
|
||||
10 -7.4709028759946739e+02 -1.9173399004644243e+03 -5.9495231666550546e+02
|
||||
11 -6.1959233740417524e+01 -9.3373357444258517e+01 8.7926257027673998e+01
|
||||
12 5.2278142223716191e+02 -6.4849088771234563e+02 -1.8113808074276363e+01
|
||||
13 -1.4714650249643290e+02 -2.5131629765603009e+01 -1.3479374072464537e+02
|
||||
14 2.2857040301009684e+02 -1.3768541975979431e+02 4.2806018886113947e+02
|
||||
15 -4.1477060277693703e+01 5.7115876564426109e+01 -5.3039366059682528e+01
|
||||
16 4.9944201304657383e+02 -4.3383072035483559e+02 2.2091297501303973e+02
|
||||
17 2.7228851840542382e+02 6.7799738753924669e+02 -1.8446508468678948e+02
|
||||
18 -1.6336221792482595e+03 -8.5972791234834551e+02 -5.2898505983077177e+02
|
||||
19 -1.8135957890859013e+03 -1.4238141052528933e+03 1.6225337159545035e+02
|
||||
20 5.0092828583367634e+01 3.0471251647078265e+01 -7.6722240263741099e+01
|
||||
21 1.1355438484696886e+03 6.7519841904221255e+02 -2.0182855479720033e+01
|
||||
22 2.6571960650017800e+01 -1.1420696745726380e+02 -4.0529746043707348e+01
|
||||
23 1.0263737261398123e+03 1.4932072283307180e+03 -3.2636823427367165e+02
|
||||
24 1.8952730712010357e+03 1.3642712022544688e+03 6.3847090965522716e+02
|
||||
25 2.1193565520738500e+02 -1.4318069871618528e+02 -2.4487119695300959e+02
|
||||
26 -6.5812362830257700e+02 -1.7188102078185152e+03 1.1859492616184705e+03
|
||||
27 3.8565010020982788e+02 -2.4238978364677456e+02 -1.4837594446360082e+01
|
||||
28 -1.6123154438363622e+02 3.2723676308792528e+02 4.0126611479067810e+02
|
||||
29 1.0251768164068674e+03 1.0954706244344939e+03 -4.9292343676448787e+01
|
||||
30 -7.7077129341229522e+01 -1.1123357160671468e+02 -5.6043482841374544e+02
|
||||
31 2.5726118886179297e+01 -6.2166714125994793e+02 3.8003828834174129e+02
|
||||
32 -1.8520137417071724e+03 -9.8551285056318022e+02 -6.9301402300522653e+01
|
||||
...
|
||||
98
unittest/force-styles/tests/atomic-pair-threebody_table.yaml
Normal file
98
unittest/force-styles/tests/atomic-pair-threebody_table.yaml
Normal file
@ -0,0 +1,98 @@
|
||||
---
|
||||
lammps_version: 4 May 2022
|
||||
tags: generated
|
||||
date_generated: Wed Jun 1 15:28:13 2022
|
||||
epsilon: 1e-05
|
||||
skip_tests: single
|
||||
prerequisites: ! |
|
||||
pair threebody/table
|
||||
pair table
|
||||
pre_commands: ! |
|
||||
variable units index real
|
||||
variable newton_pair delete
|
||||
variable newton_pair index on
|
||||
shell cp ${input_dir}/1-1-1.table .
|
||||
shell cp ${input_dir}/1-1-2.table .
|
||||
post_commands: ! ""
|
||||
input_file: in.metal
|
||||
pair_style: hybrid/overlay table linear 1200 threebody/table
|
||||
pair_coeff: ! |
|
||||
* * table ${input_dir}/table_CG_CG.txt VOTCA
|
||||
* * threebody/table ${input_dir}/spce2.3b type1 type2
|
||||
extract: ! ""
|
||||
natoms: 32
|
||||
init_vdwl: 1491.9850663210582
|
||||
init_coul: 0
|
||||
init_stress: ! |2-
|
||||
2.1388163370760823e+04 2.1664558645983379e+04 1.4729243404366314e+04 5.6495516964437775e+03 5.1637900223635859e+02 -2.2491014848350428e+03
|
||||
init_forces: ! |2
|
||||
1 -3.4809429741393029e+02 4.6567597414239913e+02 9.4441973687110405e+01
|
||||
2 -1.8412192720214428e+01 -1.6122507911305391e+02 1.1798397229543718e+02
|
||||
3 -3.6959552927057359e+02 1.7366664628134174e+03 -5.2433878744101696e+02
|
||||
4 1.9704833162904933e+01 -2.8473480842310366e+01 -7.6632873700899410e+01
|
||||
5 1.2286660791793993e+03 1.3189646599149826e+03 -2.5750328829062335e+02
|
||||
6 2.3230773636573508e+02 1.3236909769358112e+01 -1.6536673989911372e+02
|
||||
7 -5.9250555047871524e+02 6.4419772822168966e+02 -7.3421471775369668e+01
|
||||
8 -2.7993451614515635e+02 -1.3398848087050882e+02 3.4786818228776917e+02
|
||||
9 6.0102535250319256e+01 -8.2958743522890487e+01 -4.2554780357129808e+01
|
||||
10 -7.7091710125156442e+02 -1.9316366702146124e+03 -5.7574889508217711e+02
|
||||
11 -6.6141742740805213e+01 -8.9045678804585251e+01 9.2760444861105910e+01
|
||||
12 5.3526455218407943e+02 -6.5780369404995440e+02 -9.4845914108869884e+00
|
||||
13 -1.6564083090865202e+02 -2.1519923677640854e+01 -1.3123092115579615e+02
|
||||
14 2.3337854905367016e+02 -1.1893053800382607e+02 4.3706653618942192e+02
|
||||
15 -4.7292738629250245e+01 6.1031002391688908e+01 -4.2739580007555375e+01
|
||||
16 5.0083107755149183e+02 -4.2750321084667428e+02 2.3013161197871258e+02
|
||||
17 2.9344728675986164e+02 6.8063155134388398e+02 -1.9478515772339574e+02
|
||||
18 -1.6545067282255825e+03 -8.7660521680902912e+02 -5.2486018536431391e+02
|
||||
19 -1.7992841953748712e+03 -1.4223424241054529e+03 1.3975194023223264e+02
|
||||
20 5.3624371129881432e+01 4.2727424976681945e+01 -7.2478104483156997e+01
|
||||
21 1.0897088707455639e+03 7.2603975137317627e+02 -5.1120443430894568e+01
|
||||
22 2.9564358575254730e+01 -1.1955500923091164e+02 -5.6658561557696522e+01
|
||||
23 1.0024095663866029e+03 1.4815830194767184e+03 -3.4241061954729582e+02
|
||||
24 1.8958818608698684e+03 1.3513089573990835e+03 6.4474764645157461e+02
|
||||
25 2.1916799984257568e+02 -1.3480959959762640e+02 -2.5143909195778633e+02
|
||||
26 -6.7819370861387029e+02 -1.7332731917983826e+03 1.1759045104066886e+03
|
||||
27 3.9534539412579926e+02 -2.5831579543312483e+02 -3.0068663848303565e+01
|
||||
28 -1.7049475634836011e+02 3.0557653380225258e+02 3.9667516538156866e+02
|
||||
29 1.0124807267829628e+03 1.0704993768753102e+03 -6.6827000765975839e+01
|
||||
30 -6.0958426550901464e+01 -1.2048001317378247e+02 -5.4659311407722760e+02
|
||||
31 2.7988221050038256e+01 -6.0346016941066136e+02 3.8578368661473195e+02
|
||||
32 -1.8079021293566357e+03 -9.7620852873268325e+02 -2.6848072654532874e+01
|
||||
run_vdwl: 1491.7709826398068
|
||||
run_coul: 0
|
||||
run_stress: ! |2-
|
||||
2.1425064502585417e+04 2.1669158423280089e+04 1.4747818832883342e+04 5.6707077802984231e+03 5.4904652273389740e+02 -2.2288016407219948e+03
|
||||
run_forces: ! |2
|
||||
1 -3.5191405522221208e+02 4.6719496439671821e+02 9.5160203317446630e+01
|
||||
2 -1.8514912253517174e+01 -1.6167539395717191e+02 1.1813103212569202e+02
|
||||
3 -3.7131956856014324e+02 1.7323797805692438e+03 -5.3074350641727199e+02
|
||||
4 1.9443697338922082e+01 -2.8705869694043574e+01 -7.6619749545143492e+01
|
||||
5 1.2260941309164493e+03 1.3152677694284498e+03 -2.5516802038260235e+02
|
||||
6 2.3322089165638701e+02 1.2862869405899247e+01 -1.6667928243690179e+02
|
||||
7 -5.9486664911064634e+02 6.4667229208938795e+02 -7.3770515588557785e+01
|
||||
8 -2.8115980510834169e+02 -1.3411610929266357e+02 3.4879265918581427e+02
|
||||
9 6.0619815294130184e+01 -8.3407515617420799e+01 -4.2999027435430953e+01
|
||||
10 -7.6429388783007221e+02 -1.9278828213213337e+03 -5.7135063514687909e+02
|
||||
11 -6.6016349831195143e+01 -8.9273528569753211e+01 9.2448231530143218e+01
|
||||
12 5.3785966085411701e+02 -6.5959880855240010e+02 -9.4946135547062909e+00
|
||||
13 -1.6641765862714641e+02 -2.2476690068776783e+01 -1.3101354808256795e+02
|
||||
14 2.3412385380018750e+02 -1.1604834858772836e+02 4.3551831586668749e+02
|
||||
15 -4.7170704261276846e+01 6.1314362285436715e+01 -4.2456307464845409e+01
|
||||
16 4.9810572379063296e+02 -4.2579585263982955e+02 2.2832018077001956e+02
|
||||
17 2.8312588438935751e+02 6.6818248758640414e+02 -2.0387201756671962e+02
|
||||
18 -1.6531905535811466e+03 -8.7325660040293678e+02 -5.2578287557408930e+02
|
||||
19 -1.7717540076086393e+03 -1.4078190114507104e+03 1.5867151421748673e+02
|
||||
20 5.3592891283041261e+01 4.2790401518478163e+01 -7.2482601253922383e+01
|
||||
21 1.0838596264146349e+03 7.2319003723887010e+02 -4.8376931553139812e+01
|
||||
22 2.9966891580030588e+01 -1.1952124544144279e+02 -5.6490752521580006e+01
|
||||
23 1.0114992761256929e+03 1.4877360943225701e+03 -3.3815548489754946e+02
|
||||
24 1.8960185852616901e+03 1.3500190426973886e+03 6.4578239684187895e+02
|
||||
25 2.1972099306244061e+02 -1.3453156341222555e+02 -2.5160650947118697e+02
|
||||
26 -6.7764066859903267e+02 -1.7320367725724675e+03 1.1718183638063283e+03
|
||||
27 3.9544661983654925e+02 -2.5857342971810368e+02 -3.0379689818002142e+01
|
||||
28 -1.7960213685989305e+02 3.0633168671228430e+02 3.8807838388686309e+02
|
||||
29 1.0106463589881000e+03 1.0687745018889480e+03 -6.6861303703586529e+01
|
||||
30 -6.0852844150139362e+01 -1.2007219497990148e+02 -5.4687005523315872e+02
|
||||
31 2.6341847548417515e+01 -6.0380405513895437e+02 3.8656694437683996e+02
|
||||
32 -1.8149728033842284e+03 -9.8411584459831852e+02 -2.8109959690485834e+01
|
||||
...
|
||||
18
unittest/force-styles/tests/spce.sw
Normal file
18
unittest/force-styles/tests/spce.sw
Normal file
@ -0,0 +1,18 @@
|
||||
type
|
||||
type
|
||||
type
|
||||
1 #epsilon in kcal/mol
|
||||
1 #sigma in dimensionless
|
||||
3.7 # a in Ang
|
||||
1.0 #lambda dimensionless
|
||||
0.8 #gamma in Ang
|
||||
0.0 #costheta0 dimensionless
|
||||
0 #two body part A=0
|
||||
0 #two body part B=0
|
||||
0 #two body part p=0
|
||||
0 #two body part q=0
|
||||
0.0 # use the standard Stillinger-Weber cutoff
|
||||
table_CG_CG_CG.txt
|
||||
VOTCA
|
||||
linear
|
||||
1001
|
||||
71
unittest/force-styles/tests/spce2.3b
Normal file
71
unittest/force-styles/tests/spce2.3b
Normal file
@ -0,0 +1,71 @@
|
||||
type1
|
||||
type1
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type1
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type2
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type1
|
||||
type2
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type1
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type1
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type2
|
||||
type1
|
||||
3.7 # cut in Ang
|
||||
1-1-2.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
type2
|
||||
type2
|
||||
type2
|
||||
3.7 # cut in Ang
|
||||
1-1-1.table
|
||||
ENTRY1
|
||||
linear
|
||||
12
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
1203
unittest/force-styles/tests/table_CG_CG.txt
Normal file
1203
unittest/force-styles/tests/table_CG_CG.txt
Normal file
File diff suppressed because it is too large
Load Diff
1004
unittest/force-styles/tests/table_CG_CG_CG.txt
Normal file
1004
unittest/force-styles/tests/table_CG_CG_CG.txt
Normal file
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user