lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3468 from akohlmey/collected-small-changes 

Collected small changes
This commit is contained in:
Axel Kohlmeyer 2022-10-03 09:37:29 -04:00 committed by GitHub
parent f136bcbef7 d7ff6337ed
commit 2970a2140f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
40 changed files with 980 additions and 546 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
doc/src/Build_manual.rst
View File

@ -216,9 +216,20 @@ be multiple tests run automatically:
- A test that only standard, printable ASCII text characters are used. - A test that only standard, printable ASCII text characters are used.
This runs the command ``env LC_ALL=C grep -n '[^ -~]' src/*.rst`` and This runs the command ``env LC_ALL=C grep -n '[^ -~]' src/*.rst`` and
thus prints all offending lines with filename and line number thus prints all offending lines with filename and line number
prepended to the screen. Special characters like the Angstrom prepended to the screen. Special characters like greek letters
:math:`\mathrm{\mathring{A}}` should be typeset with embedded math (:math:`\alpha~~\sigma~~\epsilon`), super- or subscripts
(like this ``:math:`\mathrm{\mathring{A}}```\ ). (:math:`x^2~~\mathrm{U}_{LJ}`), mathematical expressions
(:math:`\frac{1}{2}\mathrm{N}~~x\to\infty`), or the Angstrom symbol
(:math:`\AA`) should be typeset with embedded LaTeX (like this
``:math:`\alpha \sigma \epsilon```, ``:math:`x^2 \mathrm{E}_{LJ}```,
``:math:`\frac{1}{2}\mathrm{N} x\to\infty```, or ``:math:`\AA```\ ).
- Embedded LaTeX is rendered in HTML output with `MathJax
<https://www.mathjax.org/>`_ and in PDF output by passing the embedded
text to LaTeX. Some care has to be taken, though, since there are
limitations which macros and features can be used in either mode, so
it is recommended to always check whether any new or changed
documentation does translate and render correctly with either output.
- A test whether all styles are documented and listed in their - A test whether all styles are documented and listed in their
respective overview pages. A typical output with warnings looks like this: respective overview pages. A typical output with warnings looks like this:

24
doc/src/Speed_packages.rst
View File

@ -117,33 +117,15 @@ script.
with all its accelerator packages installed. Note however that the with all its accelerator packages installed. Note however that the
INTEL and KOKKOS packages require you to choose one of their INTEL and KOKKOS packages require you to choose one of their
hardware options when building for a specific platform. I.e. CPU or hardware options when building for a specific platform. I.e. CPU or
Phi option for the INTEL package. Or the OpenMP, Cuda, or Phi Phi option for the INTEL package. Or the OpenMP, CUDA, HIP, SYCL,
option for the KOKKOS package. or Phi option for the KOKKOS package. Or the OpenCL, HIP, or CUDA
option for the GPU package.
These are the exceptions. You cannot build a single executable with: These are the exceptions. You cannot build a single executable with:
* both the INTEL Phi and KOKKOS Phi options * both the INTEL Phi and KOKKOS Phi options
* the INTEL Phi or Kokkos Phi option, and the GPU package * the INTEL Phi or Kokkos Phi option, and the GPU package
See the examples/accelerate/README and make.list files for sample
Make.py commands that build LAMMPS with any or all of the accelerator
packages. As an example, here is a command that builds with all the
GPU related packages installed (GPU, KOKKOS with Cuda), including
settings to build the needed auxiliary GPU libraries for Kepler GPUs:
.. code-block:: bash
Make.py -j 16 -p omp gpu kokkos -cc nvcc wrap=mpi -gpu mode=double arch=35 -kokkos cuda arch=35 lib-all file mpi
The examples/accelerate directory also has input scripts that can be
used with all of the accelerator packages. See its README file for
details.
Likewise, the bench directory has FERMI and KEPLER and PHI
sub-directories with Make.py commands and input scripts for using all
the accelerator packages on various machines. See the README files in
those directories.
As mentioned above, the `Benchmark page <https://www.lammps.org/bench.html>`_ of the LAMMPS website gives As mentioned above, the `Benchmark page <https://www.lammps.org/bench.html>`_ of the LAMMPS website gives
performance results for the various accelerator packages for several performance results for the various accelerator packages for several
of the standard LAMMPS benchmark problems, as a function of problem of the standard LAMMPS benchmark problems, as a function of problem

10
doc/src/bond_mm3.rst
View File

@ -26,14 +26,14 @@ as defined in :ref:`(Allinger) <mm3-allinger1989>`
.. math:: .. math::
E = K (r - r_0)^2 \left[ 1 - 2.55(r-r_0) + (7/12) 2.55^2(r-r_0)^2 \right] E = K (r - r_0)^2 \left[ 1 - 2.55(r-r_0) + \frac{7}{12} 2.55^2(r-r_0)^2 \right]
where :math:`r_0` is the equilibrium value of the bond, and :math:`K` is a where :math:`r_0` is the equilibrium value of the bond, and :math:`K` is a
prefactor. The anharmonic prefactors have units angstrom\^(-n): prefactor. The anharmonic prefactors have units :math:`\AA^{-n}`:
-2.55 angstrom\^(-1) and (7/12)2.55\^2 angstrom\^(-2). The code takes :math:`-2.55 \AA^{-1}` and :math:`\frac{7}{12} 2.55^2 \AA^{-2}`. The code takes
care of the necessary unit conversion for these factors internally. care of the necessary unit conversion for these factors internally.
Note that the MM3 papers contains an error in Eq (1): Note that the MM3 papers contain an error in Eq (1):
(7/12)2.55 should be replaced with (7/12)2.55\^2 :math:`\frac{7}{12} 2.55` should be replaced with :math:`\frac{7}{12} 2.55^2`
The following coefficients must be defined for each bond type via the The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in :doc:`bond_coeff <bond_coeff>` command as in the example above, or in

19
doc/src/bond_quartic.rst
View File

@ -28,11 +28,18 @@ The *quartic* bond style uses the potential
.. math:: .. math::
E = K (r - R_c)^ 2 (r - R_c - B_1) (r - R_c - B_2) + U_0 + 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right] + \epsilon E & = E_q + E_{LJ} \\
E_q & = K (r - R_c)^ 2 (r - R_c - B_1) (r - R_c - B_2) + U_0 \\
E_{LJ} & = \left\{ \begin{array} {l@{\quad:\quad}l}
4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right] + \epsilon & r < 2^{\frac{1}{6}}, \epsilon = 1, \sigma = 1 \\
0 & r >= 2^{\frac{1}{6}}
\end{array} \right.
to define a bond that can be broken as the simulation proceeds (e.g. to define a bond that can be broken as the simulation proceeds (e.g.
due to a polymer being stretched). The :math:`\sigma` and :math:`\epsilon` used in the due to a polymer being stretched). The :math:`\sigma` and
LJ portion of the formula are both set equal to 1.0 by LAMMPS. :math:`\epsilon` used in the LJ portion of the formula are both set
equal to 1.0 by LAMMPS and the LJ portion is cut off at its minimum,
i.e. at :math:`r_c = 2^{\frac{1}{6}}`.
The following coefficients must be defined for each bond type via the The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in :doc:`bond_coeff <bond_coeff>` command as in the example above, or in
@ -46,9 +53,9 @@ or :doc:`read_restart <read_restart>` commands:
* :math:`U_0` (energy) * :math:`U_0` (energy)
This potential was constructed to mimic the FENE bond potential for This potential was constructed to mimic the FENE bond potential for
coarse-grained polymer chains. When monomers with :math:`\sigma = \epsilon = 1.0` coarse-grained polymer chains. When monomers with :math:`\sigma =
are used, the following choice of parameters gives a quartic potential that \epsilon = 1.0` are used, the following choice of parameters gives a
looks nearly like the FENE potential: quartic potential that looks nearly like the FENE potential:
.. math:: .. math::

3
doc/src/compute_chunk_atom.rst
View File

@ -602,8 +602,7 @@ be used. For non-orthogonal (triclinic) simulation boxes, only the
*reduced* option may be used. *reduced* option may be used.
A *box* value selects standard distance units as defined by the A *box* value selects standard distance units as defined by the
:doc:`units <units>` command (e.g., :math:`\mathrm{\mathring A}` :doc:`units <units>` command (e.g., :math:`\AA` for units = *real* or *metal*).
for units = *real* or *metal*).
A *lattice* value means the distance units are in lattice spacings. A *lattice* value means the distance units are in lattice spacings.
The :doc:`lattice <lattice>` command must have been previously used to The :doc:`lattice <lattice>` command must have been previously used to
define the lattice spacing. A *reduced* value means normalized define the lattice spacing. A *reduced* value means normalized

2
doc/src/compute_displace_atom.rst
View File

@ -95,7 +95,7 @@ something like the following commands:
refresh c_dsp delay 100 refresh c_dsp delay 100
The :doc:`dump_modify thresh <dump_modify>` command will only output The :doc:`dump_modify thresh <dump_modify>` command will only output
atoms that have displaced more than :math:`0.6~\mathrm{\mathring A}` on each atoms that have displaced more than :math:`0.6~\AA` on each
snapshot (assuming metal units). The dump_modify *refresh* option triggers a snapshot (assuming metal units). The dump_modify *refresh* option triggers a
call to this compute at the end of every dump. call to this compute at the end of every dump.

4
doc/src/compute_entropy_atom.rst
View File

@ -97,13 +97,13 @@ by the corresponding volume. This option can be useful when dealing with
inhomogeneous systems such as those that have surfaces. inhomogeneous systems such as those that have surfaces.
Here are typical input parameters for fcc aluminum (lattice Here are typical input parameters for fcc aluminum (lattice
constant :math:`4.05~\mathrm{\mathring A}`), constant :math:`4.05~\AA`),
.. parsed-literal:: .. parsed-literal::
compute 1 all entropy/atom 0.25 5.7 avg yes 3.7 compute 1 all entropy/atom 0.25 5.7 avg yes 3.7
and for bcc sodium (lattice constant 4.23 Angstroms), and for bcc sodium (lattice constant :math:`4.23~\AA`),
.. parsed-literal:: .. parsed-literal::

3
doc/src/compute_rigid_local.rst
View File

@ -109,8 +109,7 @@ The *mass* attribute is the total mass of the rigid body.
There are two options for outputting the coordinates of the center of There are two options for outputting the coordinates of the center of
mass (COM) of the body. The *x*, *y*, *z* attributes write the COM mass (COM) of the body. The *x*, *y*, *z* attributes write the COM
"unscaled", in the appropriate distance :doc:`units <units>` "unscaled", in the appropriate distance :doc:`units <units>`
(:math:`\mathrm{\mathring A}`, (:math:`\AA`, :math:`\sigma`, etc). Use *xu*, *yu*, *zu* if you want the COM "unwrapped" by
sigma, etc). Use *xu*, *yu*, *zu* if you want the COM "unwrapped" by
the image flags for each body. Unwrapped means that if the body the image flags for each body. Unwrapped means that if the body
COM has passed through a periodic boundary one or more times, the value COM has passed through a periodic boundary one or more times, the value
is generated what the COM coordinate would be if it had not been is generated what the COM coordinate would be if it had not been

2
doc/src/compute_saed.rst
View File

@ -86,7 +86,7 @@ will defined using the *c* values for the spacing along each reciprocal
lattice axis. Note that manual mapping of the reciprocal space mesh is lattice axis. Note that manual mapping of the reciprocal space mesh is
good for comparing diffraction results from multiple simulations; however good for comparing diffraction results from multiple simulations; however
it can reduce the likelihood that Bragg reflections will be satisfied it can reduce the likelihood that Bragg reflections will be satisfied
unless small spacing parameters (:math:`<0.05~\mathrm{\mathring A}^-1`) unless small spacing parameters (:math:`<0.05~\AA^-1`)
are implemented. Meshes with manual spacing do not require a periodic are implemented. Meshes with manual spacing do not require a periodic
boundary. boundary.

2
doc/src/compute_temp_ramp.rst
View File

@ -58,7 +58,7 @@ constant, and :math:`T` is the absolute temperature.
The *units* keyword determines the meaning of the distance units used The *units* keyword determines the meaning of the distance units used
for coordinates (*clo*, *chi*) and velocities (*vlo*, *vhi*). A *box* value for coordinates (*clo*, *chi*) and velocities (*vlo*, *vhi*). A *box* value
selects standard distance units as defined by the :doc:`units <units>` selects standard distance units as defined by the :doc:`units <units>`
command (e.g., :math:`\mathrm{\mathring{A}}` for units = real or metal). A command (e.g., :math:`\AA` for units = real or metal). A
*lattice* value means the distance units are in lattice spacings (i.e., *lattice* value means the distance units are in lattice spacings (i.e.,
velocity in lattice spacings per unit time). The :doc:`lattice <lattice>` velocity in lattice spacings per unit time). The :doc:`lattice <lattice>`
command must have been previously used to define the lattice spacing. command must have been previously used to define the lattice spacing.

2
doc/src/compute_xrd.rst
View File

@ -91,7 +91,7 @@ reciprocal lattice axis. Note that manual mapping of the reciprocal
space mesh is good for comparing diffraction results from multiple space mesh is good for comparing diffraction results from multiple
simulations; however, it can reduce the likelihood that Bragg simulations; however, it can reduce the likelihood that Bragg
reflections will be satisfied unless small spacing parameters reflections will be satisfied unless small spacing parameters
(:math:`< 0.05~\mathrm{\mathring{A}}^{-1}`) are implemented. (:math:`< 0.05~\AA^{-1}`) are implemented.
Meshes with manual spacing do not require a periodic boundary. Meshes with manual spacing do not require a periodic boundary.
The limits of the reciprocal lattice mesh are determined by range of The limits of the reciprocal lattice mesh are determined by range of

2
doc/src/create_atoms.rst
View File

@ -464,7 +464,7 @@ The *units* keyword determines the meaning of the distance units used
to specify the coordinates of the one particle created by the *single* to specify the coordinates of the one particle created by the *single*
style, or the overlap distance *Doverlap* by the *overlap* keyword. A style, or the overlap distance *Doverlap* by the *overlap* keyword. A
*box* value selects standard distance units as defined by the *box* value selects standard distance units as defined by the
:doc:`units <units>` command (e.g., :math:`\mathrm{\mathring{A}}` for :doc:`units <units>` command (e.g., :math:`\AA` for
units = *real* or *metal*\ . A *lattice* value means the distance units are in units = *real* or *metal*\ . A *lattice* value means the distance units are in
lattice spacings. lattice spacings.

2
doc/src/displace_atoms.rst
View File

@ -104,7 +104,7 @@ atom's rotation.
Distance units for displacements and the origin point of the *rotate* Distance units for displacements and the origin point of the *rotate*
style are determined by the setting of *box* or *lattice* for the style are determined by the setting of *box* or *lattice* for the
*units* keyword. *Box* means distance units as defined by the *units* keyword. *Box* means distance units as defined by the
:doc:`units <units>` command (e.g., :math:`\mathrm{\mathring A}` for :doc:`units <units>` command (e.g., :math:`\AA` for
*real* or *metal* units). *Lattice* means distance units are in lattice *real* or *metal* units). *Lattice* means distance units are in lattice
spacings. The :doc:`lattice <lattice>` command must have been previously used spacings. The :doc:`lattice <lattice>` command must have been previously used
to define the lattice spacing. to define the lattice spacing.

2
doc/src/dump.rst
View File

@ -693,7 +693,7 @@ charge.
There are several options for outputting atom coordinates. The *x*, There are several options for outputting atom coordinates. The *x*,
*y*, and *z* attributes write atom coordinates "unscaled", in the *y*, and *z* attributes write atom coordinates "unscaled", in the
appropriate distance :doc:`units <units>` (:math:`\mathrm{\mathring A}`, appropriate distance :doc:`units <units>` (:math:`\AA`,
:math:`\sigma`, etc.). Use *xs*, *ys*, and *zs* if you want the :math:`\sigma`, etc.). Use *xs*, *ys*, and *zs* if you want the
coordinates "scaled" to the box size so that each value is 0.0 to 1.0. coordinates "scaled" to the box size so that each value is 0.0 to 1.0.
If the simulation box is triclinic (tilted), then all atom coords will If the simulation box is triclinic (tilted), then all atom coords will

6
doc/src/dump_modify.rst
View File

@ -632,7 +632,7 @@ calculates the displacement of each atom from its reference position.
The "4" index is the scalar displacement; 1, 2, and 3 are the :math:`xyz` The "4" index is the scalar displacement; 1, 2, and 3 are the :math:`xyz`
components of the displacement. The :doc:`dump_modify thresh <dump_modify>` components of the displacement. The :doc:`dump_modify thresh <dump_modify>`
command will cause only atoms that have displaced more than command will cause only atoms that have displaced more than
:math:`0.6~\mathrm{\mathring A}` to be output on a given snapshot (assuming :math:`0.6~\AA` to be output on a given snapshot (assuming
metal units). However, note that when an atom is output, we also need to metal units). However, note that when an atom is output, we also need to
update the reference position for that atom to its new coordinates. So that it update the reference position for that atom to its new coordinates. So that it
will not be output in every snapshot thereafter. That reference position is will not be output in every snapshot thereafter. That reference position is
@ -675,7 +675,7 @@ value of *yes* means atom coords are written in normalized units from
0.0 to 1.0 in each box dimension. If the simulation box is triclinic 0.0 to 1.0 in each box dimension. If the simulation box is triclinic
(tilted), then all atom coords will still be between 0.0 and 1.0. A (tilted), then all atom coords will still be between 0.0 and 1.0. A
value of *no* means they are written in absolute distance units value of *no* means they are written in absolute distance units
(e.g., :math:`\mathrm{\mathring A}` or :math:`\sigma`). (e.g., :math:`\AA` or :math:`\sigma`).
Using this keyword will reset all custom header names set with Using this keyword will reset all custom header names set with
*dump_modify colname* to their respective default values. *dump_modify colname* to their respective default values.
@ -687,7 +687,7 @@ when writing to XTC files. By default, they are initialized for
whatever :doc:`units <units>` style is being used, to write out whatever :doc:`units <units>` style is being used, to write out
coordinates in nanometers and time in picoseconds. For example, for *real* coordinates in nanometers and time in picoseconds. For example, for *real*
units, LAMMPS defines *sfactor* = 0.1 and *tfactor* = 0.001, since the units, LAMMPS defines *sfactor* = 0.1 and *tfactor* = 0.001, since the
:math:`\mathrm{\mathring A}` and fs used by *real* units are 0.1 nm and :math:`\AA` and fs used by *real* units are 0.1 nm and
0.001 ps, respectively. If you are using a units system with distance and time 0.001 ps, respectively. If you are using a units system with distance and time
units far from nm and ps, you may wish to write XTC files with units far from nm and ps, you may wish to write XTC files with
different units, since the compression algorithm used in XTC files is different units, since the compression algorithm used in XTC files is

4
doc/src/fix_acks2_reaxff.rst
View File

@ -71,7 +71,7 @@ potential in eV, *gamma*, the valence orbital exponent, and *bcut*, the
bond cutoff distance. Note that these 4 quantities are also in the bond cutoff distance. Note that these 4 quantities are also in the
ReaxFF potential file, except that eta is defined here as twice the eta ReaxFF potential file, except that eta is defined here as twice the eta
value in the ReaxFF file. Note that unlike the rest of LAMMPS, the units value in the ReaxFF file. Note that unlike the rest of LAMMPS, the units
of this fix are hard-coded to be :math:`\mathrm{\mathring{A}}`, eV, and of this fix are hard-coded to be :math:`\AA`, eV, and
electronic charge. electronic charge.
The optional *maxiter* keyword allows changing the max number The optional *maxiter* keyword allows changing the max number
@ -111,7 +111,7 @@ LAMMPS was built with that package. See the :doc:`Build package
This fix does not correctly handle interactions involving multiple This fix does not correctly handle interactions involving multiple
periodic images of the same atom. Hence, it should not be used for periodic images of the same atom. Hence, it should not be used for
periodic cell dimensions less than :math:`10~\mathrm{\mathring{A}}`. periodic cell dimensions less than :math:`10~\AA`.
This fix may be used in combination with :doc:`fix efield <fix_efield>` This fix may be used in combination with :doc:`fix efield <fix_efield>`
and will apply the external electric field during charge equilibration, and will apply the external electric field during charge equilibration,

2
doc/src/fix_append_atoms.rst
View File

@ -79,7 +79,7 @@ measured from zhi and is set with the *extent* argument.
The *units* keyword determines the meaning of the distance units used The *units* keyword determines the meaning of the distance units used
to define a wall position, but only when a numeric constant is used. to define a wall position, but only when a numeric constant is used.
A *box* value selects standard distance units as defined by the A *box* value selects standard distance units as defined by the
:doc:`units <units>` command (e.g., :math:`\mathrm{\mathring A}` :doc:`units <units>` command (e.g., :math:`\AA`
for units = real or metal. for units = real or metal.
A *lattice* value means the distance units are in lattice spacings. A *lattice* value means the distance units are in lattice spacings.
The :doc:`lattice <lattice>` command must have been previously used to The :doc:`lattice <lattice>` command must have been previously used to

6
doc/src/fix_bocs.rst
View File

@ -59,7 +59,7 @@ Note both the COMMA and the SPACE separating the volume's
value and its corresponding pressure correction. The volumes in the file value and its corresponding pressure correction. The volumes in the file
must be uniformly spaced. Both the volumes and the pressure corrections must be uniformly spaced. Both the volumes and the pressure corrections
should be provided in the proper units, e.g. if you are using *units real*, should be provided in the proper units, e.g. if you are using *units real*,
the volumes should all be in cubic angstroms, and the pressure corrections the volumes should all be in cubic Angstroms, and the pressure corrections
should all be in atmospheres. Furthermore, the table should start/end at a should all be in atmospheres. Furthermore, the table should start/end at a
volume considerably smaller/larger than you expect your system to sample volume considerably smaller/larger than you expect your system to sample
during the simulation. If the system ever reaches a volume outside of the during the simulation. If the system ever reaches a volume outside of the
@ -72,8 +72,8 @@ With the *analytic* option, the arguments are as follows:
... analytic V_avg N_particles N_coeff Coeff_1 Coeff_2 ... Coeff_N ... analytic V_avg N_particles N_coeff Coeff_1 Coeff_2 ... Coeff_N
Note that *V_avg* and *Coeff_i* should all be in the proper units, e.g. if you Note that *V_avg* and *Coeff_i* should all be in the proper units, e.g. if you
are using *units real*, *V_avg* should be in cubic angstroms, and the are using *units real*, *V_avg* should be in cubic Angstroms, and the
coefficients should all be in atmospheres \* cubic angstroms. coefficients should all be in atmospheres \* cubic Angstroms.
---------- ----------

4
doc/src/fix_imd.rst
View File

@ -88,8 +88,8 @@ to send "unwrapped" coordinates to the IMD client that undo the
wrapping back of coordinates into the principle unit cell, as done by wrapping back of coordinates into the principle unit cell, as done by
default in LAMMPS. The *fscale* keyword allows to apply a scaling default in LAMMPS. The *fscale* keyword allows to apply a scaling
factor to forces transmitted by the IMD client. The IMD protocols factor to forces transmitted by the IMD client. The IMD protocols
stipulates that forces are transferred in kcal/mol/angstrom under the stipulates that forces are transferred in kcal/mol/Angstrom under the
assumption that coordinates are given in angstrom. For LAMMPS runs assumption that coordinates are given in Angstrom. For LAMMPS runs
with different units or as a measure to tweak the forces generated by with different units or as a measure to tweak the forces generated by
the manipulation of the IMD client, this option allows to make the manipulation of the IMD client, this option allows to make
adjustments. adjustments.

2
doc/src/fix_qeq_reaxff.rst
View File

@ -124,7 +124,7 @@ LAMMPS was built with that package. See the :doc:`Build package
This fix does not correctly handle interactions involving multiple This fix does not correctly handle interactions involving multiple
periodic images of the same atom. Hence, it should not be used for periodic images of the same atom. Hence, it should not be used for
periodic cell dimensions less than 10 angstroms. periodic cell dimensions less than 10 Angstroms.
This fix may be used in combination with :doc:`fix efield <fix_efield>` This fix may be used in combination with :doc:`fix efield <fix_efield>`
and will apply the external electric field during charge equilibration, and will apply the external electric field during charge equilibration,

2
doc/src/kspace_style.rst
View File

@ -301,7 +301,7 @@ and for mixed periodic and non-periodic boundaries.
MSM is most competitive versus Ewald and PPPM when only relatively MSM is most competitive versus Ewald and PPPM when only relatively
low accuracy forces, about 1e-4 relative error or less accurate, low accuracy forces, about 1e-4 relative error or less accurate,
are needed. Note that use of a larger Coulombic cutoff (i.e. 15 are needed. Note that use of a larger Coulombic cutoff (i.e. 15
angstroms instead of 10 angstroms) provides better MSM accuracy for Angstroms instead of 10 Angstroms) provides better MSM accuracy for
both the real space and grid computed forces. both the real space and grid computed forces.
Currently calculation of the full pressure tensor in MSM is expensive. Currently calculation of the full pressure tensor in MSM is expensive.

6
doc/src/pair_airebo.rst
View File

@ -102,7 +102,7 @@ a few fitted spline values are slightly different. For most cases the
statistical averages as the original REBO potential from which it was statistical averages as the original REBO potential from which it was
derived. The :math:`E^{\text{REBO}}` term in the AIREBO potential gives the model its derived. The :math:`E^{\text{REBO}}` term in the AIREBO potential gives the model its
reactive capabilities and only describes short-ranged C-C, C-H and H-H reactive capabilities and only describes short-ranged C-C, C-H and H-H
interactions (:math:`r < 2` Angstroms). These interactions have strong interactions (:math:`r < 2 \AA`). These interactions have strong
coordination-dependence through a bond order parameter, which adjusts coordination-dependence through a bond order parameter, which adjusts
the attraction between the I,J atoms based on the position of other the attraction between the I,J atoms based on the position of other
nearby atoms and thus has 3- and 4-body dependence. nearby atoms and thus has 3- and 4-body dependence.
@ -116,9 +116,9 @@ interactions is determined by the *cutoff* argument to the pair_style
command which is a scale factor. For each type pair (C-C, C-H, H-H) command which is a scale factor. For each type pair (C-C, C-H, H-H)
the cutoff is obtained by multiplying the scale factor by the sigma the cutoff is obtained by multiplying the scale factor by the sigma
value defined in the potential file for that type pair. In the value defined in the potential file for that type pair. In the
standard AIREBO potential, :math:`\sigma_{CC} = 3.4` Angstroms, so with a scale standard AIREBO potential, :math:`\sigma_{CC} = 3.4 \AA`, so with a scale
factor of 3.0 (the argument in pair_style), the resulting :math:`E^{\text{LJ}}` cutoff factor of 3.0 (the argument in pair_style), the resulting :math:`E^{\text{LJ}}` cutoff
would be 10.2 Angstroms. would be :math:`10.2 \AA`.
By default, the longer-ranged interaction is smoothly switched off By default, the longer-ranged interaction is smoothly switched off
between 2.16 and 3.0 :math:`\sigma`. By specifying *cutoff_min* in addition between 2.16 and 3.0 :math:`\sigma`. By specifying *cutoff_min* in addition

2
doc/src/pair_charmm.rst
View File

@ -144,7 +144,7 @@ artifacts.
conversion factor used internally in the code, from the LAMMPS value conversion factor used internally in the code, from the LAMMPS value
to the CHARMM value, as if it were effectively a parameter of the to the CHARMM value, as if it were effectively a parameter of the
force field. This is because the CHARMM code uses a slightly force field. This is because the CHARMM code uses a slightly
different value for the this conversion factor in :doc:`real units <units>` (Kcal/mole), namely CHARMM = 332.0716, LAMMPS = different value for the this conversion factor in :doc:`real units <units>` (kcal/mol), namely CHARMM = 332.0716, LAMMPS =
332.06371. This is to enable more precise agreement by LAMMPS with 332.06371. This is to enable more precise agreement by LAMMPS with
the CHARMM force field energies and forces, when using one of these the CHARMM force field energies and forces, when using one of these
two CHARMM pair styles. two CHARMM pair styles.

16
doc/src/pair_smtbq.rst
View File

@ -182,19 +182,19 @@ For each cations (metal):
* Potential parameter: * Potential parameter:
- If type of potential is 'second_moment' : A (eV), *p*, - If type of potential is 'second_moment' : A (eV), *p*,
:math:`\zeta^0` (eV) and *q*, :math:`r_{c1} (\mathrm{\mathring{A}})`, :math:`r_{c2} :math:`\zeta^0` (eV) and *q*, :math:`r_{c1} (\AA)`, :math:`r_{c2}
(\mathrm{\mathring{A}})` and :math:`r_0 (\mathrm{\mathring{A}})` (\AA)` and :math:`r_0 (\AA)`
- If type of potential is 'buck' : *C* (eV) and :math:`\rho (\mathrm{\mathring{A}})` - If type of potential is 'buck' : *C* (eV) and :math:`\rho (\AA)`
- If type of potential is 'buckPlusAttr' : *C* (eV) and :math:`\rho - If type of potential is 'buckPlusAttr' : *C* (eV) and :math:`\rho
(\mathrm{\mathring{A}})` *D* (eV), *B* :math:`(\mathrm{\mathring{A}}^{-1})`, :math:`r^{OO}_1 (\mathrm{\mathring{A}})` and (\AA)` *D* (eV), *B* :math:`(\AA^{-1})`, :math:`r^{OO}_1 (\AA)` and
:math:`r^{OO}_2 (\mathrm{\mathring{A}})` :math:`r^{OO}_2 (\AA)`
* Divider line * Divider line
4) Tables parameters: 4) Tables parameters:
* Cutoff radius for the Coulomb interaction (:math:`R_{coul}`) * Cutoff radius for the Coulomb interaction (:math:`R_{coul}`)
* Starting radius (:math:`r_{min} = 1,18845 \mathrm{\mathring{A}}`) and increments * Starting radius (:math:`r_{min} = 1,18845 \AA`) and increments
(:math:`dr = 0.001 \mathrm{\mathring{A}}`) for creating the potential table. (:math:`dr = 0.001 \AA`) for creating the potential table.
* Divider line * Divider line
5) Rick model parameter: 5) Rick model parameter:
@ -208,7 +208,7 @@ For each cations (metal):
6) Coordination parameter: 6) Coordination parameter:
* First (:math:`r_{1n}`) and second (:math:`r_{2n}`) neighbor distances * First (:math:`r_{1n}`) and second (:math:`r_{2n}`) neighbor distances
in angstrom in Angstrom
* Divider line * Divider line
7) Charge initialization mode: 7) Charge initialization mode:

2
doc/src/pair_zbl.rst
View File

@ -73,7 +73,7 @@ be included in a pair_coeff command.
The numerical values of the exponential decay constants in the The numerical values of the exponential decay constants in the
screening function depend on the unit of distance. In the above screening function depend on the unit of distance. In the above
equation they are given for units of angstroms. LAMMPS will equation they are given for units of Angstroms. LAMMPS will
automatically convert these values to the distance unit of the automatically convert these values to the distance unit of the
specified LAMMPS :doc:`units <units>` setting. The values of Z should specified LAMMPS :doc:`units <units>` setting. The values of Z should
always be given as multiples of a proton's charge, e.g. 29.0 for always be given as multiples of a proton's charge, e.g. 29.0 for

4
doc/src/run_style.rst
View File

@ -221,7 +221,7 @@ impropers, and dihedrals can be computed on this innermost 0.5 fs
step. The outermost timestep cannot be greater than 4.0 fs without step. The outermost timestep cannot be greater than 4.0 fs without
risking energy drift. Smooth switching of forces between the levels risking energy drift. Smooth switching of forces between the levels
of the rRESPA hierarchy is also necessary to avoid drift, and a 1-2 of the rRESPA hierarchy is also necessary to avoid drift, and a 1-2
angstrom "healing distance" (the distance between the outer and inner Angstrom "healing distance" (the distance between the outer and inner
cutoffs) works reasonably well. We thus recommend the following cutoffs) works reasonably well. We thus recommend the following
settings for use of the *respa* style without SHAKE in biomolecular settings for use of the *respa* style without SHAKE in biomolecular
simulations: simulations:
@ -277,7 +277,7 @@ Even a LJ system can benefit from rRESPA if the interactions are
divided by the inner, middle and outer keywords. A 2-fold or more divided by the inner, middle and outer keywords. A 2-fold or more
speedup can be obtained while maintaining good energy conservation. speedup can be obtained while maintaining good energy conservation.
In real units, for a pure LJ fluid at liquid density, with a sigma of In real units, for a pure LJ fluid at liquid density, with a sigma of
3.0 angstroms, and epsilon of 0.1 Kcal/mol, the following settings 3.0 Angstroms, and epsilon of 0.1 kcal/mol, the following settings
seem to work well: seem to work well:
.. code-block:: LAMMPS .. code-block:: LAMMPS

10
doc/src/units.rst
View File

@ -30,7 +30,7 @@ and dump files. Typically, this command is used at the very beginning
of an input script. of an input script.
For all units except *lj*, LAMMPS uses physical constants from For all units except *lj*, LAMMPS uses physical constants from
www.physics.nist.gov. For the definition of Kcal in real units, www.physics.nist.gov. For the definition of kcal in real units,
LAMMPS uses the thermochemical calorie = 4.184 J. LAMMPS uses the thermochemical calorie = 4.184 J.
The choice you make for units simply sets some internal conversion The choice you make for units simply sets some internal conversion
@ -102,17 +102,17 @@ For style *real*, these are the units:
* mass = grams/mole * mass = grams/mole
* distance = Angstroms * distance = Angstroms
* time = femtoseconds * time = femtoseconds
* energy = Kcal/mole * energy = kcal/mol
* velocity = Angstroms/femtosecond * velocity = Angstroms/femtosecond
* force = Kcal/mole-Angstrom * force = (kcal/mol)/Angstrom
* torque = Kcal/mole * torque = kcal/mol
* temperature = Kelvin * temperature = Kelvin
* pressure = atmospheres * pressure = atmospheres
* dynamic viscosity = Poise * dynamic viscosity = Poise
* charge = multiple of electron charge (1.0 is a proton) * charge = multiple of electron charge (1.0 is a proton)
* dipole = charge\*Angstroms * dipole = charge\*Angstroms
* electric field = volts/Angstrom * electric field = volts/Angstrom
* density = gram/cm\^dim * density = g/cm\^dim
For style *metal*, these are the units: For style *metal*, these are the units:

11
doc/src/variable.rst
View File

@ -319,11 +319,12 @@ in the input script, or if the script is read again in a loop. The other
difference is that *string* performs variable substitution even if the difference is that *string* performs variable substitution even if the
string parameter is quoted. string parameter is quoted.
For the *format* style, an equal-style variable is specified along For the *format* style, an equal-style or compatible variable is
with a C-style format string, e.g. "%f" or "%.10g", which must be specified along with a C-style format string, e.g. "%f" or "%.10g",
appropriate for formatting a double-precision floating-point value. which must be appropriate for formatting a double-precision
The default format is "%.15g". This variable style allows an floating-point value and may not have extra characters. The default
equal-style variable to be formatted precisely when it is evaluated. format is "%.15g". This variable style allows an equal-style variable
to be formatted precisely when it is evaluated.
Note that if you simply wish to print a variable value with desired Note that if you simply wish to print a variable value with desired
precision to the screen or logfile via the :doc:`print <print>` or precision to the screen or logfile via the :doc:`print <print>` or

14
doc/utils/sphinx-config/conf.py.in
View File

@ -159,7 +159,6 @@ pygments_style = 'default'
# If true, keep warnings as "system message" paragraphs in the built documents. # If true, keep warnings as "system message" paragraphs in the built documents.
#keep_warnings = False #keep_warnings = False
# -- Options for HTML output ---------------------------------------------- # -- Options for HTML output ----------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for # The theme to use for HTML and HTML Help pages. See the documentation for
@ -264,8 +263,18 @@ else:
html_math_renderer = 'mathjax' html_math_renderer = 'mathjax'
# use relative path for mathjax, so it is looked for in the # use relative path for mathjax, so it is looked for in the
# html tree and the manual becomes readable when offline # html tree and the manual becomes readable when offline
mathjax_path = 'mathjax/es5/tex-mml-chtml.js' mathjax_path = 'mathjax/es5/tex-mml-chtml.js'
# hack to enable use of \AA in :math:
rst_prolog = r"""
.. only:: html
:math:`\renewcommand{\AA}{\text{Å}}`
"""
# -- Options for LaTeX output --------------------------------------------- # -- Options for LaTeX output ---------------------------------------------
latex_elements = { latex_elements = {
@ -278,6 +287,7 @@ latex_elements = {
# Additional stuff for the LaTeX preamble. # Additional stuff for the LaTeX preamble.
'preamble': r''' 'preamble': r'''
\setcounter{tocdepth}{2} \setcounter{tocdepth}{2}
\renewcommand{\AA}{\mbox{\textrm{\r{A}}}}
''' '''
} }

6
src/compute_chunk_atom.cpp
View File

@ -906,7 +906,11 @@ void ComputeChunkAtom::assign_chunk_ids()
// update region if necessary // update region if necessary
if (regionflag) region->prematch(); if (regionflag) {
region = domain->get_region_by_id(idregion);
if (!region) error->all(FLERR, "Region {} for compute chunk/atom does not exist", idregion);
region->prematch();
}
// exclude = 1 if atom is not assigned to a chunk // exclude = 1 if atom is not assigned to a chunk
// exclude atoms not in group or not in optional region // exclude atoms not in group or not in optional region

156
src/compute_chunk_spread_atom.cpp
View File

@ -32,8 +32,7 @@ using namespace LAMMPS_NS;
ComputeChunkSpreadAtom:: ComputeChunkSpreadAtom::
ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) : ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), idchunk(nullptr)
idchunk(nullptr), ids(nullptr), which(nullptr), argindex(nullptr), value2index(nullptr)
{ {
if (narg < 5) error->all(FLERR,"Illegal compute chunk/spread/atom command"); if (narg < 5) error->all(FLERR,"Illegal compute chunk/spread/atom command");
@ -54,32 +53,27 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
// parse values // parse values
which = new int[nargnew]; values.clear();
argindex = new int[nargnew]; for (iarg = 0; iarg < nargnew; ++iarg) {
ids = new char*[nargnew];
value2index = new int[nargnew];
nvalues = 0;
for (iarg = 0; iarg < nargnew; iarg++) {
ids[nvalues] = nullptr;
ArgInfo argi(arg[iarg], ArgInfo::COMPUTE|ArgInfo::FIX); ArgInfo argi(arg[iarg], ArgInfo::COMPUTE|ArgInfo::FIX);
which[nvalues] = argi.get_type(); value_t val;
argindex[nvalues] = argi.get_index1(); val.which = argi.get_type();
ids[nvalues] = argi.copy_name(); val.argindex = argi.get_index1();
val.id = argi.get_name();
val.val.c = nullptr;
if ((which[nvalues] == ArgInfo::UNKNOWN) || (which[nvalues] == ArgInfo::NONE) if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE)
|| (argi.get_dim() > 1)) || (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute chunk/spread/atom command"); error->all(FLERR,"Illegal compute chunk/spread/atom argument: {}", arg[iarg]);
nvalues++; values.push_back(val);
} }
// if wildcard expansion occurred, free earg memory from expand_args() // if wildcard expansion occurred, free earg memory from expand_args()
if (expand) { if (expand) {
for (int i = 0; i < nargnew; i++) delete [] earg[i]; for (int i = 0; i < nargnew; i++) delete[] earg[i];
memory->sfree(earg); memory->sfree(earg);
} }
@ -87,38 +81,45 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
// for compute, must calculate per-chunk values, i.e. style ends in "/chunk" // for compute, must calculate per-chunk values, i.e. style ends in "/chunk"
// for fix, assume a global vector or array is per-chunk data // for fix, assume a global vector or array is per-chunk data
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
if (which[i] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
auto icompute = modify->get_compute_by_id(ids[i]); auto icompute = modify->get_compute_by_id(val.id);
if (!icompute) if (!icompute)
error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", ids[i]); error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", val.id);
if (!utils::strmatch(icompute->style,"/chunk$")) if (!utils::strmatch(icompute->style,"/chunk$"))
error->all(FLERR,"Compute for compute chunk/spread/atom " error->all(FLERR,"Compute chunk/spread/atom compute {} does not calculate per-chunk values",
"does not calculate per-chunk values"); val.id);
if (argindex[i] == 0) { if (val.argindex == 0) {
if (!icompute->vector_flag) if (!icompute->vector_flag)
error->all(FLERR,"Compute chunk/spread/atom compute does not calculate global vector"); error->all(FLERR,"Compute chunk/spread/atom compute {} does not calculate global vector",
val.id);
} else { } else {
if (!icompute->array_flag) if (!icompute->array_flag)
error->all(FLERR,"Compute chunk/spread/atom compute does not calculate global array"); error->all(FLERR,"Compute chunk/spread/atom compute {} does not calculate global array",
if (argindex[i] > icompute->size_array_cols) val.id);
error->all(FLERR,"Compute chunk/spread/atom compute array is accessed out-of-range"); if (val.argindex > icompute->size_array_cols)
error->all(FLERR,"Compute chunk/spread/atom compute {} array is accessed out-of-range",
val.id);
} }
val.val.c = icompute;
} else if (which[i] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]); auto ifix = modify->get_fix_by_id(val.id);
if (ifix) if (ifix)
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", ids[i]); error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", val.id);
if (argindex[i] == 0) { if (val.argindex == 0) {
if (!ifix->vector_flag) if (!ifix->vector_flag)
error->all(FLERR,"Compute chunk/spread/atom fix does not calculate global vector"); error->all(FLERR,"Compute chunk/spread/atom {} fix does not calculate global vector",
val.id);
} else { } else {
if (!ifix->array_flag) if (!ifix->array_flag)
error->all(FLERR,"Compute chunk/spread/atom fix does not calculate global array"); error->all(FLERR,"Compute chunk/spread/atom {} fix does not calculate global array",
if (argindex[i] > ifix->size_array_cols) val.id);
error->all(FLERR,"Compute chunk/spread/atom fix array is accessed out-of-range"); if (val.argindex > ifix->size_array_cols)
error->all(FLERR,"Compute chunk/spread/atom fix {} array is accessed out-of-range",
val.id);
} }
} }
} }
@ -126,8 +127,8 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
// this compute produces a peratom vector or array // this compute produces a peratom vector or array
peratom_flag = 1; peratom_flag = 1;
if (nvalues == 1) size_peratom_cols = 0; if (values.size() == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues; else size_peratom_cols = values.size();
// per-atom vector or array // per-atom vector or array
@ -140,13 +141,7 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
ComputeChunkSpreadAtom::~ComputeChunkSpreadAtom() ComputeChunkSpreadAtom::~ComputeChunkSpreadAtom()
{ {
delete [] idchunk; delete[] idchunk;
delete [] which;
delete [] argindex;
for (int i = 0; i < nvalues; i++) delete [] ids[i];
delete [] ids;
delete [] value2index;
memory->destroy(vector_atom); memory->destroy(vector_atom);
memory->destroy(array_atom); memory->destroy(array_atom);
@ -160,18 +155,16 @@ void ComputeChunkSpreadAtom::init()
// set indices of all computes,fixes,variables // set indices of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) { for (auto &val : values) {
if (which[m] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]); val.val.c = modify->get_compute_by_id(val.id);
if (icompute < 0) if (!val.val.c)
error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", ids[m]); error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", val.id);
value2index[m] = icompute;
} else if (which[m] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]); val.val.f = modify->get_fix_by_id(val.id);
if (ifix < 0) if (!val.val.f)
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", ids[m]); error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", val.id);
value2index[m] = ifix;
} }
} }
} }
@ -182,7 +175,8 @@ void ComputeChunkSpreadAtom::init_chunk()
{ {
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk)); cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk));
if (!cchunk) if (!cchunk)
error->all(FLERR,"Chunk/atom compute does not exist for compute chunk/spread/atom {}", idchunk); error->all(FLERR,"Chunk/atom compute {} does not exist for compute chunk/spread/atom "
"or is of invalid style", idchunk);
if (strcmp(cchunk->style,"chunk/atom") != 0) if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute chunk/spread/atom {} does not use chunk/atom compute", idchunk); error->all(FLERR,"Compute chunk/spread/atom {} does not use chunk/atom compute", idchunk);
} }
@ -196,14 +190,14 @@ void ComputeChunkSpreadAtom::compute_peratom()
// grow local vector_atom or array_atom if necessary // grow local vector_atom or array_atom if necessary
if (atom->nmax > nmax) { if (atom->nmax > nmax) {
if (nvalues == 1) { if (values.size() == 1) {
memory->destroy(vector_atom); memory->destroy(vector_atom);
nmax = atom->nmax; nmax = atom->nmax;
memory->create(vector_atom,nmax,"chunk/spread/atom:vector_atom"); memory->create(vector_atom,nmax,"chunk/spread/atom:vector_atom");
} else { } else {
memory->destroy(array_atom); memory->destroy(array_atom);
nmax = atom->nmax; nmax = atom->nmax;
memory->create(array_atom,nmax,nvalues,"chunk/spread/atom:array_atom"); memory->create(array_atom,nmax,values.size(),"chunk/spread/atom:array_atom");
} }
} }
@ -221,35 +215,35 @@ void ComputeChunkSpreadAtom::compute_peratom()
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int i,m,n,index,nstride; int index,nstride;
double *ptr; double *ptr;
for (m = 0; m < nvalues; m++) { int m = 0;
n = value2index[m]; for (auto &val : values) {
// copy compute/fix values into vector_atom or array_atom // copy compute/fix values into vector_atom or array_atom
// nstride between values for each atom // nstride between values for each atom
if (nvalues == 1) { if (values.size() == 1) {
ptr = vector_atom; ptr = vector_atom;
nstride = 1; nstride = 1;
} else { } else {
ptr = &array_atom[0][m]; ptr = &array_atom[0][m];
nstride = nvalues; nstride = values.size();
} }
// invoke compute if not previously invoked // invoke compute if not previously invoked
if (which[m] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n]; Compute *compute = val.val.c;
if (argindex[m] == 0) { if (val.argindex == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector(); compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR; compute->invoked_flag |= Compute::INVOKED_VECTOR;
} }
double *cvector = compute->vector; double *cvector = compute->vector;
for (i = 0; i < nlocal; i++, ptr += nstride) { for (int i = 0; i < nlocal; i++, ptr += nstride) {
*ptr = 0.0; *ptr = 0.0;
if (!(mask[i] & groupbit)) continue; if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1; index = ichunk[i]-1;
@ -262,9 +256,9 @@ void ComputeChunkSpreadAtom::compute_peratom()
compute->compute_array(); compute->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY; compute->invoked_flag |= Compute::INVOKED_ARRAY;
} }
int icol = argindex[m]-1; int icol = val.argindex-1;
double **carray = compute->array; double **carray = compute->array;
for (i = 0; i < nlocal; i++, ptr += nstride) { for (int i = 0; i < nlocal; i++, ptr += nstride) {
*ptr = 0.0; *ptr = 0.0;
if (!(mask[i] & groupbit)) continue; if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1; index = ichunk[i]-1;
@ -277,15 +271,15 @@ void ComputeChunkSpreadAtom::compute_peratom()
// are assuming the fix global vector/array is per-chunk data // are assuming the fix global vector/array is per-chunk data
// check if index exceeds fix output length/rows // check if index exceeds fix output length/rows
} else if (which[m] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
auto &fix = modify->get_fix_list()[n]; Fix *fix = val.val.f;
if (update->ntimestep % fix->global_freq) if (update->ntimestep % fix->global_freq)
error->all(FLERR,"Fix used in compute chunk/spread/atom not " error->all(FLERR,"Fix {} used in compute chunk/spread/atom not computed at compatible time",
"computed at compatible time"); val.id);
if (argindex[m] == 0) { if (val.argindex == 0) {
int nfix = fix->size_vector; int nfix = fix->size_vector;
for (i = 0; i < nlocal; i++, ptr += nstride) { for (int i = 0; i < nlocal; i++, ptr += nstride) {
*ptr = 0.0; *ptr = 0.0;
if (!(mask[i] & groupbit)) continue; if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1; index = ichunk[i]-1;
@ -294,9 +288,9 @@ void ComputeChunkSpreadAtom::compute_peratom()
} }
} else { } else {
int icol = argindex[m]-1; int icol = val.argindex-1;
int nfix = fix->size_array_rows; int nfix = fix->size_array_rows;
for (i = 0; i < nlocal; i++, ptr += nstride) { for (int i = 0; i < nlocal; i++, ptr += nstride) {
*ptr = 0.0; *ptr = 0.0;
if (!(mask[i] & groupbit)) continue; if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1; index = ichunk[i]-1;
@ -305,6 +299,7 @@ void ComputeChunkSpreadAtom::compute_peratom()
} }
} }
} }
++m;
} }
} }
@ -314,6 +309,7 @@ void ComputeChunkSpreadAtom::compute_peratom()
double ComputeChunkSpreadAtom::memory_usage() double ComputeChunkSpreadAtom::memory_usage()
{ {
double bytes = (double)nmax*nvalues * sizeof(double); double bytes = (double)nmax*values.size() * sizeof(double);
bytes += values.size() * sizeof(value_t);
return bytes; return bytes;
} }

17
src/compute_chunk_spread_atom.h
View File

@ -33,13 +33,20 @@ class ComputeChunkSpreadAtom : public Compute {
double memory_usage() override; double memory_usage() override;
protected: protected:
int mode, nvalues; struct value_t {
char *idchunk; int which;
char **ids; int argindex;
int *which, *argindex, *value2index; std::string id;
union {
class Compute *c;
class Fix *f;
} val;
};
std::vector<value_t> values;
int nmax; char *idchunk;
class ComputeChunkAtom *cchunk; class ComputeChunkAtom *cchunk;
int nmax;
void init_chunk(); void init_chunk();
}; };

633
src/fix_ave_time.cpp
View File

@ -19,6 +19,7 @@
#include "fix_ave_time.h" #include "fix_ave_time.h"
#include "arg_info.h" #include "arg_info.h"
#include "comm.h"
#include "compute.h" #include "compute.h"
#include "error.h" #include "error.h"
#include "input.h" #include "input.h"
@ -32,22 +33,18 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{ONE,RUNNING,WINDOW}; enum{ ONE, RUNNING, WINDOW };
enum{SCALAR,VECTOR}; enum{ SCALAR, VECTOR };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Fix(lmp, narg, arg),
nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr), nvalues(0), fp(nullptr), offlist(nullptr), format(nullptr), format_user(nullptr),
offcol(nullptr), varlen(nullptr), ids(nullptr),
fp(nullptr), offlist(nullptr), format(nullptr), format_user(nullptr),
vector(nullptr), vector_total(nullptr), vector_list(nullptr), vector(nullptr), vector_total(nullptr), vector_list(nullptr),
column(nullptr), array(nullptr), array_total(nullptr), array_list(nullptr) column(nullptr), array(nullptr), array_total(nullptr), array_list(nullptr)
{ {
if (narg < 7) error->all(FLERR,"Illegal fix ave/time command"); if (narg < 7) utils::missing_cmd_args(FLERR, "fix ave/time", error);
MPI_Comm_rank(world,&me);
nevery = utils::inumeric(FLERR,arg[3],false,lmp); nevery = utils::inumeric(FLERR,arg[3],false,lmp);
nrepeat = utils::inumeric(FLERR,arg[4],false,lmp); nrepeat = utils::inumeric(FLERR,arg[4],false,lmp);
@ -62,18 +59,17 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
// then read options so know mode = SCALAR/VECTOR before re-reading values // then read options so know mode = SCALAR/VECTOR before re-reading values
nvalues = 0; nvalues = 0;
int iarg = 6; int iarg = 6;
while (iarg < narg) { while (iarg < narg) {
if ((strncmp(arg[iarg],"c_",2) == 0) || if (utils::strmatch(arg[iarg],"^[cfv]_")) {
(strncmp(arg[iarg],"f_",2) == 0) ||
(strncmp(arg[iarg],"v_",2) == 0)) {
nvalues++; nvalues++;
iarg++; iarg++;
} else break; } else break;
} }
if (nvalues == 0)
error->all(FLERR,"No values from computes, fixes, or variables used in fix ave/time command");
if (nvalues == 0) error->all(FLERR,"No values in fix ave/time command"); // parse optional keywords
options(iarg,narg,arg); options(iarg,narg,arg);
@ -83,7 +79,6 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
int expand = 0; int expand = 0;
char **earg; char **earg;
nvalues = utils::expand_args(FLERR,nvalues,&arg[6],mode,earg,lmp); nvalues = utils::expand_args(FLERR,nvalues,&arg[6],mode,earg,lmp);
keyword.resize(nvalues);
key2col.clear(); key2col.clear();
if (earg != &arg[6]) expand = 1; if (earg != &arg[6]) expand = 1;
@ -91,132 +86,126 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
// parse values // parse values
which = new int[nvalues]; values.clear();
argindex = new int[nvalues];
value2index = new int[nvalues];
offcol = new int[nvalues];
varlen = new int[nvalues];
ids = new char*[nvalues];
for (int i = 0; i < nvalues; i++) { for (int i = 0; i < nvalues; i++) {
ArgInfo argi(arg[i]); ArgInfo argi(arg[i]);
keyword[i] = arg[i];
value_t val;
val.keyword = arg[i];
key2col[arg[i]] = i; key2col[arg[i]] = i;
if ((argi.get_type() == ArgInfo::NONE) if ((argi.get_type() == ArgInfo::NONE)
|| (argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_type() == ArgInfo::UNKNOWN)
|| (argi.get_dim() > 1)) || (argi.get_dim() > 1))
error->all(FLERR,"Invalid fix ave/time command"); error->all(FLERR,"Invalid fix ave/time argument: {}", arg[i]);
which[i] = argi.get_type(); val.which = argi.get_type();
argindex[i] = argi.get_index1(); val.argindex = argi.get_index1();
ids[i] = argi.copy_name(); val.varlen = 0;
val.offcol = 0;
val.id = argi.get_name();
val.val.c = nullptr;
values.push_back(val);
} }
if (nvalues != (int)values.size())
error->all(FLERR, "Could not parse value data consistently for fix ave/time");
// set off columns now that nvalues is finalized // set off columns now that nvalues is finalized
for (int i = 0; i < nvalues; i++) offcol[i] = 0;
for (int i = 0; i < noff; i++) { for (int i = 0; i < noff; i++) {
if (offlist[i] < 1 || offlist[i] > nvalues) if (offlist[i] < 1 || offlist[i] > nvalues)
error->all(FLERR,"Invalid fix ave/time off column"); error->all(FLERR,"Invalid fix ave/time off column: {}", offlist[i]);
offcol[offlist[i]-1] = 1; values[offlist[i]-1].offcol = 1;
} }
// setup and error check // setup and error check
// for fix inputs, check that fix frequency is acceptable // for fix inputs, check that fix frequency is acceptable
// set variable_length if any compute is variable length // set variable_length if any compute is variable length
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0) if (nevery <= 0)
error->all(FLERR,"Illegal fix ave/time command"); error->all(FLERR,"Illegal fix ave/time nevery value: {}", nevery);
if (nrepeat <= 0)
error->all(FLERR,"Illegal fix ave/time nrepeat value: {}", nrepeat);
if (nfreq <= 0)
error->all(FLERR,"Illegal fix ave/time nfreq value: {}", nfreq);
if (nfreq % nevery || nrepeat*nevery > nfreq) if (nfreq % nevery || nrepeat*nevery > nfreq)
error->all(FLERR,"Illegal fix ave/time command"); error->all(FLERR,"Inconsistent fix ave/time nevery/nrepeat/nfreq values");
if (ave != RUNNING && overwrite) if (ave != RUNNING && overwrite)
error->all(FLERR,"Illegal fix ave/time command"); error->all(FLERR,"Fix ave/time overwrite keyword requires ave running setting");
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
varlen[i] = 0;
if (which[i] == ArgInfo::COMPUTE && mode == SCALAR) { if ((val.which == ArgInfo::COMPUTE) && (mode == SCALAR)) {
int icompute = modify->find_compute(ids[i]); val.val.c = modify->get_compute_by_id(val.id);
if (icompute < 0) if (!val.val.c) error->all(FLERR,"Compute ID {} for fix ave/time does not exist", val.id);
error->all(FLERR,"Compute ID for fix ave/time does not exist"); if (val.argindex == 0 && (val.val.c->scalar_flag == 0))
if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0) error->all(FLERR,"Fix ave/time compute {} does not calculate a scalar", val.id);
error->all(FLERR,"Fix ave/time compute does not calculate a scalar"); if (val.argindex && (val.val.c->vector_flag == 0))
if (argindex[i] && modify->compute[icompute]->vector_flag == 0) error->all(FLERR,"Fix ave/time compute {} does not calculate a vector", val.id);
error->all(FLERR,"Fix ave/time compute does not calculate a vector"); if (val.argindex && (val.argindex > val.val.c->size_vector) &&
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector && (val.val.c->size_vector_variable == 0))
modify->compute[icompute]->size_vector_variable == 0) error->all(FLERR, "Fix ave/time compute {} vector is accessed out-of-range", val.id);
error->all(FLERR, if (val.argindex && val.val.c->size_vector_variable) val.varlen = 1;
"Fix ave/time compute vector is accessed out-of-range");
if (argindex[i] && modify->compute[icompute]->size_vector_variable)
varlen[i] = 1;
} else if (which[i] == ArgInfo::COMPUTE && mode == VECTOR) { } else if ((val.which == ArgInfo::COMPUTE) && (mode == VECTOR)) {
int icompute = modify->find_compute(ids[i]); val.val.c = modify->get_compute_by_id(val.id);
if (icompute < 0) if (!val.val.c) error->all(FLERR,"Compute ID {} for fix ave/time does not exist", val.id);
error->all(FLERR,"Compute ID for fix ave/time does not exist"); if ((val.argindex == 0) && (val.val.c->vector_flag == 0))
if (argindex[i] == 0 && modify->compute[icompute]->vector_flag == 0) error->all(FLERR,"Fix ave/time compute {} does not calculate a vector", val.id);
error->all(FLERR,"Fix ave/time compute does not calculate a vector"); if (val.argindex && (val.val.c->array_flag == 0))
if (argindex[i] && modify->compute[icompute]->array_flag == 0) error->all(FLERR,"Fix ave/time compute {} does not calculate an array", val.id);
error->all(FLERR,"Fix ave/time compute does not calculate an array"); if (val.argindex && (val.argindex > val.val.c->size_array_cols))
if (argindex[i] && error->all(FLERR,"Fix ave/time compute {} array is accessed out-of-range", val.id);
argindex[i] > modify->compute[icompute]->size_array_cols) if ((val.argindex == 0) && (val.val.c->size_vector_variable)) val.varlen = 1;
error->all(FLERR,"Fix ave/time compute array is accessed out-of-range"); if (val.argindex && (val.val.c->size_array_rows_variable)) val.varlen = 1;
if (argindex[i] == 0 && modify->compute[icompute]->size_vector_variable)
varlen[i] = 1;
if (argindex[i] && modify->compute[icompute]->size_array_rows_variable)
varlen[i] = 1;
} else if (which[i] == ArgInfo::FIX && mode == SCALAR) { } else if ((val.which == ArgInfo::FIX) && (mode == SCALAR)) {
int ifix = modify->find_fix(ids[i]); val.val.f = modify->get_fix_by_id(val.id);
if (ifix < 0) if (!val.val.f) error->all(FLERR,"Fix ID {} for fix ave/time does not exist", val.id);
error->all(FLERR,"Fix ID for fix ave/time does not exist"); if ((val.argindex == 0) && (val.val.f->scalar_flag == 0))
if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0) error->all(FLERR,"Fix ave/time fix {} does not calculate a scalar", val.id);
error->all(FLERR,"Fix ave/time fix does not calculate a scalar"); if (val.argindex && (val.val.f->vector_flag == 0))
if (argindex[i] && modify->fix[ifix]->vector_flag == 0) error->all(FLERR,"Fix ave/time fix {} does not calculate a vector", val.id);
error->all(FLERR,"Fix ave/time fix does not calculate a vector"); if (val.argindex && (val.val.f->size_vector_variable))
if (argindex[i] && modify->fix[ifix]->size_vector_variable) error->all(FLERR,"Fix ave/time fix {} vector cannot be variable length", val.id);
error->all(FLERR,"Fix ave/time fix vector cannot be variable length"); if (val.argindex && (val.argindex > val.val.f->size_vector))
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector) error->all(FLERR,"Fix ave/time fix {} vector is accessed out-of-range", val.id);
error->all(FLERR,"Fix ave/time fix vector is accessed out-of-range"); if (nevery % val.val.f->global_freq)
if (nevery % modify->fix[ifix]->global_freq) error->all(FLERR, "Fix {} for fix ave/time not computed at compatible time", val.id);
error->all(FLERR,
"Fix for fix ave/time not computed at compatible time");
} else if (which[i] == ArgInfo::FIX && mode == VECTOR) { } else if ((val.which == ArgInfo::FIX) && (mode == VECTOR)) {
int ifix = modify->find_fix(ids[i]); val.val.f = modify->get_fix_by_id(val.id);
if (ifix < 0) if (!val.val.f) error->all(FLERR,"Fix ID {} for fix ave/time does not exist", val.id);
error->all(FLERR,"Fix ID for fix ave/time does not exist"); if ((val.argindex == 0) && (val.val.f->vector_flag == 0))
if (argindex[i] == 0 && modify->fix[ifix]->vector_flag == 0) error->all(FLERR,"Fix ave/time fix {} does not calculate a vector", val.id);
error->all(FLERR,"Fix ave/time fix does not calculate a vector"); if (val.argindex && (val.val.f->array_flag == 0))
if (argindex[i] && modify->fix[ifix]->array_flag == 0) error->all(FLERR,"Fix ave/time fix {} does not calculate an array", val.id);
error->all(FLERR,"Fix ave/time fix does not calculate an array"); if (val.argindex && (val.val.f->size_array_rows_variable))
if (argindex[i] && modify->fix[ifix]->size_array_rows_variable) error->all(FLERR,"Fix ave/time fix {} array cannot be variable length", val.id);
error->all(FLERR,"Fix ave/time fix array cannot be variable length"); if (val.argindex && (val.argindex > val.val.f->size_array_cols))
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_array_cols) error->all(FLERR,"Fix ave/time fix {} array is accessed out-of-range", val.id);
error->all(FLERR,"Fix ave/time fix array is accessed out-of-range"); if (nevery % val.val.f->global_freq)
if (nevery % modify->fix[ifix]->global_freq) error->all(FLERR, "Fix {} for fix ave/time not computed at compatible time", val.id);
error->all(FLERR,
"Fix for fix ave/time not computed at compatible time");
} else if (which[i] == ArgInfo::VARIABLE && mode == SCALAR) { } else if ((val.which == ArgInfo::VARIABLE) && (mode == SCALAR)) {
int ivariable = input->variable->find(ids[i]); int ivariable = input->variable->find(val.id.c_str());
if (ivariable < 0) if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/time does not exist"); error->all(FLERR,"Variable name {} for fix ave/time does not exist", val.id);
if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0) if ((val.argindex == 0) && (input->variable->equalstyle(ivariable) == 0))
error->all(FLERR,"Fix ave/time variable is not equal-style variable"); error->all(FLERR,"Fix ave/time variable {} is not equal-style variable", val.id);
if (argindex[i] && input->variable->vectorstyle(ivariable) == 0) if ((val.argindex) && (input->variable->vectorstyle(ivariable) == 0))
error->all(FLERR,"Fix ave/time variable is not vector-style variable"); error->all(FLERR,"Fix ave/time variable {} is not vector-style variable", val.id);
} else if (which[i] == ArgInfo::VARIABLE && mode == VECTOR) { } else if ((val.which == ArgInfo::VARIABLE) && (mode == VECTOR)) {
int ivariable = input->variable->find(ids[i]); int ivariable = input->variable->find(val.id.c_str());
if (ivariable < 0) if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/time does not exist"); error->all(FLERR,"Variable name {} for fix ave/time does not exist", val.id);
if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0) if ((val.argindex == 0) && (input->variable->vectorstyle(ivariable) == 0))
error->all(FLERR,"Fix ave/time variable is not vector-style variable"); error->all(FLERR,"Fix ave/time variable {} is not vector-style variable", val.id);
if (argindex[i]) if (val.argindex)
error->all(FLERR,"Fix ave/time mode vector variable cannot be indexed"); error->all(FLERR,"Fix ave/time mode vector variable {} cannot be indexed", val.id);
varlen[i] = 1; val.varlen = 1;
} }
} }
@ -225,9 +214,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
all_variable_length = 1; all_variable_length = 1;
any_variable_length = 0; any_variable_length = 0;
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
if (varlen[i] == 0) all_variable_length = 0; if (val.varlen == 0) all_variable_length = 0;
if (varlen[i]) any_variable_length = 1; if (val.varlen) any_variable_length = 1;
} }
// if VECTOR mode, check that all columns are same length // if VECTOR mode, check that all columns are same length
@ -245,19 +234,18 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
// only if nrepeat > 1 or ave = RUNNING/WINDOW, // only if nrepeat > 1 or ave = RUNNING/WINDOW,
// so that locking spans multiple timesteps // so that locking spans multiple timesteps
if (any_variable_length && if (any_variable_length && ((nrepeat > 1) || (ave == RUNNING) || (ave == WINDOW))) {
(nrepeat > 1 || ave == RUNNING || ave == WINDOW)) { for (auto &val : values) {
for (int i = 0; i < nvalues; i++) if (val.varlen && val.which == ArgInfo::COMPUTE) val.val.c->lock_enable();
if (varlen[i] && which[i] == ArgInfo::COMPUTE) lockforever = 0;
modify->get_compute_by_id(ids[i])->lock_enable(); }
lockforever = 0;
} }
// print file comment lines // print file comment lines
// for mode = VECTOR, cannot use arg to print // for mode = VECTOR, cannot use arg to print
// since array args may have been expanded to multiple vectors // since array args may have been expanded to multiple vectors
if (fp && me == 0) { if (fp && comm->me == 0) {
clearerr(fp); clearerr(fp);
if (title1) fprintf(fp,"%s\n",title1); if (title1) fprintf(fp,"%s\n",title1);
else fprintf(fp,"# Time-averaged data for fix %s\n",id); else fprintf(fp,"# Time-averaged data for fix %s\n",id);
@ -300,8 +288,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
if (mode == SCALAR) { if (mode == SCALAR) {
vector = new double[nvalues]; vector = new double[nvalues];
vector_total = new double[nvalues]; vector_total = new double[nvalues];
if (ave == WINDOW) if (ave == WINDOW) memory->create(vector_list,nwindow,nvalues,"ave/time:vector_list");
memory->create(vector_list,nwindow,nvalues,"ave/time:vector_list");
} else allocate_arrays(); } else allocate_arrays();
// this fix produces either a global scalar or vector or array // this fix produces either a global scalar or vector or array
@ -314,17 +301,16 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
if (mode == SCALAR) { if (mode == SCALAR) {
if (nvalues == 1) { if (nvalues == 1) {
scalar_flag = 1; scalar_flag = 1;
if (which[0] == ArgInfo::COMPUTE) { auto &val = values[0];
Compute *compute = modify->compute[modify->find_compute(ids[0])]; if (val.which == ArgInfo::COMPUTE) {
if (argindex[0] == 0) extscalar = compute->extscalar; if (val.argindex == 0) extscalar = val.val.c->extscalar;
else if (compute->extvector >= 0) extscalar = compute->extvector; else if (val.val.c->extvector >= 0) extscalar = val.val.c->extvector;
else extscalar = compute->extlist[argindex[0]-1]; else extscalar = val.val.c->extlist[val.argindex-1];
} else if (which[0] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
Fix *fix = modify->fix[modify->find_fix(ids[0])]; if (val.argindex == 0) extscalar = val.val.f->extscalar;
if (argindex[0] == 0) extscalar = fix->extscalar; else if (val.val.f->extvector >= 0) extscalar = val.val.f->extvector;
else if (fix->extvector >= 0) extscalar = fix->extvector; else extscalar = val.val.f->extlist[val.argindex-1];
else extscalar = fix->extlist[argindex[0]-1]; } else if (val.which == ArgInfo::VARIABLE) {
} else if (which[0] == ArgInfo::VARIABLE) {
extscalar = 0; extscalar = 0;
} }
@ -333,47 +319,46 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
size_vector = nrows = nvalues; size_vector = nrows = nvalues;
extvector = -1; extvector = -1;
extlist = new int[nvalues]; extlist = new int[nvalues];
for (int i = 0; i < nvalues; i++) { int i = 0;
if (which[i] == ArgInfo::COMPUTE) { for (auto &val : values) {
Compute *compute = modify->compute[modify->find_compute(ids[i])]; if (val.which == ArgInfo::COMPUTE) {
if (argindex[i] == 0) extlist[i] = compute->extscalar; if (val.argindex == 0) extlist[i] = val.val.c->extscalar;
else if (compute->extvector >= 0) extlist[i] = compute->extvector; else if (val.val.c->extvector >= 0) extlist[i] = val.val.c->extvector;
else extlist[i] = compute->extlist[argindex[i]-1]; else extlist[i] = val.val.c->extlist[val.argindex-1];
} else if (which[i] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
Fix *fix = modify->fix[modify->find_fix(ids[i])]; if (val.argindex == 0) extlist[i] = val.val.f->extscalar;
if (argindex[i] == 0) extlist[i] = fix->extscalar; else if (val.val.f->extvector >= 0) extlist[i] = val.val.f->extvector;
else if (fix->extvector >= 0) extlist[i] = fix->extvector; else extlist[i] = val.val.f->extlist[val.argindex-1];
else extlist[i] = fix->extlist[argindex[i]-1]; } else if (val.which == ArgInfo::VARIABLE) {
} else if (which[i] == ArgInfo::VARIABLE) {
extlist[i] = 0; extlist[i] = 0;
} }
++i;
} }
} }
} else { } else {
if (nvalues == 1) { if (nvalues == 1) {
auto &val = values[0];
vector_flag = 1; vector_flag = 1;
size_vector = nrows; size_vector = nrows;
if (all_variable_length) size_vector_variable = 1; if (all_variable_length) size_vector_variable = 1;
if (which[0] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[modify->find_compute(ids[0])]; if (val.argindex == 0) {
if (argindex[0] == 0) { extvector = val.val.c->extvector;
extvector = compute->extvector;
if (extvector == -1) { if (extvector == -1) {
extlist = new int[nrows]; extlist = new int[nrows];
for (int i = 0; i < nrows; i++) extlist[i] = compute->extlist[i]; for (int i = 0; i < nrows; i++) extlist[i] = val.val.c->extlist[i];
} }
} else extvector = compute->extarray; } else extvector = val.val.c->extarray;
} else if (which[0] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
Fix *fix = modify->fix[modify->find_fix(ids[0])]; if (val.argindex == 0) {
if (argindex[0] == 0) { extvector = val.val.f->extvector;
extvector = fix->extvector;
if (extvector == -1) { if (extvector == -1) {
extlist = new int[nrows]; extlist = new int[nrows];
for (int i = 0; i < nrows; i++) extlist[i] = fix->extlist[i]; for (int i = 0; i < nrows; i++) extlist[i] = val.val.f->extlist[i];
} }
} else extvector = fix->extarray; } else extvector = val.val.f->extarray;
} else if (which[0] == ArgInfo::VARIABLE) { } else if (val.which == ArgInfo::VARIABLE) {
extlist = new int[nrows]; extlist = new int[nrows];
for (int i = 0; i < nrows; i++) extlist[i] = 0; for (int i = 0; i < nrows; i++) extlist[i] = 0;
} }
@ -383,26 +368,25 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
size_array_rows = nrows; size_array_rows = nrows;
size_array_cols = nvalues; size_array_cols = nvalues;
if (all_variable_length) size_array_rows_variable = 1; if (all_variable_length) size_array_rows_variable = 1;
int value; int extvalue = 0;
for (int i = 0; i < nvalues; i++) { extarray = -2;
if (which[i] == ArgInfo::COMPUTE) { for (auto &val : values) {
Compute *compute = modify->compute[modify->find_compute(ids[i])]; if (val.which == ArgInfo::COMPUTE) {
if (argindex[i] == 0) value = compute->extvector; if (val.argindex == 0) extvalue = val.val.c->extvector;
else value = compute->extarray; else extvalue = val.val.c->extarray;
} else if (which[i] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
Fix *fix = modify->fix[modify->find_fix(ids[i])]; if (val.argindex == 0) extvalue = val.val.f->extvector;
if (argindex[i] == 0) value = fix->extvector; else extvalue = val.val.f->extarray;
else value = fix->extarray; } else if (val.which == ArgInfo::VARIABLE) {
} else if (which[i] == ArgInfo::VARIABLE) { extvalue = 0;
value = 0;
} }
if (value == -1) if (extvalue == -1)
error->all(FLERR,"Fix ave/time cannot set output array " error->all(FLERR,"Fix ave/time cannot set output array intensive/extensive "
"intensive/extensive from these inputs"); "from these inputs");
if (i == 0) extarray = value; if (extarray < -1) extarray = extvalue;
else if (value != extarray) else if (extvalue != extarray)
error->all(FLERR,"Fix ave/time cannot set output array " error->all(FLERR,"Fix ave/time cannot set output array intensive/extensive "
"intensive/extensive from these inputs"); "from these inputs");
} }
} }
} }
@ -434,32 +418,23 @@ FixAveTime::~FixAveTime()
{ {
// decrement lock counter in compute chunk/atom, it if still exists // decrement lock counter in compute chunk/atom, it if still exists
if (any_variable_length && if (any_variable_length && ((nrepeat > 1) || (ave == RUNNING) || (ave == WINDOW))) {
(nrepeat > 1 || ave == RUNNING || ave == WINDOW)) { for (auto &val : values) {
for (int i = 0; i < nvalues; i++) if (val.varlen) {
if (varlen[i]) { auto icompute = modify->get_compute_by_id(val.id);
int icompute = modify->find_compute(ids[i]); if (icompute) {
if (icompute >= 0) { if ((ave == RUNNING) || (ave == WINDOW))
if (ave == RUNNING || ave == WINDOW) icompute->unlock(this);
modify->compute[icompute]->unlock(this); icompute->lock_disable();
modify->compute[icompute]->lock_disable();
} }
} }
}
} }
delete[] format_user; delete[] format_user;
delete[] which;
delete[] argindex;
delete[] value2index;
delete[] offcol;
delete[] varlen;
for (int i = 0; i < nvalues; i++) delete[] ids[i];
delete[] ids;
delete[] extlist; delete[] extlist;
if (fp && me == 0) { if (fp && comm->me == 0) {
if (yaml_flag) fputs("...\n", fp); if (yaml_flag) fputs("...\n", fp);
fclose(fp); fclose(fp);
} }
@ -485,24 +460,21 @@ int FixAveTime::setmask()
void FixAveTime::init() void FixAveTime::init()
{ {
// set current indices for all computes,fixes,variables // update indices/pointers for all computes,fixes,variables
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
if (which[i] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[i]); val.val.c = modify->get_compute_by_id(val.id);
if (icompute < 0) if (!val.val.c)
error->all(FLERR,"Compute ID for fix ave/time does not exist"); error->all(FLERR,"Compute ID {} for fix ave/time does not exist", val.id);
value2index[i] = icompute; } else if (val.which == ArgInfo::FIX) {
} else if (which[i] == ArgInfo::FIX) { val.val.f = modify->get_fix_by_id(val.id);
int ifix = modify->find_fix(ids[i]); if (!val.val.f)
if (ifix < 0) error->all(FLERR,"Fix ID {} for fix ave/time does not exist", val.id);
error->all(FLERR,"Fix ID for fix ave/time does not exist"); } else if (val.which == ArgInfo::VARIABLE) {
value2index[i] = ifix; val.val.v = input->variable->find(val.id.c_str());
} else if (which[i] == ArgInfo::VARIABLE) { if (val.val.v < 0)
int ivariable = input->variable->find(ids[i]); error->all(FLERR,"Variable name {} for fix ave/time does not exist", val.id);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/time does not exist");
value2index[i] = ivariable;
} }
} }
@ -542,9 +514,6 @@ void FixAveTime::end_of_step()
void FixAveTime::invoke_scalar(bigint ntimestep) void FixAveTime::invoke_scalar(bigint ntimestep)
{ {
int i,m;
double scalar;
// zero if first sample within single Nfreq epoch // zero if first sample within single Nfreq epoch
// if any input is variable length, initialize current length // if any input is variable length, initialize current length
// check for exceeding length is done below // check for exceeding length is done below
@ -556,7 +525,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
modify->addstep_compute(ntimestep+nevery); modify->addstep_compute(ntimestep+nevery);
modify->addstep_compute(ntimestep+nfreq); modify->addstep_compute(ntimestep+nfreq);
} }
for (i = 0; i < nvalues; i++) vector[i] = 0.0; for (int i = 0; i < nvalues; i++) vector[i] = 0.0;
} }
// accumulate results of computes,fixes,variables to local copy // accumulate results of computes,fixes,variables to local copy
@ -564,56 +533,57 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
modify->clearstep_compute(); modify->clearstep_compute();
for (i = 0; i < nvalues; i++) { int i = 0;
m = value2index[i]; double scalar = 0.0;
for (auto &val : values) {
// invoke compute if not previously invoked // invoke compute if not previously invoked
// insure no out-of-range access to variable-length compute vector // insure no out-of-range access to variable-length compute vector
if (which[i] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[m];
if (argindex[i] == 0) { if (val.argindex == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) {
compute->compute_scalar(); val.val.c->compute_scalar();
compute->invoked_flag |= Compute::INVOKED_SCALAR; val.val.c->invoked_flag |= Compute::INVOKED_SCALAR;
} }
scalar = compute->scalar; scalar = val.val.c->scalar;
} else { } else {
if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector(); val.val.c->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR; val.val.c->invoked_flag |= Compute::INVOKED_VECTOR;
} }
if (varlen[i] && compute->size_vector < argindex[i]) scalar = 0.0; if (val.varlen && (val.val.c->size_vector < val.argindex)) scalar = 0.0;
else scalar = compute->vector[argindex[i]-1]; else scalar = val.val.c->vector[val.argindex-1];
} }
// access fix fields, guaranteed to be ready // access fix fields, guaranteed to be ready
} else if (which[i] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
if (argindex[i] == 0) if (val.argindex == 0)
scalar = modify->fix[m]->compute_scalar(); scalar = val.val.f->compute_scalar();
else else
scalar = modify->fix[m]->compute_vector(argindex[i]-1); scalar = val.val.f->compute_vector(val.argindex-1);
// evaluate equal-style or vector-style variable // evaluate equal-style or vector-style variable
// insure no out-of-range access to vector-style variable // insure no out-of-range access to vector-style variable
} else if (which[i] == ArgInfo::VARIABLE) { } else if (val.which == ArgInfo::VARIABLE) {
if (argindex[i] == 0) if (val.argindex == 0)
scalar = input->variable->compute_equal(m); scalar = input->variable->compute_equal(val.val.v);
else { else {
double *varvec; double *varvec;
int nvec = input->variable->compute_vector(m,&varvec); int nvec = input->variable->compute_vector(val.val.v,&varvec);
if (nvec < argindex[i]) scalar = 0.0; if (nvec < val.argindex) scalar = 0.0;
else scalar = varvec[argindex[i]-1]; else scalar = varvec[val.argindex-1];
} }
} }
// add value to vector or just set directly if offcol is set // add value to vector or just set directly if offcol is set
if (offcol[i]) vector[i] = scalar; if (val.offcol) vector[i] = scalar;
else vector[i] += scalar; else vector[i] += scalar;
++i;
} }
// done if irepeat < nrepeat // done if irepeat < nrepeat
@ -634,7 +604,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
double repeat = nrepeat; double repeat = nrepeat;
for (i = 0; i < nvalues; i++) for (i = 0; i < nvalues; i++)
if (offcol[i] == 0) vector[i] /= repeat; if (values[i].offcol == 0) vector[i] /= repeat;
// if ave = ONE, only single Nfreq timestep value is needed // if ave = ONE, only single Nfreq timestep value is needed
// if ave = RUNNING, combine with all previous Nfreq timestep values // if ave = RUNNING, combine with all previous Nfreq timestep values
@ -667,18 +637,18 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
// insure any columns with offcol set are effectively set to last value // insure any columns with offcol set are effectively set to last value
for (i = 0; i < nvalues; i++) for (i = 0; i < nvalues; i++)
if (offcol[i]) vector_total[i] = norm*vector[i]; if (values[i].offcol) vector_total[i] = norm*vector[i];
// output result to file // output result to file
if (fp && me == 0) { if (fp && comm->me == 0) {
clearerr(fp); clearerr(fp);
if (overwrite) platform::fseek(fp,filepos); if (overwrite) platform::fseek(fp,filepos);
if (yaml_flag) { if (yaml_flag) {
if (!yaml_header || overwrite) { if (!yaml_header || overwrite) {
yaml_header = true; yaml_header = true;
fputs("keywords: ['Step', ", fp); fputs("keywords: ['Step', ", fp);
for (const auto &k : keyword) fmt::print(fp, "'{}', ", k); for (const auto &val : values) fmt::print(fp, "'{}', ", val.keyword);
fputs("]\ndata:\n", fp); fputs("]\ndata:\n", fp);
} }
fmt::print(fp, " - [{}, ", ntimestep); fmt::print(fp, " - [{}, ", ntimestep);
@ -704,8 +674,6 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
void FixAveTime::invoke_vector(bigint ntimestep) void FixAveTime::invoke_vector(bigint ntimestep)
{ {
int i,j,m;
// first sample within single Nfreq epoch // first sample within single Nfreq epoch
// zero out arrays that accumulate over many samples, but not across epochs // zero out arrays that accumulate over many samples, but not across epochs
// invoke setup_chunks() to determine current nchunk // invoke setup_chunks() to determine current nchunk
@ -735,22 +703,20 @@ void FixAveTime::invoke_vector(bigint ntimestep)
} }
int lockforever_flag = 0; int lockforever_flag = 0;
for (i = 0; i < nvalues; i++) { for (auto &val : values) {
if (!varlen[i] || which[i] != ArgInfo::COMPUTE) continue; if (!val.varlen || (val.which != ArgInfo::COMPUTE)) continue;
if (nrepeat > 1 && ave == ONE) { if ((nrepeat > 1) && (ave == ONE)) {
Compute *compute = modify->compute[value2index[i]]; val.val.c->lock(this,ntimestep,ntimestep+static_cast<bigint>(nrepeat-1)*nevery);
compute->lock(this,ntimestep,ntimestep+static_cast<bigint>(nrepeat-1)*nevery); } else if (((ave == RUNNING) || (ave == WINDOW)) && !lockforever) {
} else if ((ave == RUNNING || ave == WINDOW) && !lockforever) { val.val.c->lock(this,update->ntimestep,-1);
Compute *compute = modify->compute[value2index[i]];
compute->lock(this,update->ntimestep,-1);
lockforever_flag = 1; lockforever_flag = 1;
} }
} }
if (lockforever_flag) lockforever = 1; if (lockforever_flag) lockforever = 1;
} }
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) array[i][j] = 0.0; for (int j = 0; j < nvalues; j++) array[i][j] = 0.0;
} }
// accumulate results of computes,fixes,variables to local copy // accumulate results of computes,fixes,variables to local copy
@ -758,69 +724,67 @@ void FixAveTime::invoke_vector(bigint ntimestep)
modify->clearstep_compute(); modify->clearstep_compute();
for (j = 0; j < nvalues; j++) { int j = 0;
m = value2index[j]; for (auto &val : values) {
// invoke compute if not previously invoked // invoke compute if not previously invoked
if (which[j] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[m]; if (val.argindex == 0) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) {
if (argindex[j] == 0) { val.val.c->compute_vector();
if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { val.val.c->invoked_flag |= Compute::INVOKED_VECTOR;
compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR;
} }
double *cvector = compute->vector; double *cvector = val.val.c->vector;
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
column[i] = cvector[i]; column[i] = cvector[i];
} else { } else {
if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { if (!(val.val.c->invoked_flag & Compute::INVOKED_ARRAY)) {
compute->compute_array(); val.val.c->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY; val.val.c->invoked_flag |= Compute::INVOKED_ARRAY;
} }
double **carray = compute->array; double **carray = val.val.c->array;
int icol = argindex[j]-1; int icol = val.argindex-1;
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
column[i] = carray[i][icol]; column[i] = carray[i][icol];
} }
// access fix fields, guaranteed to be ready // access fix fields, guaranteed to be ready
} else if (which[j] == ArgInfo::FIX) { } else if (val.which == ArgInfo::FIX) {
Fix *fix = modify->fix[m]; if (val.argindex == 0)
if (argindex[j] == 0) for (int i = 0; i < nrows; i++)
for (i = 0; i < nrows; i++) column[i] = val.val.f->compute_vector(i);
column[i] = fix->compute_vector(i);
else { else {
int icol = argindex[j]-1; int icol = val.argindex-1;
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
column[i] = fix->compute_array(i,icol); column[i] = val.val.f->compute_array(i,icol);
} }
// evaluate vector-style variable // evaluate vector-style variable
// insure nvec = nrows, else error // insure nvec = nrows, else error
// could be different on this timestep than when column_length(1) set nrows // could be different on this timestep than when column_length(1) set nrows
} else if (which[j] == ArgInfo::VARIABLE) { } else if (val.which == ArgInfo::VARIABLE) {
double *varvec; double *varvec;
int nvec = input->variable->compute_vector(m,&varvec); int nvec = input->variable->compute_vector(val.val.v,&varvec);
if (nvec != nrows) if (nvec != nrows)
error->all(FLERR,"Fix ave/time vector-style variable changed length"); error->all(FLERR,"Fix ave/time vector-style variable {} changed length", val.id);
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
column[i] = varvec[i]; column[i] = varvec[i];
} }
// add columns of values to array or just set directly if offcol is set // add columns of values to array or just set directly if offcol is set
if (offcol[j]) { if (val.offcol) {
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
array[i][j] = column[i]; array[i][j] = column[i];
} else { } else {
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
array[i][j] += column[i]; array[i][j] += column[i];
} }
++j;
} }
// done if irepeat < nrepeat // done if irepeat < nrepeat
@ -840,38 +804,37 @@ void FixAveTime::invoke_vector(bigint ntimestep)
// unlock any variable length computes at end of Nfreq epoch // unlock any variable length computes at end of Nfreq epoch
// do not unlock if ave = RUNNING or WINDOW // do not unlock if ave = RUNNING or WINDOW
if (any_variable_length && nrepeat > 1 && ave == ONE) { if (any_variable_length && (nrepeat > 1) && (ave == ONE)) {
for (i = 0; i < nvalues; i++) { for (auto &val : values) {
if (!varlen[i]) continue; if (!val.varlen) continue;
Compute *compute = modify->compute[value2index[i]]; if ((val.which == ArgInfo::COMPUTE) && val.val.c) val.val.c->unlock(this);
compute->unlock(this);
} }
} }
// average the final result for the Nfreq timestep // average the final result for the Nfreq timestep
double repeat = nrepeat; double repeat = nrepeat;
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) for (int j = 0; j < nvalues; j++)
if (offcol[j] == 0) array[i][j] /= repeat; if (values[j].offcol == 0) array[i][j] /= repeat;
// if ave = ONE, only single Nfreq timestep value is needed // if ave = ONE, only single Nfreq timestep value is needed
// if ave = RUNNING, combine with all previous Nfreq timestep values // if ave = RUNNING, combine with all previous Nfreq timestep values
// if ave = WINDOW, combine with nwindow most recent Nfreq timestep values // if ave = WINDOW, combine with nwindow most recent Nfreq timestep values
if (ave == ONE) { if (ave == ONE) {
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) array_total[i][j] = array[i][j]; for (int j = 0; j < nvalues; j++) array_total[i][j] = array[i][j];
norm = 1; norm = 1;
} else if (ave == RUNNING) { } else if (ave == RUNNING) {
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) array_total[i][j] += array[i][j]; for (int j = 0; j < nvalues; j++) array_total[i][j] += array[i][j];
norm++; norm++;
} else if (ave == WINDOW) { } else if (ave == WINDOW) {
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) { for (int j = 0; j < nvalues; j++) {
array_total[i][j] += array[i][j]; array_total[i][j] += array[i][j];
if (window_limit) array_total[i][j] -= array_list[iwindow][i][j]; if (window_limit) array_total[i][j] -= array_list[iwindow][i][j];
array_list[iwindow][i][j] = array[i][j]; array_list[iwindow][i][j] = array[i][j];
@ -888,32 +851,32 @@ void FixAveTime::invoke_vector(bigint ntimestep)
// insure any columns with offcol set are effectively set to last value // insure any columns with offcol set are effectively set to last value
for (i = 0; i < nrows; i++) for (int i = 0; i < nrows; i++)
for (j = 0; j < nvalues; j++) for (int j = 0; j < nvalues; j++)
if (offcol[j]) array_total[i][j] = norm*array[i][j]; if (values[j].offcol) array_total[i][j] = norm*array[i][j];
// output result to file // output result to file
if (fp && me == 0) { if (fp && comm->me == 0) {
if (overwrite) platform::fseek(fp,filepos); if (overwrite) platform::fseek(fp,filepos);
if (yaml_flag) { if (yaml_flag) {
if (!yaml_header || overwrite) { if (!yaml_header || overwrite) {
yaml_header = true; yaml_header = true;
fputs("keywords: [", fp); fputs("keywords: [", fp);
for (const auto &k : keyword) fmt::print(fp, "'{}', ", k); for (const auto &val : values) fmt::print(fp, "'{}', ", val.keyword);
fputs("]\ndata:\n", fp); fputs("]\ndata:\n", fp);
} }
fmt::print(fp, " {}:\n", ntimestep); fmt::print(fp, " {}:\n", ntimestep);
for (i = 0; i < nrows; i++) { for (int i = 0; i < nrows; i++) {
fputs(" - [", fp); fputs(" - [", fp);
for (j = 0; j < nvalues; j++) fmt::print(fp,"{}, ",array_total[i][j]/norm); for (int j = 0; j < nvalues; j++) fmt::print(fp,"{}, ",array_total[i][j]/norm);
fputs("]\n", fp); fputs("]\n", fp);
} }
} else { } else {
fmt::print(fp,"{} {}\n",ntimestep,nrows); fmt::print(fp,"{} {}\n",ntimestep,nrows);
for (i = 0; i < nrows; i++) { for (int i = 0; i < nrows; i++) {
fprintf(fp,"%d",i+1); fprintf(fp,"%d",i+1);
for (j = 0; j < nvalues; j++) fprintf(fp,format,array_total[i][j]/norm); for (int j = 0; j < nvalues; j++) fprintf(fp,format,array_total[i][j]/norm);
fprintf(fp,"\n"); fprintf(fp,"\n");
} }
} }
@ -938,18 +901,16 @@ int FixAveTime::column_length(int dynamic)
if (!dynamic) { if (!dynamic) {
length = 0; length = 0;
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
if (varlen[i]) continue; if (val.varlen) continue;
if (which[i] == ArgInfo::COMPUTE) { if (val.which == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[i]); if (val.argindex == 0)
if (argindex[i] == 0) lengthone = val.val.c->size_vector;
lengthone = modify->compute[icompute]->size_vector; else lengthone = val.val.c->size_array_rows;
else lengthone = modify->compute[icompute]->size_array_rows; } else if (val.which == ArgInfo::FIX) {
} else if (which[i] == ArgInfo::FIX) { if (val.argindex == 0) lengthone = val.val.f->size_vector;
int ifix = modify->find_fix(ids[i]); else lengthone = val.val.f->size_array_rows;
if (argindex[i] == 0) lengthone = modify->fix[ifix]->size_vector; } else if (val.which == ArgInfo::VARIABLE) {
else lengthone = modify->fix[ifix]->size_array_rows;
} else if (which[i] == ArgInfo::VARIABLE) {
// variables are always varlen = 1, so dynamic // variables are always varlen = 1, so dynamic
} }
if (length == 0) length = lengthone; if (length == 0) length = lengthone;
@ -965,15 +926,13 @@ int FixAveTime::column_length(int dynamic)
if (dynamic) { if (dynamic) {
length = 0; length = 0;
for (int i = 0; i < nvalues; i++) { for (auto &val : values) {
if (varlen[i] == 0) continue; if (val.varlen == 0) continue;
m = value2index[i]; if (val.which == ArgInfo::COMPUTE) {
if (which[i] == ArgInfo::COMPUTE) { lengthone = val.val.c->lock_length();
Compute *compute = modify->compute[m]; } else if (val.which == ArgInfo::VARIABLE) {
lengthone = compute->lock_length();
} else if (which[i] == ArgInfo::VARIABLE) {
double *varvec; double *varvec;
lengthone = input->variable->compute_vector(m,&varvec); lengthone = input->variable->compute_vector(val.val.v,&varvec);
} }
if (mode == SCALAR) continue; if (mode == SCALAR) continue;
if (all_variable_length) { if (all_variable_length) {
@ -1036,7 +995,7 @@ int FixAveTime::modify_param(int narg, char **arg)
int icol = -1; int icol = -1;
if (utils::is_integer(arg[1])) { if (utils::is_integer(arg[1])) {
icol = utils::inumeric(FLERR, arg[1], false, lmp); icol = utils::inumeric(FLERR, arg[1], false, lmp);
if (icol < 0) icol = keyword.size() + icol + 1; if (icol < 0) icol = values.size() + icol + 1;
icol--; icol--;
} else { } else {
try { try {
@ -1045,9 +1004,9 @@ int FixAveTime::modify_param(int narg, char **arg)
icol = -1; icol = -1;
} }
} }
if ((icol < 0) || (icol >= (int) keyword.size())) if ((icol < 0) || (icol >= (int) values.size()))
error->all(FLERR, "Thermo_modify colname column {} invalid", arg[1]); error->all(FLERR, "Thermo_modify colname column {} invalid", arg[1]);
keyword[icol] = arg[2]; values[icol].keyword = arg[2];
return 3; return 3;
} }
return 0; return 0;
@ -1081,7 +1040,7 @@ void FixAveTime::options(int iarg, int narg, char **arg)
if (strcmp(arg[iarg],"file") == 0) { if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/time command"); if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/time command");
yaml_flag = utils::strmatch(arg[iarg+1],"\\.[yY][aA]?[mM][lL]$"); yaml_flag = utils::strmatch(arg[iarg+1],"\\.[yY][aA]?[mM][lL]$");
if (me == 0) { if (comm->me == 0) {
fp = fopen(arg[iarg+1],"w"); fp = fopen(arg[iarg+1],"w");
if (fp == nullptr) if (fp == nullptr)
error->one(FLERR,"Cannot open fix ave/time file {}: {}", error->one(FLERR,"Cannot open fix ave/time file {}: {}",
@ -1140,7 +1099,7 @@ void FixAveTime::options(int iarg, int narg, char **arg)
delete[] title3; delete[] title3;
title3 = utils::strdup(arg[iarg+1]); title3 = utils::strdup(arg[iarg+1]);
iarg += 2; iarg += 2;
} else error->all(FLERR,"Illegal fix ave/time command"); } else error->all(FLERR,"Unknown fix ave/time command option {}", arg[iarg]);
} }
} }

25
src/fix_ave_time.h
View File

@ -40,12 +40,24 @@ class FixAveTime : public Fix {
double compute_array(int, int) override; double compute_array(int, int) override;
private: private:
int me, nvalues; struct value_t {
int nrepeat, nfreq, irepeat; int which; // type of data: COMPUTE, FIX, VARIABLE
int argindex; // 1-based index if data is vector, else 0
int varlen; // 1 if value is from variable-length compute
int offcol;
std::string id; // compute/fix/variable ID
std::string keyword; // column keyword in output
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
std::vector<value_t> values;
int nvalues, nrepeat, nfreq, irepeat;
bigint nvalid, nvalid_last; bigint nvalid, nvalid_last;
int *which, *argindex, *value2index, *offcol;
int *varlen; // 1 if value is from variable-length compute
char **ids;
FILE *fp; FILE *fp;
int nrows; int nrows;
int any_variable_length; int any_variable_length;
@ -61,7 +73,6 @@ class FixAveTime : public Fix {
bigint filepos; bigint filepos;
std::map<std::string, int> key2col; std::map<std::string, int> key2col;
std::vector<std::string> keyword;
int norm, iwindow, window_limit; int norm, iwindow, window_limit;
double *vector; double *vector;
@ -79,8 +90,6 @@ class FixAveTime : public Fix {
void allocate_arrays(); void allocate_arrays();
bigint nextvalid(); bigint nextvalid();
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
#endif #endif
#endif #endif

7
src/fix_vector.cpp
View File

@ -25,9 +25,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum { ONE, RUNNING, WINDOW };
enum { SCALAR, VECTOR };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) : FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
@ -38,6 +35,8 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
nevery = utils::inumeric(FLERR, arg[3], false, lmp); nevery = utils::inumeric(FLERR, arg[3], false, lmp);
if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery); if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery);
// parse values
values.clear(); values.clear();
for (int iarg = 4; iarg < narg; iarg++) { for (int iarg = 4; iarg < narg; iarg++) {
ArgInfo argi(arg[iarg]); ArgInfo argi(arg[iarg]);
@ -48,7 +47,7 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
val.id = argi.get_name(); val.id = argi.get_name();
val.val.c = nullptr; val.val.c = nullptr;
if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_type() == ArgInfo::NONE) || if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) ||
(argi.get_dim() > 1)) (argi.get_dim() > 1))
error->all(FLERR, "Invalid fix vector argument: {}", arg[iarg]); error->all(FLERR, "Invalid fix vector argument: {}", arg[iarg]);

6
src/utils.cpp
View File

@ -1815,7 +1815,7 @@ re_t re_compile(re_ctx_t context, const char *pattern)
/* Private functions: */ /* Private functions: */
static int matchdigit(char c) static int matchdigit(char c)
{ {
return ((c >= '0') && (c <= '9')); return isdigit(c);
} }
static int matchint(char c) static int matchint(char c)
@ -1830,12 +1830,12 @@ static int matchfloat(char c)
static int matchalpha(char c) static int matchalpha(char c)
{ {
return ((c >= 'a') && (c <= 'z')) || ((c >= 'A') && (c <= 'Z')); return isalpha(c);
} }
static int matchwhitespace(char c) static int matchwhitespace(char c)
{ {
return ((c == ' ') || (c == '\t') || (c == '\n') || (c == '\r') || (c == '\f') || (c == '\v')); return isspace(c);
} }
static int matchalphanum(char c) static int matchalphanum(char c)

