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

Merge remote-tracking branch 'github/develop' into fortran-further-tinkering

This commit is contained in:
Axel Kohlmeyer 2022-10-03 09:37:52 -04:00
parent 658df50e97 2970a2140f
commit d711a51657
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
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.
This runs the command ``env LC_ALL=C grep -n '[^ -~]' src/*.rst`` and
thus prints all offending lines with filename and line number
prepended to the screen. Special characters like the Angstrom
:math:`\mathrm{\mathring{A}}` should be typeset with embedded math
(like this ``:math:`\mathrm{\mathring{A}}```\ ).
prepended to the screen. Special characters like greek letters
(:math:`\alpha~~\sigma~~\epsilon`), super- or subscripts
(: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
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
INTEL and KOKKOS packages require you to choose one of their
hardware options when building for a specific platform. I.e. CPU or
Phi option for the INTEL package. Or the OpenMP, Cuda, or Phi
option for the KOKKOS package.
Phi option for the INTEL package. Or the OpenMP, CUDA, HIP, SYCL,
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:
* both the INTEL Phi and KOKKOS Phi options
* 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
performance results for the various accelerator packages for several
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::
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
prefactor. The anharmonic prefactors have units angstrom\^(-n):
-2.55 angstrom\^(-1) and (7/12)2.55\^2 angstrom\^(-2). The code takes
prefactor. The anharmonic prefactors have units :math:`\AA^{-n}`:
: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.
Note that the MM3 papers contains an error in Eq (1):
(7/12)2.55 should be replaced with (7/12)2.55\^2
Note that the MM3 papers contain an error in Eq (1):
: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
: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::
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.
due to a polymer being stretched). The :math:`\sigma` and :math:`\epsilon` used in the
LJ portion of the formula are both set equal to 1.0 by LAMMPS.
due to a polymer being stretched). The :math:`\sigma` and
: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
: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)
This potential was constructed to mimic the FENE bond potential for
coarse-grained polymer chains. When monomers with :math:`\sigma = \epsilon = 1.0`
are used, the following choice of parameters gives a quartic potential that
looks nearly like the FENE potential:
coarse-grained polymer chains. When monomers with :math:`\sigma =
\epsilon = 1.0` are used, the following choice of parameters gives a
quartic potential that looks nearly like the FENE potential:
.. 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.
A *box* value selects standard distance units as defined by the
:doc:`units <units>` command (e.g., :math:`\mathrm{\mathring A}`
for units = *real* or *metal*).
:doc:`units <units>` command (e.g., :math:`\AA` for units = *real* or *metal*).
A *lattice* value means the distance units are in lattice spacings.
The :doc:`lattice <lattice>` command must have been previously used to
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
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
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.
Here are typical input parameters for fcc aluminum (lattice
constant :math:`4.05~\mathrm{\mathring A}`),
constant :math:`4.05~\AA`),
.. parsed-literal::
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::

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
mass (COM) of the body. The *x*, *y*, *z* attributes write the COM
"unscaled", in the appropriate distance :doc:`units <units>`
(:math:`\mathrm{\mathring A}`,
sigma, etc). Use *xu*, *yu*, *zu* if you want the COM "unwrapped" by
(:math:`\AA`, :math:`\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
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

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
good for comparing diffraction results from multiple simulations; however
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
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
for coordinates (*clo*, *chi*) and velocities (*vlo*, *vhi*). A *box* value
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.,
velocity in lattice spacings per unit time). The :doc:`lattice <lattice>`
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
simulations; however, it can reduce the likelihood that Bragg
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.
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*
style, or the overlap distance *Doverlap* by the *overlap* keyword. A
*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
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*
style are determined by the setting of *box* or *lattice* for 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
spacings. The :doc:`lattice <lattice>` command must have been previously used
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*,
*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
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

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`
components of the displacement. The :doc:`dump_modify thresh <dump_modify>`
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
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
@ -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
(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
(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
*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
coordinates in nanometers and time in picoseconds. For example, for *real*
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
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

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
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
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.
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
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>`
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
to define a wall position, but only when a numeric constant is used.
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.
A *lattice* value means the distance units are in lattice spacings.
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
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*,
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
volume considerably smaller/larger than you expect your system to sample
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
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
coefficients should all be in atmospheres \* cubic angstroms.
are using *units real*, *V_avg* should be in cubic Angstroms, and the
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
default in LAMMPS. The *fscale* keyword allows to apply a scaling
factor to forces transmitted by the IMD client. The IMD protocols
stipulates that forces are transferred in kcal/mol/angstrom under the
assumption that coordinates are given in angstrom. For LAMMPS runs
stipulates that forces are transferred in kcal/mol/Angstrom under the
assumption that coordinates are given in Angstrom. For LAMMPS runs
with different units or as a measure to tweak the forces generated by
the manipulation of the IMD client, this option allows to make
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
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>`
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
low accuracy forces, about 1e-4 relative error or less accurate,
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.
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
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
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
the attraction between the I,J atoms based on the position of other
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)
the cutoff is obtained by multiplying the scale factor by the sigma
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
would be 10.2 Angstroms.
would be :math:`10.2 \AA`.
By default, the longer-ranged interaction is smoothly switched off
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
to the CHARMM value, as if it were effectively a parameter of the
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
the CHARMM force field energies and forces, when using one of these
two CHARMM pair styles.

16
doc/src/pair_smtbq.rst
View File

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

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
risking energy drift. Smooth switching of forces between the levels
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
settings for use of the *respa* style without SHAKE in biomolecular
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
speedup can be obtained while maintaining good energy conservation.
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:
.. 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.
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.
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
* distance = Angstroms
* time = femtoseconds
* energy = Kcal/mole
* energy = kcal/mol
* velocity = Angstroms/femtosecond
* force = Kcal/mole-Angstrom
* torque = Kcal/mole
* force = (kcal/mol)/Angstrom
* torque = kcal/mol
* temperature = Kelvin
* pressure = atmospheres
* dynamic viscosity = Poise
* charge = multiple of electron charge (1.0 is a proton)
* dipole = charge\*Angstroms
* electric field = volts/Angstrom
* density = gram/cm\^dim
* density = g/cm\^dim
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
string parameter is quoted.
For the *format* style, an equal-style variable is specified along
with a C-style format string, e.g. "%f" or "%.10g", which must be
appropriate for formatting a double-precision floating-point value.
The default format is "%.15g". This variable style allows an
equal-style variable to be formatted precisely when it is evaluated.
For the *format* style, an equal-style or compatible variable is
specified along with a C-style format string, e.g. "%f" or "%.10g",
which must be appropriate for formatting a double-precision
floating-point value and may not have extra characters. The default
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
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.
#keep_warnings = False
# -- Options for HTML output ----------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
@ -264,8 +263,18 @@ else:
html_math_renderer = 'mathjax'
# 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'
# hack to enable use of \AA in :math:
rst_prolog = r"""
.. only:: html
:math:`\renewcommand{\AA}{\text{Å}}`
"""
# -- Options for LaTeX output ---------------------------------------------
latex_elements = {
@ -278,6 +287,7 @@ latex_elements = {
# Additional stuff for the LaTeX preamble.
'preamble': r'''
\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
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 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(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
idchunk(nullptr), ids(nullptr), which(nullptr), argindex(nullptr), value2index(nullptr)
Compute(lmp, narg, arg), idchunk(nullptr)
{
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
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
nvalues = 0;
for (iarg = 0; iarg < nargnew; iarg++) {
ids[nvalues] = nullptr;
values.clear();
for (iarg = 0; iarg < nargnew; ++iarg) {
ArgInfo argi(arg[iarg], ArgInfo::COMPUTE|ArgInfo::FIX);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
value_t val;
val.which = argi.get_type();
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))
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 (expand) {
for (int i = 0; i < nargnew; i++) delete [] earg[i];
for (int i = 0; i < nargnew; i++) delete[] earg[i];
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 fix, assume a global vector or array is per-chunk data
for (int i = 0; i < nvalues; i++) {
if (which[i] == ArgInfo::COMPUTE) {
auto icompute = modify->get_compute_by_id(ids[i]);
for (auto &val : values) {
if (val.which == ArgInfo::COMPUTE) {
auto icompute = modify->get_compute_by_id(val.id);
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$"))
error->all(FLERR,"Compute for compute chunk/spread/atom "
"does not calculate per-chunk values");
error->all(FLERR,"Compute chunk/spread/atom compute {} does not calculate per-chunk values",
val.id);
if (argindex[i] == 0) {
if (val.argindex == 0) {
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 {
if (!icompute->array_flag)
error->all(FLERR,"Compute chunk/spread/atom compute does not calculate global array");
if (argindex[i] > icompute->size_array_cols)
error->all(FLERR,"Compute chunk/spread/atom compute array is accessed out-of-range");
error->all(FLERR,"Compute chunk/spread/atom compute {} does not calculate global array",
val.id);
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) {
auto ifix = modify->get_fix_by_id(ids[i]);
} else if (val.which == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(val.id);
if (ifix)
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", ids[i]);
if (argindex[i] == 0) {
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", val.id);
if (val.argindex == 0) {
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 {
if (!ifix->array_flag)
error->all(FLERR,"Compute chunk/spread/atom fix does not calculate global array");
if (argindex[i] > ifix->size_array_cols)
error->all(FLERR,"Compute chunk/spread/atom fix array is accessed out-of-range");
error->all(FLERR,"Compute chunk/spread/atom {} fix does not calculate global array",
val.id);
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
peratom_flag = 1;
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
if (values.size() == 1) size_peratom_cols = 0;
else size_peratom_cols = values.size();
// per-atom vector or array
@ -140,13 +141,7 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
ComputeChunkSpreadAtom::~ComputeChunkSpreadAtom()
{
delete [] idchunk;
delete [] which;
delete [] argindex;
for (int i = 0; i < nvalues; i++) delete [] ids[i];
delete [] ids;
delete [] value2index;
delete[] idchunk;
memory->destroy(vector_atom);
memory->destroy(array_atom);
@ -160,18 +155,16 @@ void ComputeChunkSpreadAtom::init()
// set indices of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", ids[m]);
value2index[m] = icompute;
for (auto &val : values) {
if (val.which == ArgInfo::COMPUTE) {
val.val.c = modify->get_compute_by_id(val.id);
if (!val.val.c)
error->all(FLERR,"Compute ID {} for compute chunk/spread/atom does not exist", val.id);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", ids[m]);
value2index[m] = ifix;
} else if (val.which == ArgInfo::FIX) {
val.val.f = modify->get_fix_by_id(val.id);
if (!val.val.f)
error->all(FLERR,"Fix ID {} for compute chunk/spread/atom does not exist", val.id);
}
}
}
@ -182,7 +175,8 @@ void ComputeChunkSpreadAtom::init_chunk()
{
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk));
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)
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
if (atom->nmax > nmax) {
if (nvalues == 1) {
if (values.size() == 1) {
memory->destroy(vector_atom);
nmax = atom->nmax;
memory->create(vector_atom,nmax,"chunk/spread/atom:vector_atom");
} else {
memory->destroy(array_atom);
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 nlocal = atom->nlocal;
int i,m,n,index,nstride;
int index,nstride;
double *ptr;
for (m = 0; m < nvalues; m++) {
n = value2index[m];
int m = 0;
for (auto &val : values) {
// copy compute/fix values into vector_atom or array_atom
// nstride between values for each atom
if (nvalues == 1) {
if (values.size() == 1) {
ptr = vector_atom;
nstride = 1;
} else {
ptr = &array_atom[0][m];
nstride = nvalues;
nstride = values.size();
}
// invoke compute if not previously invoked
if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n];
if (val.which == ArgInfo::COMPUTE) {
Compute *compute = val.val.c;
if (argindex[m] == 0) {
if (val.argindex == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_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;
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -262,9 +256,9 @@ void ComputeChunkSpreadAtom::compute_peratom()
compute->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY;
}
int icol = argindex[m]-1;
int icol = val.argindex-1;
double **carray = compute->array;
for (i = 0; i < nlocal; i++, ptr += nstride) {
for (int i = 0; i < nlocal; i++, ptr += nstride) {
*ptr = 0.0;
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -277,15 +271,15 @@ void ComputeChunkSpreadAtom::compute_peratom()
// are assuming the fix global vector/array is per-chunk data
// check if index exceeds fix output length/rows
} else if (which[m] == ArgInfo::FIX) {
auto &fix = modify->get_fix_list()[n];
} else if (val.which == ArgInfo::FIX) {
Fix *fix = val.val.f;
if (update->ntimestep % fix->global_freq)
error->all(FLERR,"Fix used in compute chunk/spread/atom not "
"computed at compatible time");
error->all(FLERR,"Fix {} used in compute chunk/spread/atom not computed at compatible time",
val.id);
if (argindex[m] == 0) {
if (val.argindex == 0) {
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;
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -294,9 +288,9 @@ void ComputeChunkSpreadAtom::compute_peratom()
}
} else {
int icol = argindex[m]-1;
int icol = val.argindex-1;
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;
if (!(mask[i] & groupbit)) continue;
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 bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double)nmax*values.size() * sizeof(double);
bytes += values.size() * sizeof(value_t);
return bytes;
}

17
src/compute_chunk_spread_atom.h
View File

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

633
src/fix_ave_time.cpp
View File

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

7
src/fix_vector.cpp
View File

@ -25,9 +25,6 @@
using namespace LAMMPS_NS;
using namespace FixConst;
enum { ONE, RUNNING, WINDOW };
enum { SCALAR, VECTOR };
/* ---------------------------------------------------------------------- */
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);
if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery);
// parse values
values.clear();
for (int iarg = 4; iarg < narg; iarg++) {
ArgInfo argi(arg[iarg]);
@ -48,7 +47,7 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
val.id = argi.get_name();
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))
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: */
static int matchdigit(char c)
{
return ((c >= '0') && (c <= '9'));
return isdigit(c);
}
static int matchint(char c)
@ -1830,12 +1830,12 @@ static int matchfloat(char c)
static int matchalpha(char c)
{
return ((c >= 'a') && (c <= 'z')) || ((c >= 'A') && (c <= 'Z'));
return isalpha(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)

49
src/variable.cpp
View File

@ -389,19 +389,37 @@ void Variable::set(int narg, char **arg)
// 3rd is filled on retrieval
} 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 (find(arg[0]) >= 0) return;
if (nvar == maxvar) grow();
style[nvar] = FORMAT;
num[nvar] = 3;
which[nvar] = 0;
pad[nvar] = 0;
if (!utils::strmatch(arg[3],"%[0-9 ]*\\.[0-9]+[efgEFG]"))
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)");
int ivar = find(arg[0]);
int jvar = find(arg[2]);
if (jvar < 0)
error->all(FLERR, "Variable {}: format variable {} does not exist", arg[0], arg[2]);
if (!equalstyle(jvar))
error->all(FLERR, "Variable {}: format variable {} has incompatible style", arg[0], arg[2]);
if (ivar >= 0) {
if (style[ivar] != FORMAT)
error->all(FLERR,"Cannot redefine variable as a different style");
if (!utils::strmatch(arg[3], validfmt))
error->all(FLERR,"Incorrect conversion in format string");
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
// 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];
} else if (style[ivar] == FORMAT) {
int jvar = find(data[ivar][0]);
if (jvar == -1) return nullptr;
if (!equalstyle(jvar)) return nullptr;
if (jvar < 0)
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);
sprintf(data[ivar][2],data[ivar][1],answer);
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)
add_test(NAME ResetIDs COMMAND test_reset_ids)
add_executable(test_compute_global test_compute_global.cpp)
target_compile_definitions(test_compute_global PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(test_compute_global PRIVATE lammps GTest::GMock)
add_test(NAME ComputeGlobal COMMAND test_compute_global)
if(PKG_MOLECULE)
add_executable(test_compute_global test_compute_global.cpp)
target_compile_definitions(test_compute_global PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
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)
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.*",
command("variable ten4 uloop 2"););
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.*",
command("variable ten@12 index one two three"););
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.*",
command("variable one atom x"););
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.*",
command("variable ten1 atomfile test_variable.xxx"););
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("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
int main(int argc, char **argv)