49
src/variable.cpp
View File

@ -389,19 +389,37 @@ void Variable::set(int narg, char **arg)
// 3rd is filled on retrieval // 3rd is filled on retrieval
} else if (strcmp(arg[1],"format") == 0) { } else if (strcmp(arg[1],"format") == 0) {
constexpr char validfmt[] = "^% ?-?[0-9]*\\.?[0-9]*[efgEFG]$";
if (narg != 4) error->all(FLERR,"Illegal variable command: expected 4 arguments but found {}", narg); if (narg != 4) error->all(FLERR,"Illegal variable command: expected 4 arguments but found {}", narg);
if (find(arg[0]) >= 0) return; int ivar = find(arg[0]);
if (nvar == maxvar) grow(); int jvar = find(arg[2]);
style[nvar] = FORMAT; if (jvar < 0)
num[nvar] = 3; error->all(FLERR, "Variable {}: format variable {} does not exist", arg[0], arg[2]);
which[nvar] = 0; if (!equalstyle(jvar))
pad[nvar] = 0; error->all(FLERR, "Variable {}: format variable {} has incompatible style", arg[0], arg[2]);
if (!utils::strmatch(arg[3],"%[0-9 ]*\\.[0-9]+[efgEFG]")) if (ivar >= 0) {
error->all(FLERR,"Incorrect conversion in format string"); if (style[ivar] != FORMAT)
data[nvar] = new char*[num[nvar]]; error->all(FLERR,"Cannot redefine variable as a different style");
copy(2,&arg[2],data[nvar]); if (!utils::strmatch(arg[3], validfmt))
data[nvar][2] = new char[VALUELENGTH]; error->all(FLERR,"Incorrect conversion in format string");
strcpy(data[nvar][2],"(undefined)"); delete[] data[ivar][0];
delete[] data[ivar][1];
data[ivar][0] = utils::strdup(arg[2]);
data[ivar][1] = utils::strdup(arg[3]);
replaceflag = 1;
} else {
if (nvar == maxvar) grow();
style[nvar] = FORMAT;
num[nvar] = 3;
which[nvar] = 0;
pad[nvar] = 0;
if (!utils::strmatch(arg[3], validfmt))
error->all(FLERR,"Incorrect conversion in format string");
data[nvar] = new char*[num[nvar]];
copy(2,&arg[2],data[nvar]);
data[nvar][2] = new char[VALUELENGTH];
strcpy(data[nvar][2],"(undefined)");
}
// EQUAL // EQUAL
// replace pre-existing var if also style EQUAL (allows it to be reset) // replace pre-existing var if also style EQUAL (allows it to be reset)
@ -940,8 +958,11 @@ char *Variable::retrieve(const char *name)
str = data[ivar][1]; str = data[ivar][1];
} else if (style[ivar] == FORMAT) { } else if (style[ivar] == FORMAT) {
int jvar = find(data[ivar][0]); int jvar = find(data[ivar][0]);
if (jvar == -1) return nullptr; if (jvar < 0)
if (!equalstyle(jvar)) return nullptr; error->all(FLERR, "Variable {}: format variable {} does not exist", names[ivar],data[ivar][0]);
if (!equalstyle(jvar))
error->all(FLERR, "Variable {}: format variable {} has incompatible style",
names[ivar],data[ivar][0]);
double answer = compute_equal(jvar); double answer = compute_equal(jvar);
sprintf(data[ivar][2],data[ivar][1],answer); sprintf(data[ivar][2],data[ivar][1],answer);
str = data[ivar][2]; str = data[ivar][2];

15
unittest/commands/CMakeLists.txt
View File

@ -58,10 +58,17 @@ target_compile_definitions(test_reset_ids PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CU
target_link_libraries(test_reset_ids PRIVATE lammps GTest::GMock) target_link_libraries(test_reset_ids PRIVATE lammps GTest::GMock)
add_test(NAME ResetIDs COMMAND test_reset_ids) add_test(NAME ResetIDs COMMAND test_reset_ids)
add_executable(test_compute_global test_compute_global.cpp) if(PKG_MOLECULE)
target_compile_definitions(test_compute_global PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR}) add_executable(test_compute_global test_compute_global.cpp)
target_link_libraries(test_compute_global PRIVATE lammps GTest::GMock) target_compile_definitions(test_compute_global PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME ComputeGlobal COMMAND test_compute_global) target_link_libraries(test_compute_global PRIVATE lammps GTest::GMock)
add_test(NAME ComputeGlobal COMMAND test_compute_global)
add_executable(test_compute_chunk test_compute_chunk.cpp)
target_compile_definitions(test_compute_chunk PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(test_compute_chunk PRIVATE lammps GTest::GMock)
add_test(NAME ComputeChunk COMMAND test_compute_chunk)
endif()
add_executable(test_mpi_load_balancing test_mpi_load_balancing.cpp) add_executable(test_mpi_load_balancing test_mpi_load_balancing.cpp)
target_link_libraries(test_mpi_load_balancing PRIVATE lammps GTest::GMock) target_link_libraries(test_mpi_load_balancing PRIVATE lammps GTest::GMock)

342
unittest/commands/test_compute_chunk.cpp Normal file
View File

@ -0,0 +1,342 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "../testing/core.h"
#include "info.h"
#include "lammps.h"
#include "library.h"
#include "utils.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include <cstdio>
#include <mpi.h>
// whether to print verbose output (i.e. not capturing LAMMPS screen output).
bool verbose = false;
// we compare floating point numbers with 8 digits precision after the decimal point
static constexpr double EPSILON = 1.0e-8;
namespace LAMMPS_NS {
#define STRINGIFY(val) XSTR(val)
#define XSTR(val) #val
class ComputeChunkTest : public LAMMPSTest {
protected:
void SetUp() override
{
testbinary = "ComputeChunkTest";
LAMMPSTest::SetUp();
if (info->has_style("atom", "full")) {
BEGIN_HIDE_OUTPUT();
command("variable input_dir index \"" STRINGIFY(TEST_INPUT_FOLDER) "\"");
command("include \"${input_dir}/in.fourmol\"");
command("group allwater molecule 3:6");
command("region half block 0.0 INF INF INF INF INF");
command("compute tags all property/atom id");
command("compute bin1d all chunk/atom bin/1d x lower 3.0 units box");
command("compute bin2d all chunk/atom bin/2d x lower 3.0 y lower 3.0 units box");
command("compute bin3d all chunk/atom bin/3d x lower 3.0 y lower 3.0 z lower 3.0 units "
"box");
command("compute binsph all chunk/atom bin/sphere 0.0 0.0 0.0 0.01 6.01 6 units box");
command("compute bincyl all chunk/atom bin/cylinder z lower 3.0 1.0 1.0 0.01 6.01 6 "
"units box");
command("compute mols all chunk/atom molecule");
command("compute types all chunk/atom type");
END_HIDE_OUTPUT();
}
}
double get_scalar(const char *id)
{
return *(double *)lammps_extract_compute(lmp, id, LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR);
}
double *get_vector(const char *id)
{
return (double *)lammps_extract_compute(lmp, id, LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR);
}
double **get_array(const char *id)
{
return (double **)lammps_extract_compute(lmp, id, LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY);
}
double *get_peratom(const char *id)
{
return (double *)lammps_extract_compute(lmp, id, LMP_STYLE_ATOM, LMP_TYPE_VECTOR);
}
};
static constexpr int chunk1d[] = {0, 2, 3, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 3,
4, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 2, 2, 2};
static constexpr int chalf1d[] = {0, 0, 3, 0, 0, 0, 3, 3, 3, 3, 3, 3, 4, 4, 3,
4, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 0, 0, 0};
static constexpr int chunk2d[] = {0, 9, 14, 8, 9, 8, 13, 13, 13, 13, 13, 12, 18, 18, 13,
18, 12, 12, 14, 14, 14, 17, 17, 17, 14, 14, 14, 7, 7, 7};
static constexpr int chunk3d[] = {0, 43, 68, 38, 43, 38, 63, 62, 63, 63, 63, 58, 88, 88, 62,
88, 58, 59, 67, 67, 67, 82, 82, 82, 69, 69, 69, 34, 34, 34};
static constexpr int chunksph[] = {0, 3, 4, 2, 3, 2, 2, 3, 2, 2, 3, 4, 4, 5, 4,
4, 4, 4, 6, 6, 6, 6, 6, 6, 5, 5, 6, 6, 6, 5};
static constexpr int chunkcyl[] = {0, 8, 13, 8, 13, 8, 8, 7, 8, 8, 13, 18, 13, 18, 12,
13, 18, 19, 12, 7, 17, 27, 27, 27, 14, 14, 19, 29, 29, 29};
static constexpr int chunkmol[] = {0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6};
static constexpr int chunktyp[] = {0, 3, 2, 1, 2, 2, 1, 4, 3, 2, 1, 2, 1, 2, 2,
2, 1, 4, 4, 2, 2, 5, 2, 2, 5, 2, 2, 5, 2, 2};
TEST_F(ComputeChunkTest, ChunkAtom)
{
if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
BEGIN_HIDE_OUTPUT();
command("pair_style lj/cut/coul/cut 10.0");
command("pair_coeff * * 0.01 3.0");
command("bond_style harmonic");
command("bond_coeff * 100.0 1.5");
command("dump 1 all custom 1 compute_chunk_atom.lammpstrj "
"id c_bin1d c_bin2d c_bin3d c_binsph c_bincyl c_mols c_types c_tags");
command("run 0 post no");
END_HIDE_OUTPUT();
const int natoms = lammps_get_natoms(lmp);
EXPECT_EQ(get_scalar("bin1d"), 5);
EXPECT_EQ(get_scalar("bin2d"), 25);
EXPECT_EQ(get_scalar("bin3d"), 125);
EXPECT_EQ(get_scalar("binsph"), 6);
EXPECT_EQ(get_scalar("bincyl"), 30);
EXPECT_EQ(get_scalar("mols"), 6);
EXPECT_EQ(get_scalar("types"), 5);
auto cbin1d = get_peratom("bin1d");
auto cbin2d = get_peratom("bin2d");
auto cbin3d = get_peratom("bin3d");
auto cbinsph = get_peratom("binsph");
auto cbincyl = get_peratom("bincyl");
auto cmols = get_peratom("mols");
auto ctypes = get_peratom("types");
auto tag = get_peratom("tags");
for (int i = 0; i < natoms; ++i) {
EXPECT_EQ(cbin1d[i], chunk1d[(int)tag[i]]);
EXPECT_EQ(cbin2d[i], chunk2d[(int)tag[i]]);
EXPECT_EQ(cbin3d[i], chunk3d[(int)tag[i]]);
EXPECT_EQ(cbinsph[i], chunksph[(int)tag[i]]);
EXPECT_EQ(cbincyl[i], chunkcyl[(int)tag[i]]);
EXPECT_EQ(cmols[i], chunkmol[(int)tag[i]]);
EXPECT_EQ(ctypes[i], chunktyp[(int)tag[i]]);
}
BEGIN_HIDE_OUTPUT();
command("uncompute bin1d");
command("compute bin1d all chunk/atom bin/1d x lower 0.2 units reduced region half");
command("uncompute bin3d");
command("compute bin3d all chunk/atom bin/3d x lower 3.0 y lower 3.0 z lower 3.0 "
"compress yes units box");
END_HIDE_OUTPUT();
EXPECT_EQ(get_scalar("bin1d"), 5);
EXPECT_EQ(get_scalar("bin3d"), 12);
cbin1d = get_peratom("bin1d");
tag = get_peratom("tags");
for (int i = 0; i < natoms; ++i) {
EXPECT_EQ(cbin1d[i], chalf1d[(int)tag[i]]);
}
// cleanup
platform::unlink("compute_chunk_atom.lammpstrj");
}
TEST_F(ComputeChunkTest, PropertyChunk)
{
if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
BEGIN_HIDE_OUTPUT();
command("pair_style lj/cut/coul/cut 10.0");
command("pair_coeff * * 0.01 3.0");
command("bond_style harmonic");
command("bond_coeff * 100.0 1.5");
command("uncompute bin3d");
command("compute bin3d all chunk/atom bin/3d x lower 3.0 y lower 3.0 z lower 3.0 "
"compress yes units box");
command("compute prop1 all property/chunk bin1d count");
command("compute prop2 all property/chunk bin2d count");
command("compute prop3 all property/chunk bin3d id count");
command("fix hist1 all ave/time 1 1 1 c_prop1 mode vector");
command("fix hist2 all ave/time 1 1 1 c_prop2 mode vector");
command("fix hist3 all ave/time 1 1 1 c_prop3[*] mode vector");
command("run 0 post no");
END_HIDE_OUTPUT();
auto cprop1 = get_vector("prop1");
EXPECT_EQ(cprop1[0], 0);
EXPECT_EQ(cprop1[1], 7);
EXPECT_EQ(cprop1[2], 16);
EXPECT_EQ(cprop1[3], 6);
EXPECT_EQ(cprop1[4], 0);
auto cprop2 = get_vector("prop2");
int nempty = 0;
int ncount = 0;
for (int i = 0; i < 25; ++i) {
if (cprop2[i] == 0)
++nempty;
else
ncount += cprop2[i];
}
EXPECT_EQ(nempty, 17);
EXPECT_EQ(ncount, 29);
auto cprop3 = get_array("prop3");
EXPECT_EQ(cprop3[0][0], 34);
EXPECT_EQ(cprop3[1][0], 38);
EXPECT_EQ(cprop3[2][0], 43);
EXPECT_EQ(cprop3[3][0], 58);
EXPECT_EQ(cprop3[4][0], 59);
EXPECT_EQ(cprop3[5][0], 62);
EXPECT_EQ(cprop3[6][0], 63);
EXPECT_EQ(cprop3[7][0], 67);
EXPECT_EQ(cprop3[8][0], 68);
EXPECT_EQ(cprop3[9][0], 69);
EXPECT_EQ(cprop3[10][0], 82);
EXPECT_EQ(cprop3[11][0], 88);
EXPECT_EQ(cprop3[0][1], 3);
EXPECT_EQ(cprop3[1][1], 2);
EXPECT_EQ(cprop3[2][1], 2);
EXPECT_EQ(cprop3[3][1], 2);
EXPECT_EQ(cprop3[4][1], 1);
EXPECT_EQ(cprop3[5][1], 2);
EXPECT_EQ(cprop3[6][1], 4);
EXPECT_EQ(cprop3[7][1], 3);
EXPECT_EQ(cprop3[8][1], 1);
EXPECT_EQ(cprop3[9][1], 3);
EXPECT_EQ(cprop3[10][1], 3);
EXPECT_EQ(cprop3[11][1], 3);
}
TEST_F(ComputeChunkTest, ChunkComputes)
{
if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
BEGIN_HIDE_OUTPUT();
command("pair_style lj/cut/coul/cut 10.0");
command("pair_coeff * * 0.01 3.0");
command("bond_style harmonic");
command("bond_coeff * 100.0 1.5");
command("compute ang all angmom/chunk mols");
command("compute com all com/chunk mols");
command("compute dip all dipole/chunk mols geometry");
command("compute gyr all gyration/chunk mols");
command("compute mom all inertia/chunk mols");
command("compute omg all omega/chunk mols");
command("compute tmp all temp/chunk mols com yes");
command("compute trq all torque/chunk mols");
command("compute vcm all vcm/chunk mols");
command("fix hist1 all ave/time 1 1 1 c_ang[*] c_com[*] c_dip[*] c_gyr c_mom[*] c_omg[*] "
"c_trq[*] c_vcm[*] mode vector file all_chunk.dat format %15.8f");
command("fix hist2 all ave/time 1 1 1 c_tmp mode vector file vec_chunk.dat format %15.8f");
command("run 0 post no");
END_HIDE_OUTPUT();
auto cang = get_array("ang");
auto ccom = get_array("com");
auto cdip = get_array("dip");
auto cgyr = get_vector("gyr");
auto cmom = get_array("mom");
auto comg = get_array("omg");
auto ctmp = get_vector("tmp");
auto ctrq = get_array("trq");
auto cvcm = get_array("vcm");
EXPECT_NEAR(cang[0][0], -0.01906982, EPSILON);
EXPECT_NEAR(cang[0][1], -0.02814532, EPSILON);
EXPECT_NEAR(cang[0][2], -0.03357393, EPSILON);
EXPECT_NEAR(cang[5][0], 0.00767837, EPSILON);
EXPECT_NEAR(cang[5][1], -0.00303138, EPSILON);
EXPECT_NEAR(cang[5][2], 0.00740977, EPSILON);
EXPECT_NEAR(ccom[1][0], 2.27051374, EPSILON);
EXPECT_NEAR(ccom[1][1], -1.21038876, EPSILON);
EXPECT_NEAR(ccom[1][2], -0.58581655, EPSILON);
EXPECT_NEAR(ccom[5][0], -1.98284693, EPSILON);
EXPECT_NEAR(ccom[5][1], -4.17351226, EPSILON);
EXPECT_NEAR(ccom[5][2], 2.04850072, EPSILON);
EXPECT_NEAR(cmom[2][0], 4.28810281, EPSILON);
EXPECT_NEAR(cmom[2][1], 4.99562488, EPSILON);
EXPECT_NEAR(cmom[2][2], 5.34954800, EPSILON);
EXPECT_NEAR(cmom[5][0], 3.06867233, EPSILON);
EXPECT_NEAR(cmom[5][1], 5.24202887, EPSILON);
EXPECT_NEAR(cmom[5][2], 6.06478557, EPSILON);
EXPECT_NEAR(comg[3][0], -0.00349423, EPSILON);
EXPECT_NEAR(comg[3][1], -0.00025062, EPSILON);
EXPECT_NEAR(comg[3][2], -0.00323573, EPSILON);
EXPECT_NEAR(comg[5][0], 0.00437315, EPSILON);
EXPECT_NEAR(comg[5][1], 0.00029335, EPSILON);
EXPECT_NEAR(comg[5][2], 0.00268517, EPSILON);
EXPECT_NEAR(ctrq[4][0], -0.94086086, EPSILON);
EXPECT_NEAR(ctrq[4][1], 0.56227336, EPSILON);
EXPECT_NEAR(ctrq[4][2], 0.75139995, EPSILON);
EXPECT_NEAR(ctrq[5][0], -0.07066910, EPSILON);
EXPECT_NEAR(ctrq[5][1], -0.58556032, EPSILON);
EXPECT_NEAR(ctrq[5][2], -0.81513604, EPSILON);
EXPECT_NEAR(cvcm[0][0], -0.00011274, EPSILON);
EXPECT_NEAR(cvcm[0][1], 0.00071452, EPSILON);
EXPECT_NEAR(cvcm[0][2], -0.00017908, EPSILON);
EXPECT_NEAR(cvcm[5][0], -0.00063326, EPSILON);
EXPECT_NEAR(cvcm[5][1], 0.00007092, EPSILON);
EXPECT_NEAR(cvcm[5][2], 0.00045545, EPSILON);
EXPECT_NEAR(cdip[0][3], 0.35912150, EPSILON);
EXPECT_NEAR(cdip[1][3], 0.68453713, EPSILON);
EXPECT_NEAR(cdip[2][3], 0.50272643, EPSILON);
EXPECT_NEAR(cdip[3][3], 0.50845862, EPSILON);
EXPECT_NEAR(cdip[4][3], 0.49757365, EPSILON);
EXPECT_NEAR(cdip[5][3], 0.49105019, EPSILON);
EXPECT_NEAR(cgyr[0], 1.48351858, EPSILON);
EXPECT_NEAR(cgyr[1], 1.56649567, EPSILON);
EXPECT_NEAR(cgyr[2], 0.55196552, EPSILON);
EXPECT_NEAR(cgyr[3], 0.54573649, EPSILON);
EXPECT_NEAR(cgyr[4], 0.54793875, EPSILON);
EXPECT_NEAR(cgyr[5], 0.54708204, EPSILON);
EXPECT_NEAR(ctmp[0], 1.08268576, EPSILON);
EXPECT_NEAR(ctmp[1], 1.61905718, EPSILON);
EXPECT_NEAR(ctmp[2], 1.41991778, EPSILON);
EXPECT_NEAR(ctmp[3], 0.55484671, EPSILON);
EXPECT_NEAR(ctmp[4], -0.06062938, EPSILON);
EXPECT_NEAR(ctmp[5], -0.09219489, EPSILON);
}
} // namespace LAMMPS_NS
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
::testing::InitGoogleMock(&argc, argv);
if (LAMMPS_NS::platform::mpi_vendor() == "Open MPI" && !Info::has_exceptions())
std::cout << "Warning: using OpenMPI without exceptions. Death tests will be skipped\n";
// handle arguments passed via environment variable
if (const char *var = getenv("TEST_ARGS")) {
std::vector<std::string> env = LAMMPS_NS::utils::split_words(var);
for (auto arg : env) {
if (arg == "-v") {
verbose = true;
}
}
}
if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;
int rv = RUN_ALL_TESTS();
MPI_Finalize();
return rv;
}

85
unittest/commands/test_variables.cpp
View File

@ -219,7 +219,7 @@ TEST_F(VariableTest, CreateDelete)
TEST_FAILURE(".*ERROR: All universe/uloop variables must have same # of values.*", TEST_FAILURE(".*ERROR: All universe/uloop variables must have same # of values.*",
command("variable ten4 uloop 2");); command("variable ten4 uloop 2"););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*", TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable ten11 format two \"%08f\"");); command("variable ten11 format two \"%08x\""););
TEST_FAILURE(".*ERROR: Variable name 'ten@12' must have only letters, numbers, or undersc.*", TEST_FAILURE(".*ERROR: Variable name 'ten@12' must have only letters, numbers, or undersc.*",
command("variable ten@12 index one two three");); command("variable ten@12 index one two three"););
TEST_FAILURE(".*ERROR: Variable evaluation before simulation box is defined.*", TEST_FAILURE(".*ERROR: Variable evaluation before simulation box is defined.*",
@ -281,7 +281,7 @@ TEST_F(VariableTest, AtomicSystem)
TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*", TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*",
command("variable one atom x");); command("variable one atom x"););
TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*", TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*",
command("variable one vector f_press");); command("variable id vector f_press"););
TEST_FAILURE(".*ERROR on proc 0: Cannot open file variable file test_variable.xxx.*", TEST_FAILURE(".*ERROR on proc 0: Cannot open file variable file test_variable.xxx.*",
command("variable ten1 atomfile test_variable.xxx");); command("variable ten1 atomfile test_variable.xxx"););
TEST_FAILURE(".*ERROR: Variable loop: has a circular dependency.*", TEST_FAILURE(".*ERROR: Variable loop: has a circular dependency.*",
@ -638,6 +638,87 @@ TEST_F(VariableTest, Label2TypeMolecular)
ASSERT_THAT(variable->retrieve("d1"), StrEq("1")); ASSERT_THAT(variable->retrieve("d1"), StrEq("1"));
ASSERT_THAT(variable->retrieve("i1"), StrEq("1")); ASSERT_THAT(variable->retrieve("i1"), StrEq("1"));
} }
TEST_F(VariableTest, Format)
{
BEGIN_HIDE_OUTPUT();
command("variable idx index -0.625");
command("variable one equal -0.625");
command("variable two equal 1.0e-20");
command("variable three equal 1.0e10");
command("variable f1one format one \"%8.4f\"");
command("variable f1two format two %8.4f");
command("variable f2one format one %.2F");
command("variable f2two format two \"% .25F\"");
command("variable f3one format one \"%5f\"");
command("variable f3two format two %f");
command("variable e1one format one \"%14.4e\"");
command("variable e1two format two %-14.4e");
command("variable e2one format one %.2E");
command("variable e2two format two \"% .15E\"");
command("variable e3one format one \"%5e\"");
command("variable e3two format two %e");
command("variable g1one format one %14.4g");
command("variable g1two format two \"%-14.4g\"");
command("variable g2one format one %.2G");
command("variable g2two format two \"% .15G\"");
command("variable g3one format one \"%5g\"");
command("variable g3two format two \"%g\"");
END_HIDE_OUTPUT();
EXPECT_THAT(variable->retrieve("one"), StrEq("-0.625"));
EXPECT_THAT(variable->retrieve("two"), StrEq("1e-20"));
EXPECT_THAT(variable->retrieve("f1one"), StrEq(" -0.6250"));
EXPECT_THAT(variable->retrieve("f1two"), StrEq(" 0.0000"));
EXPECT_THAT(variable->retrieve("f2one"), StrEq("-0.62"));
EXPECT_THAT(variable->retrieve("f2two"), StrEq(" 0.0000000000000000000100000"));
EXPECT_THAT(variable->retrieve("f3one"), StrEq("-0.625000"));
EXPECT_THAT(variable->retrieve("f3two"), StrEq("0.000000"));
EXPECT_THAT(variable->retrieve("e1one"), StrEq(" -6.2500e-01"));
EXPECT_THAT(variable->retrieve("e1two"), StrEq("1.0000e-20 "));
EXPECT_THAT(variable->retrieve("e2one"), StrEq("-6.25E-01"));
EXPECT_THAT(variable->retrieve("e2two"), StrEq(" 9.999999999999999E-21"));
EXPECT_THAT(variable->retrieve("e3one"), StrEq("-6.250000e-01"));
EXPECT_THAT(variable->retrieve("e3two"), StrEq("1.000000e-20"));
EXPECT_THAT(variable->retrieve("g1one"), StrEq(" -0.625"));
EXPECT_THAT(variable->retrieve("g1two"), StrEq("1e-20 "));
EXPECT_THAT(variable->retrieve("g2one"), StrEq("-0.62"));
EXPECT_THAT(variable->retrieve("g2two"), StrEq(" 1E-20"));
EXPECT_THAT(variable->retrieve("g3one"), StrEq("-0.625"));
EXPECT_THAT(variable->retrieve("g3two"), StrEq("1e-20"));
BEGIN_HIDE_OUTPUT();
command("variable f1one format one \"%-8.4f\"");
command("variable two delete");
command("variable two index 12.5");
command("variable f1three format three %g");
command("variable three delete");
END_HIDE_OUTPUT();
EXPECT_THAT(variable->retrieve("f1one"), StrEq("-0.6250 "));
TEST_FAILURE(".*ERROR: Variable f1idx: format variable idx has incompatible style.*",
command("variable f1idx format idx %8.4f"););
TEST_FAILURE(".*ERROR: Variable f1two: format variable two has incompatible style.*",
variable->retrieve("f1two"););
TEST_FAILURE(".*ERROR: Variable f1idx: format variable yyy does not exist.*",
command("variable f1idx format yyy %8.4f"););
TEST_FAILURE(".*ERROR: Variable f1three: format variable three does not exist.*",
variable->retrieve("f1three"););
TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*",
command("variable f2one equal 0.5"););
TEST_FAILURE(".*ERROR: Illegal variable command.*", command("variable xxx format \"xxx\""););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable xxx format one \"xxx\""););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable xxx format one \"%d\""););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable xxx format one \"%g%g\""););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable xxx format one \"%g%5\""););
TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
command("variable xxx format one \"%g%%\""););
// TEST_FAILURE(".*ERROR: Incorrect conversion in format string.*",
// command("print \"${f1idx}\""););
}
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
int main(int argc, char **argv) int main(int argc, char **argv)