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

Merge branch 'develop' into create-atoms-mesh

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-05-07 09:33:40 -04:00
parent 71185b9f98 4737b9efb7
commit 102b61ae1b
36 changed files with 598 additions and 474 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Build_package.rst
View File

@ -220,7 +220,8 @@ These commands install/un-install sets of packages:
.. code-block:: bash
make yes-all # install all packages
make no-all # uninstall all packages
make no-all # check for changes and uninstall all packages
make no-installed # only check and uninstall installed packages
make yes-basic # install a few commonly used packages'
make no-basic # remove a few commonly used packages'
make yes-most # install most packages w/o libs'

19
doc/src/Build_windows.rst
View File

@ -99,6 +99,10 @@ It is possible to use both the integrated CMake support of the Visual
Studio IDE or use an external CMake installation (e.g. downloaded from
cmake.org) to create build files and compile LAMMPS from the command line.
Compilation via command line and unit tests are checked automatically
for the LAMMPS development branch through
`GitHub Actions <https://github.com/lammps/lammps/actions/workflows/compile-msvc.yml>`_.
.. note::
Versions of Visual Studio before version 17.1 may scan the entire
@ -111,6 +115,10 @@ Please note, that for either approach CMake will create a so-called
the command lines for building and testing LAMMPS must be adjusted
accordingly.
The LAMMPS cmake folder contains a ``CMakeSettings.json`` file with
build configurations for MSVC compilers and the MS provided Clang
compiler package in Debug and Release mode.
To support running in parallel you can compile with OpenMP enabled using
the OPENMP package or install Microsoft MPI (including the SDK) and compile
LAMMPS with MPI enabled.
@ -145,14 +153,3 @@ LAMMPS developers. We instead rely on the feedback of the users
of these pre-compiled LAMMPS packages for Windows. We will try to resolve
issues to the best of our abilities if we become aware of them. However
this is subject to time constraints and focus on HPC platforms.
.. _native:
Native Visual C++ support
^^^^^^^^^^^^^^^^^^^^^^^^^
Support for the Visual C++ compilers is currently not available. The
CMake build system is capable of creating suitable a Visual Studio
style build environment, but the LAMMPS source code itself is not
ported to fully support Visual C++. Volunteers to take on this task
are welcome.

15
doc/src/Commands_parse.rst
View File

@ -77,18 +77,19 @@ LAMMPS:
so that you do not have to define (or discard) a temporary variable,
"X" in this case.
Additionally, the "immediate" variable expression may be followed by
a colon, followed by a C-style format string, e.g. ":%f" or ":%.10g".
The format string must be appropriate for a double-precision
floating-point value. The format string is used to output the result
of the variable expression evaluation. If a format string is not
specified a high-precision "%.20g" is used as the default.
Additionally, the entire "immediate" variable expression may be
followed by a colon, followed by a C-style format string,
e.g. ":%f" or ":%.10g". The format string must be appropriate for
a double-precision floating-point value. The format string is used
to output the result of the variable expression evaluation. If a
format string is not specified, a high-precision "%.20g" is used as
the default format.
This can be useful for formatting print output to a desired precision:
.. code-block:: LAMMPS
print "Final energy per atom: $(pe/atoms:%10.3f) eV/atom"
print "Final energy per atom: $(v_ke_per_atom+v_pe_per_atom:%10.3f) eV/atom"
Note that neither the curly-bracket or immediate form of variables
can contain nested $ characters for other variables to substitute

60
doc/src/fix_print.rst
View File

@ -39,19 +39,35 @@ Description
Print a text string every N steps during a simulation run. This can
be used for diagnostic purposes or as a debugging tool to monitor some
quantity during a run. The text string must be a single argument, so
it should be enclosed in double quotes if it is more than one word.
If it contains variables it must be enclosed in double quotes to
insure they are not evaluated when the input script line is read, but
will instead be evaluated each time the string is printed.
it should be enclosed in single or double quotes if it is more than
one word. If it contains variables it must be enclosed in double
quotes to insure they are not evaluated when the input script line is
read, but will instead be evaluated each time the string is printed.
Instead of a numeric value, N can be specified as an :doc:`equal-style variable <variable>`, which should be specified as v_name, where
name is the variable name. In this case, the variable is evaluated at
the beginning of a run to determine the **next** timestep at which the
.. note::
As discussed on the :doc:`Commands parse <Commands_parse>` doc
page, the text string can use "immediate" variables, specified as
$(formula) with parenthesis, where the numeric formula has the same
syntax as equal-style variables described on the :doc:`variable
<variable>` doc page. This is a convenient way to evaluate a
formula immediately without using the variable command to define a
named variable and then use that variable in the text string. The
formula can include a trailing colon and format string which
determines the precision with which the numeric value is output.
This is also explained on the :doc:`Commands parse
<Commands_parse>` doc page.
Instead of a numeric value, N can be specified as an :doc:`equal-style
variable <variable>`, which should be specified as v_name, where name
is the variable name. In this case, the variable is evaluated at the
beginning of a run to determine the **next** timestep at which the
string will be written out. On that timestep, the variable will be
evaluated again to determine the next timestep, etc.
Thus the variable should return timestep values. See the stagger()
and logfreq() and stride() math functions for :doc:`equal-style variables <variable>`, as examples of useful functions to use in
this context. For example, the following commands will print output at
evaluated again to determine the next timestep, etc. Thus the
variable should return timestep values. See the stagger() and
logfreq() and stride() math functions for :doc:`equal-style variables
<variable>`, as examples of useful functions to use in this
context. For example, the following commands will print output at
timesteps 10,20,30,100,200,300,1000,2000,etc:
.. code-block:: LAMMPS
@ -61,12 +77,12 @@ timesteps 10,20,30,100,200,300,1000,2000,etc:
The specified group-ID is ignored by this fix.
See the :doc:`variable <variable>` command for a description of *equal*
style variables which are the most useful ones to use with the fix
print command, since they are evaluated afresh each timestep that the
fix print line is output. Equal-style variables calculate formulas
involving mathematical operations, atom properties, group properties,
thermodynamic properties, global values calculated by a
See the :doc:`variable <variable>` command for a description of
*equal* style variables which are the most useful ones to use with the
fix print command, since they are evaluated afresh each timestep that
the fix print line is output. Equal-style variables calculate
formulas involving mathematical operations, atom properties, group
properties, thermodynamic properties, global values calculated by a
:doc:`compute <compute>` or :doc:`fix <fix>`, or references to other
:doc:`variables <variable>`.
@ -92,11 +108,13 @@ where ID is replaced with the fix-ID.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`. None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix. No global or per-atom quantities are stored by
this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""

14
doc/src/print.rst
View File

@ -46,6 +46,20 @@ lines of output, the string can be enclosed in triple quotes, as in
the last example above. If the text string contains variables, they
will be evaluated and their current values printed.
.. note::
As discussed on the :doc:`Commands parse <Commands_parse>` doc
page, the text string can use "immediate" variables, specified as
$(formula) with parenthesis, where the numeric formula has the same
syntax as equal-style variables described on the :doc:`variable
<variable>` doc page. This is a convenient way to evaluate a
formula immediately without using the variable command to define a
named variable and then use that variable in the text string. The
formula can include a trailing colon and format string which
determines the precision with which the numeric value is output.
This is also explained on the :doc:`Commands parse
<Commands_parse>` doc page.
If the *file* or *append* keyword is used, a filename is specified to
which the output will be written. If *file* is used, then the
filename is overwritten if it already exists. If *append* is used,

141
doc/src/variable.rst
View File

@ -11,7 +11,7 @@ Syntax
variable name style args ...
* name = name of variable to define
* style = *delete* or *index* or *loop* or *world* or *universe* or *uloop* or *string* or *format* or *getenv* or *file* or *atomfile* or *python* or *internal* or *equal* or *vector* or *atom*
* style = *delete* or *index* or *loop* or *world* or *universe* or *uloop* or *string* or *format* or *getenv* or *file* or *atomfile* or *python* or *timer* or *internal* or *equal* or *vector* or *atom*
.. parsed-literal::
@ -42,6 +42,7 @@ Syntax
*file* arg = filename
*atomfile* arg = filename
*python* arg = function
*timer* arg = no arguments
*internal* arg = numeric value
*equal* or *vector* or *atom* args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute/fix/variable references
numbers = 0.0, 100, -5.4, 2.8e-4, etc
@ -96,6 +97,13 @@ Examples
variable str format x %.6g
variable x delete
.. code-block:: LAMMPS
variable start timer
other commands
variable stop timer
print "Elapsed time: $(v_stop-v_start:%.6f)"
Description
"""""""""""
@ -108,32 +116,38 @@ part of a new input command. For variable styles that store multiple
strings, the :doc:`next <next>` command can be used to increment which
string is assigned to the variable. Variables of style *equal* store
a formula which when evaluated produces a single numeric value which
can be output either directly (see the :doc:`print <print>`, :doc:`fix print <fix_print>`, and :doc:`run every <run>` commands) or as part
of thermodynamic output (see the :doc:`thermo_style <thermo_style>`
command), or used as input to an averaging fix (see the :doc:`fix ave/time <fix_ave_time>` command). Variables of style *vector*
store a formula which produces a vector of such values which can be
used as input to various averaging fixes, or elements of which can be
part of thermodynamic output. Variables of style *atom* store a
formula which when evaluated produces one numeric value per atom which
can be output to a dump file (see the :doc:`dump custom <dump>` command)
or used as input to an averaging fix (see the :doc:`fix ave/chunk <fix_ave_chunk>` and :doc:`fix ave/atom <fix_ave_atom>`
commands). Variables of style *atomfile* can be used anywhere in an
input script that atom-style variables are used; they get their
per-atom values from a file rather than from a formula. Variables of
style *python* can be hooked to Python functions using code you
provide, so that the variable gets its value from the evaluation of
the Python code. Variables of style *internal* are used by a few
commands which set their value directly.
can be output either directly (see the :doc:`print <print>`, :doc:`fix
print <fix_print>`, and :doc:`run every <run>` commands) or as part of
thermodynamic output (see the :doc:`thermo_style <thermo_style>`
command), or used as input to an averaging fix (see the :doc:`fix
ave/time <fix_ave_time>` command). Variables of style *vector* store
a formula which produces a vector of such values which can be used as
input to various averaging fixes, or elements of which can be part of
thermodynamic output. Variables of style *atom* store a formula which
when evaluated produces one numeric value per atom which can be output
to a dump file (see the :doc:`dump custom <dump>` command) or used as
input to an averaging fix (see the :doc:`fix ave/chunk
<fix_ave_chunk>` and :doc:`fix ave/atom <fix_ave_atom>` commands).
Variables of style *atomfile* can be used anywhere in an input script
that atom-style variables are used; they get their per-atom values
from a file rather than from a formula. Variables of style *python*
can be hooked to Python functions using code you provide, so that the
variable gets its value from the evaluation of the Python code.
Variables of style *internal* are used by a few commands which set
their value directly.
.. note::
As discussed on the :doc:`Commands parse <Commands_parse>` doc
page, an input script can use "immediate" variables, specified as
$(formula) with parenthesis, where the formula has the same syntax as
equal-style variables described on this page. This is a convenient
way to evaluate a formula immediately without using the variable
command to define a named variable and then evaluate that
variable. See below for a more detailed discussion of this feature.
$(formula) with parenthesis, where the numeric formula has the same
syntax as equal-style variables described on this page. This is a
convenient way to evaluate a formula immediately without using the
variable command to define a named variable and then evaluate that
variable. The formula can include a trailing colon and format
string which determines the precision with which the numeric value
is generated. This is also explained on the :doc:`Commands parse
<Commands_parse>` doc page.
In the discussion that follows, the "name" of the variable is the
arbitrary string that is the first argument in the variable command.
@ -160,22 +174,19 @@ simulation.
Variables of style *equal* and *vector* and *atom* can be used as
inputs to various other commands which evaluate their formulas as
needed, e.g. at different timesteps during a :doc:`run <run>`.
needed, e.g. at different timesteps during a :doc:`run <run>`. In
this context, variables of style *timer* or *internal* or *python* can
be used in place of an equal-style variable, with the following two
caveats.
Variables of style *internal* can be used in place of an equal-style
variable, except by commands that set the value stored by the
internal-style variable. Thus any command that states it can use an
equal-style variable as an argument, can also use an internal-style
variable. This means that when the command evaluates the variable, it
will use the value set (internally) by another command.
Variables of style *python* can be used in place of an equal-style
variable so long as the associated Python function, as defined by the
:doc:`python <python>` command, returns a numeric value. Thus any
command that states it can use an equal-style variable as an argument,
can also use such a python-style variable. This means that when the
LAMMPS command evaluates the variable, the Python function will be
executed.
First, internal-style variables can be used except by commands that
set the value stored by the internal variable. When the LAMMPS
command evaluates the internal-style variable, it will use the value
set (internally) by another command. Second, python-style variables
can be used so long as the associated Python function, as defined by
the :doc:`python <python>` command, returns a numeric value. When the
LAMMPS command evaluates the python-style variable, the Python
function will be executed.
.. note::
@ -271,14 +282,15 @@ N1 <= N2 and N2 >= 0 is required.
For the *world* style, one or more strings are specified. There must
be one string for each processor partition or "world". LAMMPS can be
run with multiple partitions via the :doc:`-partition command-line switch <Run_options>`. This variable command assigns one string to
run with multiple partitions via the :doc:`-partition command-line
switch <Run_options>`. This variable command assigns one string to
each world. All processors in the world are assigned the same string.
The next command cannot be used with *equal* style variables, since
there is only one value per world. This style of variable is useful
when you wish to run different simulations on different partitions, or
when performing a parallel tempering simulation (see the
:doc:`temper <temper>` command), to assign different temperatures to
different partitions.
when performing a parallel tempering simulation (see the :doc:`temper
<temper>` command), to assign different temperatures to different
partitions.
For the *universe* style, one or more strings are specified. There
must be at least as many strings as there are processor partitions or
@ -313,11 +325,12 @@ 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.
If you simply wish to print a variable value with desired precision to
the screen or logfile via the :doc:`print <print>` or :doc:`fix print <fix_print>` commands, you can also do this by specifying an
"immediate" variable with a trailing colon and format string, as part
of the string argument of those commands. This is explained on the
:doc:`Commands parse <Commands_parse>` doc page.
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
:doc:`fix print <fix_print>` commands, you can also do this by
specifying an "immediate" variable with a trailing colon and format
string, as part of the string argument of those commands. This is
explained on the :doc:`Commands parse <Commands_parse>` doc page.
For the *getenv* style, a single string is assigned to the variable
which should be the name of an environment variable. When the
@ -412,14 +425,25 @@ python-style variable can be used in place of an equal-style variable
anywhere in an input script, e.g. as an argument to another command
that allows for equal-style variables.
For the *timer* style no additional argument is specified. The value of
the variable is set by querying the current elapsed wall time of the
simulation. This is done at the point in time when the variable is
defined in the input script. If a second timer-style variable is also
defined, then a simple formula can be used to calculate the elapsed time
between the two timers, as in the example at the top of this manual
entry. As mentioned above, timer-style variables can be redefined
elsewhere in the input script, so the same pair of variables can be used
in a loop or to time a series of operations.
For the *internal* style a numeric value is provided. This value will
be assigned to the variable until a LAMMPS command sets it to a new
value. There are currently only two LAMMPS commands that require
*internal* variables as inputs, because they reset them:
:doc:`create_atoms <create_atoms>` and :doc:`fix controller <fix_controller>`. As mentioned above, an
internal-style variable can be used in place of an equal-style
variable anywhere else in an input script, e.g. as an argument to
another command that allows for equal-style variables.
:doc:`create_atoms <create_atoms>` and :doc:`fix controller
<fix_controller>`. As mentioned above, an internal-style variable can
be used in place of an equal-style variable anywhere else in an input
script, e.g. as an argument to another command that allows for
equal-style variables.
----------
@ -1383,14 +1407,15 @@ commands:
The first run is performed using one setting for the pairwise
potential defined by the :doc:`pair_style <pair_style>` and
:doc:`pair_coeff <pair_coeff>` commands. The potential energy is
evaluated on the final timestep and stored by the :doc:`compute pe <compute_pe>` compute (this is done by the
:doc:`thermo_style <thermo_style>` command). Then a pair coefficient is
changed, altering the potential energy of the system. When the
potential energy is printed via the "e" variable, LAMMPS will use the
potential energy value stored by the :doc:`compute pe <compute_pe>`
compute, thinking it is current. There are many other commands which
could alter the state of the system between runs, causing a variable
to evaluate incorrectly.
evaluated on the final timestep and stored by the :doc:`compute pe
<compute_pe>` compute (this is done by the :doc:`thermo_style
<thermo_style>` command). Then a pair coefficient is changed,
altering the potential energy of the system. When the potential
energy is printed via the "e" variable, LAMMPS will use the potential
energy value stored by the :doc:`compute pe <compute_pe>` compute,
thinking it is current. There are many other commands which could
alter the state of the system between runs, causing a variable to
evaluate incorrectly.
The solution to this issue is the same as for case (2) above, namely
perform a 0-timestep run before the variable is evaluated to insure

38
src/.gitignore vendored
View File

@ -37,6 +37,12 @@
/library_mdi.h
/mdi_engine.cpp
/mdi_engine.h
/fix_mdi_aimd.cpp
/fix_mdi_aimd.h
/mdi_command.cpp
/mdi_command.h
/mdi_plugin.cpp
/mdi_plugin.h
/fix_brownian*.cpp
/fix_brownian*.h
@ -269,6 +275,32 @@
/pair_bpm_spring.cpp
/pair_bpm_spring.h
/boundary_correction.cpp
/boundary_correction.h
/electrode_accel_interface.h
/electrode_kspace.h
/electrode_math.h
/electrode_matrix.cpp
/electrode_matrix.h
/electrode_vector.cpp
/electrode_vector.h
/ewald_electrode.cpp
/ewald_electrode.h
/fix_electrode_conp.cpp
/fix_electrode_conp.h
/fix_electrode_conq.cpp
/fix_electrode_conq.h
/fix_electrode_thermo.cpp
/fix_electrode_thermo.h
/pppm_electrode.cpp
/pppm_electrode.h
/slab_2d.cpp
/slab_2d.h
/slab_dipole.cpp
/slab_dipole.h
/wire_dipole.cpp
/wire_dipole.h
/compute_adf.cpp
/compute_adf.h
/compute_contact_atom.cpp
@ -579,6 +611,7 @@
/dihedral_table_cut.h
/dump_atom_adios.cpp
/dump_atom_adios.h
/adios_common.h
/dump_atom_gz.cpp
/dump_atom_gz.h
/dump_atom_zstd.cpp
@ -672,6 +705,8 @@
/fix_client_md.h
/fix_cmap.cpp
/fix_cmap.h
/fix_damping_cundall.cpp
/fix_damping_cundall.h
/fix_dpd_energy.cpp
/fix_dpd_energy.h
/fix_electron_stopping.cpp
@ -728,6 +763,7 @@
/fix_langevin_eff.h
/fix_latte.cpp
/fix_latte.h
/latboltz_const.h
/fix_lb_fluid.cpp
/fix_lb_fluid.h
/fix_lb_momentum.cpp
@ -889,6 +925,8 @@
/fix_ti_spring.h
/fix_tune_kspace.cpp
/fix_tune_kspace.h
/fix_viscous_sphere.cpp
/fix_viscous_sphere.h
/fix_wall_body_polygon.cpp
/fix_wall_body_polygon.h
/fix_wall_body_polyhedron.cpp

15
src/ELECTRODE/fix_electrode_thermo.cpp
View File

@ -33,6 +33,8 @@ using namespace LAMMPS_NS;
enum { CONST, EQUAL };
/* ----------------------------------------------------------------------- */
// 0 1 2 3 4
// fix fxupdate group1 electrode/thermo pot1 eta couple group2 pot2
FixElectrodeThermo::FixElectrodeThermo(LAMMPS *lmp, int narg, char **arg) :
@ -49,6 +51,15 @@ FixElectrodeThermo::FixElectrodeThermo(LAMMPS *lmp, int narg, char **arg) :
if (group_psi_var_styles[0] == CONST) delta_psi_0 = group_psi[1] - group_psi[0];
}
/* ----------------------------------------------------------------------- */
FixElectrodeThermo::~FixElectrodeThermo()
{
delete thermo_random;
}
/* ----------------------------------------------------------------------- */
void FixElectrodeThermo::compute_macro_matrices()
{
FixElectrodeConp::compute_macro_matrices();
@ -57,6 +68,8 @@ void FixElectrodeThermo::compute_macro_matrices()
(macro_capacitance[0][0] + macro_capacitance[1][1] + 2 * macro_capacitance[0][1]);
}
/* ----------------------------------------------------------------------- */
void FixElectrodeThermo::pre_update()
{
// total electrode charges after last step, required for update psi
@ -72,6 +85,8 @@ void FixElectrodeThermo::pre_update()
MPI_Allreduce(MPI_IN_PLACE, &group_q_old, NUM_GROUPS, MPI_DOUBLE, MPI_SUM, world);
}
/* ----------------------------------------------------------------------- */
void FixElectrodeThermo::update_psi()
{
double const dt = update->dt;

1
src/ELECTRODE/fix_electrode_thermo.h
View File

@ -33,6 +33,7 @@ namespace LAMMPS_NS {
class FixElectrodeThermo : public FixElectrodeConp {
public:
FixElectrodeThermo(class LAMMPS *, int, char **);
~FixElectrodeThermo() override;
void update_psi() override;
void pre_update() override;

21
src/ELECTRODE/pppm_electrode.h
View File

@ -27,27 +27,6 @@ KSpaceStyle(pppm/electrode, PPPMElectrode);
#include "electrode_kspace.h"
#include "pppm.h"
#if defined(FFT_FFTW3)
#define LMP_FFT_LIB "FFTW3"
#elif defined(FFT_MKL)
#define LMP_FFT_LIB "MKL FFT"
#elif defined(FFT_CUFFT)
#define LMP_FFT_LIB "cuFFT"
#else
#define LMP_FFT_LIB "KISS FFT"
#endif
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define LMP_FFT_PREC "single"
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define LMP_FFT_PREC "double"
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
namespace LAMMPS_NS {
class PPPMElectrode : public PPPM, public ElectrodeKSpace {

2
src/KOKKOS/comm_kokkos.h
View File

@ -37,6 +37,8 @@ class CommKokkos : public CommBrick {
~CommKokkos() override;
void init() override;
using CommBrick::forward_comm;
using CommBrick::reverse_comm;
void forward_comm(int dummy = 0) override; // forward comm of atom coords
void reverse_comm() override; // reverse comm of atom coords
void exchange() override; // move atoms to new procs

2
src/KOKKOS/comm_tiled_kokkos.h
View File

@ -26,6 +26,8 @@ class CommTiledKokkos : public CommTiled {
CommTiledKokkos(class LAMMPS *, class Comm *);
~CommTiledKokkos() override;
using CommTiled::forward_comm;
using CommTiled::reverse_comm;
void forward_comm(int dummy = 0) override; // forward comm of atom coords
void reverse_comm() override; // reverse comm of forces
void exchange() override; // move atoms to new procs

195
src/KOKKOS/compute_orientorder_atom_kokkos.cpp
View File

@ -13,7 +13,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Stan Moore (SNL)
Contributing authors: Stan Moore (SNL)
Tomas Oppelstrup (LLNL): Optimization which reduces the number
of iterations in the L,m1,m2 loops (by a factor of up to 10), and
avoids evaluation of Ylm functions of negative m
------------------------------------------------------------------------- */
#include "compute_orientorder_atom_kokkos.h"
@ -34,16 +38,15 @@
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathSpecial;
using namespace std;
using MathSpecial::factorial;
#ifdef DBL_EPSILON
#define MY_EPSILON (10.0*DBL_EPSILON)
static constexpr double MY_EPSILON = (10.0 * DBL_EPSILON);
#else
#define MY_EPSILON (10.0*2.220446049250313e-16)
static constexpr double MY_EPSILON = (10.0 * 2.220446049250313e-16);
#endif
#define QEPSILON 1.0e-6
static constexpr double QEPSILON = 1.0e-6;
/* ---------------------------------------------------------------------- */
@ -58,6 +61,20 @@ ComputeOrientOrderAtomKokkos<DeviceType>::ComputeOrientOrderAtomKokkos(LAMMPS *l
datamask_modify = EMPTY_MASK;
host_flag = (execution_space == Host);
d_qnormfac = t_sna_1d("orientorder/atom:qnormfac",nqlist);
d_qnormfac2 = t_sna_1d("orientorder/atom:qnormfac2",nqlist);
auto h_qnormfac = Kokkos::create_mirror_view(d_qnormfac);
auto h_qnormfac2 = Kokkos::create_mirror_view(d_qnormfac2);
for (int il = 0; il < nqlist; il++) {
h_qnormfac[il] = qnormfac[il];
h_qnormfac2[il] = qnormfac2[il];
}
Kokkos::deep_copy(d_qnormfac,h_qnormfac);
Kokkos::deep_copy(d_qnormfac2,h_qnormfac2);
}
/* ---------------------------------------------------------------------- */
@ -140,7 +157,7 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::compute_peratom()
chunk_offset = 0;
if (chunk_size > (int)d_ncount.extent(0)) {
d_qnm = t_sna_3c("orientorder/atom:qnm",chunk_size,nqlist,2*qmax+1);
d_qnm = t_sna_3c("orientorder/atom:qnm",chunk_size,nqlist,qmax+1);
d_ncount = t_sna_1i("orientorder/atom:ncount",chunk_size);
}
@ -332,7 +349,6 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::operator() (TagComputeOrientOrder
calc_boop2(ncount, ii);
}
/* ----------------------------------------------------------------------
select3 routine from Numerical Recipes (slightly modified)
find k smallest values in array of length n
@ -466,26 +482,14 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop1(int /*ncount*/, int ii
//d_qnm(ii,il,l).re += polar_prefactor(l, 0, costheta);
const double polar_pf = polar_prefactor(l, 0, costheta);
Kokkos::atomic_add(&(d_qnm(ii,il,l).re), polar_pf);
Kokkos::atomic_add(&(d_qnm(ii,il,0).re), polar_pf);
SNAcomplex expphim = {expphi.re,expphi.im};
for (int m = 1; m <= +l; m++) {
const double prefactor = polar_prefactor(l, m, costheta);
SNAcomplex ylm = {prefactor * expphim.re, prefactor * expphim.im};
//d_qnm(ii,il,m+l).re += ylm.re;
//d_qnm(ii,il,m+l).im += ylm.im;
Kokkos::atomic_add(&(d_qnm(ii,il,m+l).re), ylm.re);
Kokkos::atomic_add(&(d_qnm(ii,il,m+l).im), ylm.im);
if (m & 1) {
//d_qnm(ii,il,-m+l).re -= ylm.re;
//d_qnm(ii,il,-m+l).im += ylm.im;
Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), -ylm.re);
Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), ylm.im);
} else {
//d_qnm(ii,il,-m+l).re += ylm.re;
//d_qnm(ii,il,-m+l).im -= ylm.im;
Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), ylm.re);
Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), -ylm.im);
}
Kokkos::atomic_add(&(d_qnm(ii,il,m).re), ylm.re);
Kokkos::atomic_add(&(d_qnm(ii,il,m).im), ylm.im);
// Skip calculation of qnm for m<0 due to symmetry
SNAcomplex tmp;
tmp.re = expphim.re*expphi.re - expphim.im*expphi.im;
tmp.im = expphim.re*expphi.im + expphim.im*expphi.re;
@ -510,7 +514,7 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
double facn = 1.0 / ncount;
for (int il = 0; il < nqlist; il++) {
int l = d_qlist[il];
for (int m = 0; m < 2*l+1; m++) {
for (int m = 0; m < l+1; m++) {
d_qnm(ii,il,m).re *= facn;
d_qnm(ii,il,m).im *= facn;
}
@ -522,57 +526,57 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
int jj = 0;
for (int il = 0; il < nqlist; il++) {
int l = d_qlist[il];
double qnormfac = sqrt(MY_4PI/(2*l+1));
double qm_sum = 0.0;
for (int m = 0; m < 2*l+1; m++)
qm_sum += d_qnm(ii,il,m).re*d_qnm(ii,il,m).re + d_qnm(ii,il,m).im*d_qnm(ii,il,m).im;
d_qnarray(i,jj++) = qnormfac * sqrt(qm_sum);
double qm_sum = d_qnm(ii,il,0).re*d_qnm(ii,il,0).re;
for (int m = 1; m < l+1; m++)
qm_sum += 2.0*(d_qnm(ii,il,m).re*d_qnm(ii,il,m).re + d_qnm(ii,il,m).im*d_qnm(ii,il,m).im);
d_qnarray(i,jj++) = d_qnormfac(il) * sqrt(qm_sum);
}
// calculate W_l
if (wlflag) {
int idxcg_count = 0;
int nterms = 0;
int widx_count = 0;
if (wlflag || wlhatflag) {
for (int il = 0; il < nqlist; il++) {
int l = d_qlist[il];
double wlsum = 0.0;
for (int m1 = 0; m1 < 2*l+1; m1++) {
for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
int m = m1 + m2 - l;
SNAcomplex qm1qm2;
qm1qm2.re = d_qnm(ii,il,m1).re*d_qnm(ii,il,m2).re - d_qnm(ii,il,m1).im*d_qnm(ii,il,m2).im;
qm1qm2.im = d_qnm(ii,il,m1).re*d_qnm(ii,il,m2).im + d_qnm(ii,il,m1).im*d_qnm(ii,il,m2).re;
wlsum += (qm1qm2.re*d_qnm(ii,il,m).re + qm1qm2.im*d_qnm(ii,il,m).im)*d_cglist[idxcg_count];
idxcg_count++;
for (int m1 = -l; m1 <= 0; m1++) {
const int sgn = 1 - 2*(m1&1); // sgn = (-1)^m1
for (int m2 = 0; m2 <= ((-m1)>>1); m2++) {
const int m3 = -(m1 + m2);
// Loop enforces -L<=m1<=0<=m2<=m3<=L, and m1+m2+m3=0
// For even L, W3j is invariant under permutation of
// (m1,m2,m3) and (m1,m2,m3)->(-m1,-m2,-m3). The loop
// structure enforces visiting only one member of each
// such symmetry (invariance) group.
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m])
SNAcomplex Q1Q2;
Q1Q2.re = (d_qnm(ii,il,-m1).re*d_qnm(ii,il,m2).re + d_qnm(ii,il,-m1).im*d_qnm(ii,il,m2).im)*sgn;
Q1Q2.im = (d_qnm(ii,il,-m1).re*d_qnm(ii,il,m2).im - d_qnm(ii,il,-m1).im*d_qnm(ii,il,m2).re)*sgn;
const double Q1Q2Q3 = Q1Q2.re*d_qnm(ii,il,m3).re - Q1Q2.im*d_qnm(ii,il,m3).im;
const double c = d_w3jlist[widx_count++];
wlsum += Q1Q2Q3*c;
}
}
d_qnarray(i,jj++) = wlsum/sqrt(2.0*l+1.0);
d_qnarray(i,jj++) = wlsum/d_qnormfac2(il);
nterms++;
}
}
// calculate W_l_hat
if (wlhatflag) {
int idxcg_count = 0;
const int jptr = jj-nterms;
if (!wlflag) jj = jptr;
for (int il = 0; il < nqlist; il++) {
int l = d_qlist[il];
double wlsum = 0.0;
for (int m1 = 0; m1 < 2*l+1; m1++) {
for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
const int m = m1 + m2 - l;
SNAcomplex qm1qm2;
qm1qm2.re = d_qnm(ii,il,m1).re*d_qnm(ii,il,m2).re - d_qnm(ii,il,m1).im*d_qnm(ii,il,m2).im;
qm1qm2.im = d_qnm(ii,il,m1).re*d_qnm(ii,il,m2).im + d_qnm(ii,il,m1).im*d_qnm(ii,il,m2).re;
wlsum += (qm1qm2.re*d_qnm(ii,il,m).re + qm1qm2.im*d_qnm(ii,il,m).im)*d_cglist[idxcg_count];
idxcg_count++;
}
}
if (d_qnarray(i,il) < QEPSILON)
d_qnarray(i,jj++) = 0.0;
else {
const double qnormfac = sqrt(MY_4PI/(2*l+1));
const double qnfac = qnormfac/d_qnarray(i,il);
d_qnarray(i,jj++) = wlsum/sqrt(2.0*l+1.0)*(qnfac*qnfac*qnfac);
const double qnfac = d_qnormfac(il)/d_qnarray(i,il);
d_qnarray(i,jj++) = d_qnarray(i,jptr+il) * (qnfac*qnfac*qnfac) * d_qnormfac2(il);
}
}
}
@ -588,9 +592,15 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
d_qnarray(i,jj++) = 0.0;
}
else {
const double qnormfac = sqrt(MY_4PI/(2*l+1));
const double qnfac = qnormfac/d_qnarray(i,il);
for (int m = 0; m < 2*l+1; m++) {
const double qnfac = d_qnormfac(il)/d_qnarray(i,il);
for (int m = -l; m < 0; m++) {
// Computed only qnm for m>=0.
// qnm[-m] = (-1)^m * conjg(qnm[m])
const int sgn = 1 - 2*(m&1); // sgn = (-1)^m
d_qnarray(i,jj++) = d_qnm(ii,il,-m).re * qnfac * sgn;
d_qnarray(i,jj++) = -d_qnm(ii,il,-m).im * qnfac * sgn;
}
for (int m = 0; m < l+1; m++) {
d_qnarray(i,jj++) = d_qnm(ii,il,m).re * qnfac;
d_qnarray(i,jj++) = d_qnm(ii,il,m).im * qnfac;
}
@ -652,70 +662,21 @@ double ComputeOrientOrderAtomKokkos<DeviceType>::associated_legendre(int l, int
}
/* ----------------------------------------------------------------------
assign Clebsch-Gordan coefficients
using the quasi-binomial formula VMK 8.2.1(3)
specialized for case j1=j2=j=l
Initialize table of Wigner 3j symbols
------------------------------------------------------------------------- */
template<class DeviceType>
void ComputeOrientOrderAtomKokkos<DeviceType>::init_clebsch_gordan()
void ComputeOrientOrderAtomKokkos<DeviceType>::init_wigner3j()
{
double sum,dcg,sfaccg, sfac1, sfac2;
int m, aa2, bb2, cc2;
int ifac, idxcg_count;
ComputeOrientOrderAtom::init_wigner3j();
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m1 = 0; m1 < 2*l+1; m1++)
for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++)
idxcg_count++;
}
idxcg_max = idxcg_count;
d_cglist = t_sna_1d("orientorder/atom:d_cglist",idxcg_max);
auto h_cglist = Kokkos::create_mirror_view(d_cglist);
d_w3jlist = t_sna_1d("computeorientorderatom:w3jlist",widx_max);
auto h_w3jlist = Kokkos::create_mirror_view(d_w3jlist);
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m1 = 0; m1 < 2*l+1; m1++) {
aa2 = m1 - l;
for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
bb2 = m2 - l;
m = aa2 + bb2 + l;
for (int i = 0; i< widx_max; i++)
h_w3jlist(i) = w3jlist[i];
sum = 0.0;
for (int z = MAX(0, MAX(-aa2, bb2));
z <= MIN(l, MIN(l - aa2, l + bb2)); z++) {
ifac = z % 2 ? -1 : 1;
sum += ifac /
(factorial(z) *
factorial(l - z) *
factorial(l - aa2 - z) *
factorial(l + bb2 - z) *
factorial(aa2 + z) *
factorial(-bb2 + z));
}
cc2 = m - l;
sfaccg = sqrt(factorial(l + aa2) *
factorial(l - aa2) *
factorial(l + bb2) *
factorial(l - bb2) *
factorial(l + cc2) *
factorial(l - cc2) *
(2*l + 1));
sfac1 = factorial(3*l + 1);
sfac2 = factorial(l);
dcg = sqrt(sfac2*sfac2*sfac2 / sfac1);
h_cglist[idxcg_count] = sum * dcg * sfaccg;
idxcg_count++;
}
}
}
Kokkos::deep_copy(d_cglist,h_cglist);
Kokkos::deep_copy(d_w3jlist,h_w3jlist);
}
/* ----------------------------------------------------------------------

5
src/KOKKOS/compute_orientorder_atom_kokkos.h
View File

@ -71,6 +71,7 @@ class ComputeOrientOrderAtomKokkos : public ComputeOrientOrderAtom {
void init() override;
void compute_peratom() override;
t_sna_1i d_qlist;
t_sna_1d d_qnormfac,d_qnormfac2;
template<class TagStyle>
void check_team_size_for(int, int&, int);
@ -127,8 +128,8 @@ class ComputeOrientOrderAtomKokkos : public ComputeOrientOrderAtom {
KOKKOS_INLINE_FUNCTION
double associated_legendre(int, int, double) const;
void init_clebsch_gordan() override;
t_sna_1d d_cglist; // Clebsch-Gordan coeffs
void init_wigner3j();
t_sna_1d d_w3jlist; // Wigner coeffs
};
}

15
src/KOKKOS/fftdata_kokkos.h
View File

@ -23,20 +23,7 @@
#ifndef LMP_FFT_DATA_KOKKOS_H
#define LMP_FFT_DATA_KOKKOS_H
// User-settable FFT precision
// FFT_PRECISION = 1 is single-precision complex (4-byte real, 4-byte imag)
// FFT_PRECISION = 2 is double-precision complex (8-byte real, 8-byte imag)
#ifdef FFT_SINGLE
#define FFT_PRECISION 1
#define MPI_FFT_SCALAR MPI_FLOAT
typedef float FFT_SCALAR;
#else
#define FFT_PRECISION 2
#define MPI_FFT_SCALAR MPI_DOUBLE
typedef double FFT_SCALAR;
#endif
#include "lmpfftsettings.h"
// -------------------------------------------------------------------------

6
src/KOKKOS/npair_kokkos.cpp
View File

@ -31,6 +31,8 @@ namespace LAMMPS_NS {
template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI, int SIZE>
NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI,SIZE>::NPairKokkos(LAMMPS *lmp) : NPair(lmp) {
last_stencil_old = -1;
// use 1D view for scalars to reduce GPU memory operations
d_scalars = typename AT::t_int_1d("neighbor:scalars",2);
@ -114,9 +116,11 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI,SIZE>::copy_stencil_info()
NPair::copy_stencil_info();
nstencil = ns->nstencil;
if (ns->last_stencil == update->ntimestep) {
if (ns->last_stencil != last_stencil_old) {
// copy stencil to device as it may have changed
last_stencil_old = ns->last_stencil;
int maxstencil = ns->get_maxstencil();
if (maxstencil > (int)k_stencil.extent(0))

2
src/KOKKOS/npair_kokkos.h
View File

@ -137,7 +137,7 @@ class NPairKokkos : public NPair {
// data from NStencil class
int nstencil;
int nstencil,last_stencil_old;
DAT::tdual_int_1d k_stencil; // # of J neighs for each I
DAT::tdual_int_1d_3 k_stencilxyz;
};

7
src/KOKKOS/rand_pool_wrap_kokkos.h
View File

@ -48,7 +48,6 @@ class RandPoolWrap : protected Pointers {
void destroy();
void init(RanMars*, int);
KOKKOS_INLINE_FUNCTION
RandWrap get_state() const
{
#ifdef LMP_KOKKOS_GPU
@ -68,11 +67,7 @@ class RandPoolWrap : protected Pointers {
return rand_wrap;
}
KOKKOS_INLINE_FUNCTION
void free_state(RandWrap) const
{
}
void free_state(RandWrap) const {}
private:
class RanMars **random_thr;

18
src/KSPACE/fft3d.h
View File

@ -13,24 +13,8 @@
#include <mpi.h>
// User-settable FFT precision
// FFT_PRECISION = 1 is single-precision complex (4-byte real, 4-byte imag)
// FFT_PRECISION = 2 is double-precision complex (8-byte real, 8-byte imag)
#ifdef FFT_SINGLE
#define FFT_PRECISION 1
typedef float FFT_SCALAR;
#else
#define FFT_PRECISION 2
typedef double FFT_SCALAR;
#endif
// if user set FFTW, it means FFTW3
#ifdef FFT_FFTW
#define FFT_FFTW3
#endif
#include "lmpfftsettings.h"
// -------------------------------------------------------------------------

24
src/KSPACE/pppm.h
View File

@ -21,29 +21,7 @@ KSpaceStyle(pppm,PPPM);
#define LMP_PPPM_H
#include "kspace.h"
#if defined(FFT_FFTW3)
#define LMP_FFT_LIB "FFTW3"
#elif defined(FFT_MKL)
#define LMP_FFT_LIB "MKL FFT"
#elif defined(FFT_CUFFT)
#define LMP_FFT_LIB "cuFFT"
#elif defined(FFT_HIPFFT)
#define LMP_FFT_LIB "hipFFT"
#else
#define LMP_FFT_LIB "KISS FFT"
#endif
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define LMP_FFT_PREC "single"
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define LMP_FFT_PREC "double"
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#include "lmpfftsettings.h"
namespace LAMMPS_NS {

19
src/KSPACE/pppm_disp.h
View File

@ -21,24 +21,7 @@ KSpaceStyle(pppm/disp,PPPMDisp);
#define LMP_PPPM_DISP_H
#include "kspace.h"
#if defined(FFT_FFTW3)
#define LMP_FFT_LIB "FFTW3"
#elif defined(FFT_MKL)
#define LMP_FFT_LIB "MKL FFT"
#else
#define LMP_FFT_LIB "KISS FFT"
#endif
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define LMP_FFT_PREC "single"
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define LMP_FFT_PREC "double"
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#include "lmpfftsettings.h"
namespace LAMMPS_NS {

8
src/KSPACE/remap.h
View File

@ -13,13 +13,7 @@
#include <mpi.h>
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#include "lmpfftsettings.h"
// details of how to do a 3d remap

20
src/Makefile
View File

@ -71,7 +71,6 @@ PACKAGE = \
dpd-smooth \
drude \
eff \
electrode \
extra-compute \
extra-dump \
extra-fix \
@ -137,6 +136,7 @@ PACKAGE = \
yaff \
atc \
dielectric \
electrode \
ml-iap \
phonon
@ -163,7 +163,6 @@ PACKMOST = \
dpd-react \
dpd-smooth \
drude \
electrode \
eff \
extra-compute \
extra-dump \
@ -232,7 +231,8 @@ PACKLIB = \
ml-quip \
scafacos \
machdyn \
vtk
vtk \
electrode
PACKSYS = compress latboltz mpiio python
@ -281,14 +281,15 @@ help:
@echo 'make package list available packages and their dependencies'
@echo 'make package-status (ps) status of all packages'
@echo 'make package-installed (pi) list of installed packages'
@echo 'make yes-package install a single pgk in src dir'
@echo 'make yes-package install a single pkg in src dir'
@echo 'make no-package remove a single pkg from src dir'
@echo 'make yes-all install all pgks in src dir'
@echo 'make yes-all install all pkgs in src dir'
@echo 'make no-all remove all pkgs from src dir'
@echo 'make yes-basic install a few commonly used pgks'
@echo 'make no-installed (ni) quick uninstall of installed pkgs'
@echo 'make yes-basic install a few commonly used pkgs'
@echo 'make no-basic remove a few commonly used pkgs'
@echo 'make yes-most install most pgks w/o libs'
@echo 'make no-most remove most pgks w/o libs'
@echo 'make yes-most install most pkgs w/o libs'
@echo 'make no-most remove most pkgs w/o libs'
@echo 'make yes-lib install all pkgs with libs (included or ext)'
@echo 'make no-lib remove all pkgs with libs (included or ext)'
@echo 'make yes-ext install all pkgs with external libs'
@ -547,6 +548,9 @@ yes-all:
no-all:
@for p in $(PACKAGE); do $(MAKE) no-$$p; done
no-installed ni:
@for p in $(PACKAGESORTED); do $(SHELL) Package.sh $$p info && $(MAKE) no-$$p || : ; done
yes-standard yes-std:
@echo 'There are no more "standard" or "user" packages in LAMMPS'

10
src/PHONON/fix_phonon.h
View File

@ -32,15 +32,9 @@ FixStyle(phonon,FixPhonon);
#ifndef FIX_PHONON_H
#define FIX_PHONON_H
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#include "fix.h"
#include "lmpfftsettings.h"
#include <complex>
#include <map>

9
src/Package.sh
View File

@ -41,6 +41,15 @@ elif (test $2 = "installed") then
echo "Installed YES: package $1"
fi
# info, exit with true/false status depending on whether a package is installed
elif (test $2 = "info") then
if (test $installed = 1) then
exit 0
else
exit 1
fi
# update, only if installed
# perform a re-install, but only if the package is already installed

257
src/compute_orientorder_atom.cpp
View File

@ -12,8 +12,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
Axel Kohlmeyer (Temple U)
Contributing authors: Aidan Thompson (SNL), Axel Kohlmeyer (Temple U)
Tomas Oppelstrup (LLNL): Optimization which reduces the number
of iterations in the L,m1,m2 loops (by a factor of up to 10), and
avoids evaluation of Ylm functions of negative m
------------------------------------------------------------------------- */
#include "compute_orientorder_atom.h"
@ -36,21 +39,22 @@
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathSpecial;
using MathSpecial::factorial;
#ifdef DBL_EPSILON
#define MY_EPSILON (10.0 * DBL_EPSILON)
static constexpr double MY_EPSILON = (10.0 * DBL_EPSILON);
#else
#define MY_EPSILON (10.0 * 2.220446049250313e-16)
static constexpr double MY_EPSILON = (10.0 * 2.220446049250313e-16);
#endif
#define QEPSILON 1.0e-6
static constexpr double QEPSILON = 1.0e-6;
/* ---------------------------------------------------------------------- */
ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), qlist(nullptr), distsq(nullptr), nearest(nullptr), rlist(nullptr),
qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr), cglist(nullptr)
Compute(lmp, narg, arg), qlist(nullptr), qnormfac(nullptr), qnormfac2(nullptr), distsq(nullptr),
nearest(nullptr), rlist(nullptr), qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr),
w3jlist(nullptr)
{
if (narg < 3) error->all(FLERR, "Illegal compute orientorder/atom command");
@ -147,9 +151,17 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
nmax = 0;
maxneigh = 0;
memory->create(qnormfac, nqlist, "orientorder/atom:qnormfac");
memory->create(qnormfac2, nqlist, "orientorder/atom:qnormfac2");
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
qnormfac[il] = sqrt(MY_4PI / (2 * l + 1));
qnormfac2[il] = sqrt(2 * l + 1);
}
}
/* ---------------------------------------------------------------------- */
/* --------------------------------------------------------------------- */
ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
{
@ -160,9 +172,11 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
memory->destroy(rlist);
memory->destroy(nearest);
memory->destroy(qlist);
memory->destroy(qnormfac);
memory->destroy(qnormfac2);
memory->destroy(qnm_r);
memory->destroy(qnm_i);
memory->destroy(cglist);
memory->destroy(w3jlist);
}
/* ---------------------------------------------------------------------- */
@ -178,8 +192,8 @@ void ComputeOrientOrderAtom::init()
memory->destroy(qnm_r);
memory->destroy(qnm_i);
memory->create(qnm_r, nqlist, 2 * qmax + 1, "orientorder/atom:qnm_r");
memory->create(qnm_i, nqlist, 2 * qmax + 1, "orientorder/atom:qnm_i");
memory->create(qnm_r, nqlist, qmax + 1, "orientorder/atom:qnm_r");
memory->create(qnm_i, nqlist, qmax + 1, "orientorder/atom:qnm_i");
// need an occasional full neighbor list
@ -188,7 +202,7 @@ void ComputeOrientOrderAtom::init()
if ((modify->get_compute_by_style("orientorder/atom").size() > 1) && (comm->me == 0))
error->warning(FLERR, "More than one instance of compute orientorder/atom");
if (wlflag || wlhatflag) init_clebsch_gordan();
if (wlflag || wlhatflag) init_wigner3j();
}
/* ---------------------------------------------------------------------- */
@ -427,7 +441,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m = 0; m < 2 * l + 1; m++) {
for (int m = 0; m < l + 1; m++) {
qnm_r[il][m] = 0.0;
qnm_i[il][m] = 0.0;
}
@ -458,7 +472,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
// Ylm, -l <= m <= l
// sign convention: sign(Yll(0,0)) = (-1)^l
qnm_r[il][l] += polar_prefactor(l, 0, costheta);
qnm_r[il][0] += polar_prefactor(l, 0, costheta);
double expphim_r = expphi_r;
double expphim_i = expphi_i;
for (int m = 1; m <= +l; m++) {
@ -466,15 +480,9 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
double prefactor = polar_prefactor(l, m, costheta);
double ylm_r = prefactor * expphim_r;
double ylm_i = prefactor * expphim_i;
qnm_r[il][m + l] += ylm_r;
qnm_i[il][m + l] += ylm_i;
if (m & 1) {
qnm_r[il][-m + l] -= ylm_r;
qnm_i[il][-m + l] += ylm_i;
} else {
qnm_r[il][-m + l] += ylm_r;
qnm_i[il][-m + l] -= ylm_i;
}
qnm_r[il][m] += ylm_r;
qnm_i[il][m] += ylm_i;
// Skip calculation of qnm for m<0 due to symmetry
double tmp_r = expphim_r * expphi_r - expphim_i * expphi_i;
double tmp_i = expphim_r * expphi_i + expphim_i * expphi_r;
expphim_r = tmp_r;
@ -488,7 +496,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
double facn = 1.0 / ncount;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m = 0; m < 2 * l + 1; m++) {
for (int m = 0; m < l + 1; m++) {
qnm_r[il][m] *= facn;
qnm_i[il][m] *= facn;
}
@ -500,55 +508,57 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
int jj = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
double qnormfac = sqrt(MY_4PI / (2 * l + 1));
double qm_sum = 0.0;
for (int m = 0; m < 2 * l + 1; m++)
qm_sum += qnm_r[il][m] * qnm_r[il][m] + qnm_i[il][m] * qnm_i[il][m];
qn[jj++] = qnormfac * sqrt(qm_sum);
double qm_sum = qnm_r[il][0] * qnm_r[il][0];
for (int m = 1; m < l + 1; m++)
qm_sum += 2.0 * (qnm_r[il][m] * qnm_r[il][m] + qnm_i[il][m] * qnm_i[il][m]);
qn[jj++] = qnormfac[il] * sqrt(qm_sum);
}
// calculate W_l
if (wlflag) {
int idxcg_count = 0;
int nterms = 0;
int widx_count = 0;
if (wlflag || wlhatflag) {
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
double wlsum = 0.0;
for (int m1 = 0; m1 < 2 * l + 1; m1++) {
for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) {
int m = m1 + m2 - l;
double qm1qm2_r = qnm_r[il][m1] * qnm_r[il][m2] - qnm_i[il][m1] * qnm_i[il][m2];
double qm1qm2_i = qnm_r[il][m1] * qnm_i[il][m2] + qnm_i[il][m1] * qnm_r[il][m2];
wlsum += (qm1qm2_r * qnm_r[il][m] + qm1qm2_i * qnm_i[il][m]) * cglist[idxcg_count];
idxcg_count++;
for (int m1 = -l; m1 <= 0; m1++) {
const int sgn = 1 - 2 * (m1 & 1); // sgn = (-1)^m1
for (int m2 = 0; m2 <= ((-m1) >> 1); m2++) {
const int m3 = -(m1 + m2);
// Loop enforces -L <= m1 <= 0 <= m2 <= m3 <= L, and m1 + m2 + m3 = 0
// For even L, W3j is invariant under permutation of
// (m1, m2, m3) and (m1, m2, m3) -> (-m1, -m2, -m3). The loop
// structure enforces visiting only one member of each
// such symmetry (invariance) group.
// m1 <= 0, and Qlm[-m] = (-1)^m * conjg(Qlm[m])
const double Q1Q2_r =
(qnm_r[il][-m1] * qnm_r[il][m2] + qnm_i[il][-m1] * qnm_i[il][m2]) * sgn;
const double Q1Q2_i =
(qnm_r[il][-m1] * qnm_i[il][m2] - qnm_i[il][-m1] * qnm_r[il][m2]) * sgn;
const double Q1Q2Q3 = Q1Q2_r * qnm_r[il][m3] - Q1Q2_i * qnm_i[il][m3];
const double c = w3jlist[widx_count++];
wlsum += Q1Q2Q3 * c;
}
}
qn[jj++] = wlsum / sqrt(2 * l + 1);
qn[jj++] = wlsum / qnormfac2[il];
nterms++;
}
}
// calculate W_l_hat
if (wlhatflag) {
int idxcg_count = 0;
const int jptr = jj - nterms;
if (!wlflag) jj = jptr;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
double wlsum = 0.0;
for (int m1 = 0; m1 < 2 * l + 1; m1++) {
for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) {
int m = m1 + m2 - l;
double qm1qm2_r = qnm_r[il][m1] * qnm_r[il][m2] - qnm_i[il][m1] * qnm_i[il][m2];
double qm1qm2_i = qnm_r[il][m1] * qnm_i[il][m2] + qnm_i[il][m1] * qnm_r[il][m2];
wlsum += (qm1qm2_r * qnm_r[il][m] + qm1qm2_i * qnm_i[il][m]) * cglist[idxcg_count];
idxcg_count++;
}
}
if (qn[il] < QEPSILON)
qn[jj++] = 0.0;
else {
double qnormfac = sqrt(MY_4PI / (2 * l + 1));
double qnfac = qnormfac / qn[il];
qn[jj++] = wlsum / sqrt(2 * l + 1) * (qnfac * qnfac * qnfac);
double qnfac = qnormfac[il] / qn[il];
qn[jj++] = qn[jptr + il] * (qnfac * qnfac * qnfac) * qnormfac2[il];
}
}
}
@ -564,9 +574,15 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[],
qn[jj++] = 0.0;
}
else {
double qnormfac = sqrt(MY_4PI / (2 * l + 1));
double qnfac = qnormfac / qn[il];
for (int m = 0; m < 2 * l + 1; m++) {
double qnfac = qnormfac[il] / qn[il];
for (int m = -l; m < 0; m++) {
// Computed only qnm for m>=0.
// qnm[-m] = (-1)^m * conjg(qnm[m])
const int sgn = 1 - 2 * (m & 1); // sgn = (-1)^m
qn[jj++] = qnm_r[il][-m] * qnfac * sgn;
qn[jj++] = -qnm_i[il][-m] * qnfac * sgn;
}
for (int m = 0; m < l + 1; m++) {
qn[jj++] = qnm_r[il][m] * qnfac;
qn[jj++] = qnm_i[il][m] * qnfac;
}
@ -621,68 +637,95 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
}
/* ----------------------------------------------------------------------
assign Clebsch-Gordan coefficients
using the quasi-binomial formula VMK 8.2.1(3)
specialized for case j1=j2=j=l
Initialize table of Wigner 3j symbols
------------------------------------------------------------------------- */
void ComputeOrientOrderAtom::init_clebsch_gordan()
void ComputeOrientOrderAtom::init_wigner3j()
{
double sum, dcg, sfaccg, sfac1, sfac2;
int m, aa2, bb2, cc2;
int ifac, idxcg_count;
int widx_count = 0;
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m1 = 0; m1 < 2 * l + 1; m1++)
for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) idxcg_count++;
const int l = qlist[il];
for (int m1 = -l; m1 <= 0; m1++) {
for (int m2 = 0; m2 <= ((-m1) >> 1); m2++) { widx_count++; }
}
}
idxcg_max = idxcg_count;
memory->destroy(cglist);
memory->create(cglist, idxcg_max, "computeorientorderatom:cglist");
widx_max = widx_count;
memory->destroy(w3jlist);
memory->create(w3jlist, widx_max, "computeorientorderatom:w3jlist");
widx_count = 0;
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for (int m1 = 0; m1 < 2 * l + 1; m1++) {
aa2 = m1 - l;
for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) {
bb2 = m2 - l;
m = aa2 + bb2 + l;
const int l = qlist[il];
// clang-format off
for (int m1 = -l; m1 <= 0; m1++) {
for (int m2 = 0; m2 <= ((-m1) >> 1); m2++) {
const int m3 = -(m1 + m2);
// Loop enforces -L<=m1<=0<=m2<=m3<=L, and m1+m2+m3=0
sum = 0.0;
for (int z = MAX(0, MAX(-aa2, bb2));
z <= MIN(l, MIN(l - aa2, l + bb2)); z++) {
ifac = z % 2 ? -1 : 1;
sum += ifac /
(factorial(z) *
factorial(l - z) *
factorial(l - aa2 - z) *
factorial(l + bb2 - z) *
factorial(aa2 + z) *
factorial(-bb2 + z));
}
// For even L, W3j is invariant under permutation of
// (m1,m2,m3) and (m1,m2,m3)->(-m1,-m2,-m3). The loop
// structure enforces visiting only one member of each
// such symmetry (invariance) group.
cc2 = m - l;
sfaccg = sqrt(factorial(l + aa2) *
factorial(l - aa2) *
factorial(l + bb2) *
factorial(l - bb2) *
factorial(l + cc2) *
factorial(l - cc2) *
(2*l + 1));
// clang-format on
// Determine number of elements in symmetry group of (m1,m2,m3)
// Concise determination exploiting (m1,m2,m3) loop structure.
int pfac;
if (m1 == 0)
pfac = 1; // m1 = m2 = m3 = 0
else if (m2 == 0 || m2 == m3) {
// reduced group when only 3 permutations, or sign inversion
// is equivalent to permutation
pfac = 6;
} else
pfac = 12; // 6 permutations * 2 signs
sfac1 = factorial(3 * l + 1);
sfac2 = factorial(l);
dcg = sqrt(sfac2 * sfac2 * sfac2 / sfac1);
cglist[idxcg_count] = sum * dcg * sfaccg;
idxcg_count++;
w3jlist[widx_count] = w3j(l, m1, m2, m3) * pfac;
widx_count++;
}
}
}
}
/* ---------------------------------------------------------------------- */
double ComputeOrientOrderAtom::triangle_coeff(const int a, const int b, const int c)
{
return factorial(a + b - c) * factorial(a - b + c) * factorial(-a + b + c) /
factorial(a + b + c + 1);
}
/* ---------------------------------------------------------------------- */
double ComputeOrientOrderAtom::w3j(const int lmax, const int j1, const int j2, const int j3)
{
const int a = lmax, b = lmax, c = lmax;
const int alpha = j1, beta = j2, gamma = j3;
struct {
double operator()(const int a, const int b, const int c, const int alpha, const int beta,
const int t)
{
return factorial(t) * factorial(c - b + t + alpha) * factorial(c - a + t - beta) *
factorial(a + b - c - t) * factorial(a - t - alpha) * factorial(b - t + beta);
}
} x;
const double sgn = 1 - 2 * ((a - b - gamma) & 1);
const double g = sqrt(triangle_coeff(lmax, lmax, lmax)) *
sqrt(factorial(a + alpha) * factorial(a - alpha) * factorial(b + beta) * factorial(b - beta) *
factorial(c + gamma) * factorial(c - gamma));
double s = 0;
int t = 0;
while (c - b + t + alpha < 0 || c - a + t - beta < 0) t++;
// ^^ t>=-j1 ^^ t>=j2
while (1) {
if (a + b - c - t < 0) break; // t<=lmax
if (a - t - alpha < 0) break; // t<=lmax-j1
if (b - t + beta < 0) break; // t<=lmax+j2
const int m1t = 1 - 2 * (t & 1);
s += m1t / x(lmax, lmax, lmax, alpha, beta, t);
t++;
}
return sgn * g * s;
}

9
src/compute_orientorder_atom.h
View File

@ -36,6 +36,7 @@ class ComputeOrientOrderAtom : public Compute {
int iqlcomp, qlcomp, qlcompflag, wlflag, wlhatflag;
int *qlist;
int nqlist;
double *qnormfac, *qnormfac2;
protected:
int nmax, maxneigh, ncol, nnn;
@ -54,9 +55,11 @@ class ComputeOrientOrderAtom : public Compute {
double polar_prefactor(int, int, double);
double associated_legendre(int, int, double);
virtual void init_clebsch_gordan();
double *cglist; // Clebsch-Gordan coeffs
int idxcg_max;
virtual void init_wigner3j();
double triangle_coeff(const int a, const int b, const int c);
double w3j(const int L, const int j1, const int j2, const int j3);
double *w3jlist; // Wigner coeffs
int widx_max;
int chunksize;
};

1
src/dump.cpp
View File

@ -1086,6 +1086,7 @@ void Dump::modify_params(int narg, char **arg)
}
output->mode_dump[idump] = 1;
output->every_time_dump[idump] = delta;
output->next_dump[idump] = update->ntimestep;
iarg += 2;
} else if (strcmp(arg[iarg],"fileper") == 0) {

2
src/finish.cpp
View File

@ -483,7 +483,7 @@ void Finish::end(int flag)
if (neighflag) {
if (me == 0) utils::logmesg(lmp,"\n");
tmp = MAX(nneigh,nneighfull);
tmp = MAX(MAX(nneigh,nneighfull),0.0);
double nall;
MPI_Allreduce(&tmp,&nall,1,MPI_DOUBLE,MPI_SUM,world);

2
src/info.cpp
View File

@ -110,7 +110,7 @@ using namespace LAMMPS_NS;
static const char *varstyles[] = {
"index", "loop", "world", "universe", "uloop", "string", "getenv",
"file", "atomfile", "format", "equal", "atom", "vector", "python",
"internal", "(unknown)"};
"timer", "internal", "(unknown)"};
static const char *mapstyles[] = { "none", "array", "hash", "yes" };

10
src/library.cpp
View File

@ -4850,8 +4850,9 @@ int lammps_id_name(void *handle, const char *category, int idx, char *buffer, in
auto lmp = (LAMMPS *) handle;
if (strcmp(category,"compute") == 0) {
if ((idx >=0) && (idx < lmp->modify->ncompute)) {
strncpy(buffer, lmp->modify->compute[idx]->id, buf_size);
auto icompute = lmp->modify->get_compute_by_index(idx);
if (icompute) {
strncpy(buffer, icompute->id, buf_size);
return 1;
}
} else if (strcmp(category,"dump") == 0) {
@ -4860,8 +4861,9 @@ int lammps_id_name(void *handle, const char *category, int idx, char *buffer, in
return 1;
}
} else if (strcmp(category,"fix") == 0) {
if ((idx >=0) && (idx < lmp->modify->nfix)) {
strncpy(buffer, lmp->modify->fix[idx]->id, buf_size);
auto ifix = lmp->modify->get_fix_by_index(idx);
if (ifix) {
strncpy(buffer, ifix->id, buf_size);
return 1;
}
} else if (strcmp(category,"group") == 0) {

54
src/lmpfftsettings.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
// common FFT library related defines and compilation settings
#ifndef LMP_FFT_SETTINGS_H
#define LMP_FFT_SETTINGS_H
// if user set FFTW, it means FFTW3
#ifdef FFT_FFTW
#ifndef FFT_FFTW3
#define FFT_FFTW3
#endif
#endif
// set strings for library info output
#if defined(FFT_FFTW3)
#define LMP_FFT_LIB "FFTW3"
#elif defined(FFT_MKL)
#define LMP_FFT_LIB "MKL FFT"
#elif defined(FFT_CUFFT)
#define LMP_FFT_LIB "cuFFT"
#elif defined(FFT_HIPFFT)
#define LMP_FFT_LIB "hipFFT"
#else
#define LMP_FFT_LIB "KISS FFT"
#endif
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define FFT_PRECISION 1
#define LMP_FFT_PREC "single"
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define FFT_PRECISION 2
#define LMP_FFT_PREC "double"
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#endif

8
src/modify.h
View File

@ -113,7 +113,9 @@ class Modify : protected Pointers {
int find_fix(const std::string &);
// new API
Fix *get_fix_by_id(const std::string &) const;
Fix *get_fix_by_index(int idx) const { return fix[idx]; }
Fix *get_fix_by_index(int idx) const {
return ((idx >= 0) && (idx < nfix)) ? fix[idx] : nullptr;
}
const std::vector<Fix *> get_fix_by_style(const std::string &) const;
const std::vector<Fix *> &get_fix_list();
@ -127,7 +129,9 @@ class Modify : protected Pointers {
int find_compute(const std::string &);
// new API
Compute *get_compute_by_id(const std::string &) const;
Compute *get_compute_by_index(int idx) const { return compute[idx]; }
Compute *get_compute_by_index(int idx) const {
return ((idx >= 0) && (idx < ncompute)) ? compute[idx] : nullptr;
}
const std::vector<Compute *> get_compute_by_style(const std::string &) const;
const std::vector<Compute *> &get_compute_list();

10
src/output.cpp
View File

@ -495,8 +495,7 @@ void Output::calculate_next_dump(int which, int idump, bigint ntimestep)
// which = WRITE: increment nextdump by every_dump
if (which == SETUP)
next_dump[idump] =
(ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump];
next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump];
else if (which == WRITE)
next_dump[idump] += every_dump[idump];
@ -754,16 +753,14 @@ void Output::add_dump(int narg, char **arg)
if (ndump == max_dump) {
max_dump += DELTA;
dump = (Dump **)
memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump");
dump = (Dump **) memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump");
memory->grow(mode_dump,max_dump,"output:mode_dump");
memory->grow(every_dump,max_dump,"output:every_dump");
memory->grow(every_time_dump,max_dump,"output:every_time_dump");
memory->grow(next_dump,max_dump,"output:next_dump");
memory->grow(next_time_dump,max_dump,"output:next_time_dump");
memory->grow(last_dump,max_dump,"output:last_dump");
var_dump = (char **)
memory->srealloc(var_dump,max_dump*sizeof(char *),"output:var_dump");
var_dump = (char **) memory->srealloc(var_dump,max_dump*sizeof(char *),"output:var_dump");
memory->grow(ivar_dump,max_dump,"output:ivar_dump");
}
@ -784,6 +781,7 @@ void Output::add_dump(int narg, char **arg)
last_dump[ndump] = -1;
var_dump[ndump] = nullptr;
ivar_dump[ndump] = -1;
next_dump[ndump] = 0;
ndump++;
}

49
src/variable.cpp
View File

@ -500,6 +500,31 @@ void Variable::set(int narg, char **arg)
strcpy(data[nvar][1],"(undefined)");
}
// TIMER
// stores current walltime as a timestamp in seconds
// replace pre-existing var if also style TIMER (allows reset with current time)
// num = 1, for string representation of dvalue, set by retrieve()
// dvalue = numeric initialization via platform::walltime()
} else if (strcmp(arg[1],"timer") == 0) {
if (narg != 2) error->all(FLERR,"Illegal variable command");
int ivar = find(arg[0]);
if (ivar >= 0) {
if (style[ivar] != TIMER)
error->all(FLERR,"Cannot redefine variable as a different style");
dvalue[ivar] = platform::walltime();
replaceflag = 1;
} else {
if (nvar == maxvar) grow();
style[nvar] = TIMER;
num[nvar] = 1;
which[nvar] = 0;
pad[nvar] = 0;
data[nvar] = new char*[num[nvar]];
data[nvar][0] = new char[VALUELENGTH];
dvalue[nvar] = platform::walltime();
}
// INTERNAL
// replace pre-existing var if also style INTERNAL (allows it to be reset)
// num = 1, for string representation of dvalue, set by retrieve()
@ -524,6 +549,8 @@ void Variable::set(int narg, char **arg)
dvalue[nvar] = utils::numeric(FLERR,arg[2],false,lmp);
}
// unrecognized variable style
} else error->all(FLERR,"Illegal variable command");
// set name of variable, if not replacing one flagged with replaceflag
@ -610,13 +637,14 @@ int Variable::next(int narg, char **arg)
error->all(FLERR,"All variables in next command must have same style");
}
// invalid styles: STRING, EQUAL, WORLD, ATOM, VECTOR, GETENV,
// FORMAT, PYTHON, INTERNAL
// invalid styles: STRING, EQUAL, WORLD, GETENV, ATOM, VECTOR,
// FORMAT, PYTHON, TIMER, INTERNAL
int istyle = style[find(arg[0])];
if (istyle == STRING || istyle == EQUAL || istyle == WORLD ||
istyle == GETENV || istyle == ATOM || istyle == VECTOR ||
istyle == FORMAT || istyle == PYTHON || istyle == INTERNAL)
if (istyle == STRING || istyle == EQUAL ||
istyle == WORLD || istyle == GETENV || istyle == ATOM ||
istyle == VECTOR || istyle == FORMAT || istyle == PYTHON ||
istyle == TIMER || istyle == INTERNAL)
error->all(FLERR,"Invalid variable style with next command");
// if istyle = UNIVERSE or ULOOP, insure all such variables are incremented
@ -794,13 +822,15 @@ void Variable::python_command(int narg, char **arg)
}
/* ----------------------------------------------------------------------
return 1 if variable is EQUAL or INTERNAL or PYTHON numeric style, 0 if not
return 1 if variable is EQUAL style, 0 if not
TIMER, INTERNAL, PYTHON qualify as EQUAL style
this is checked before call to compute_equal() to return a double
------------------------------------------------------------------------- */
int Variable::equalstyle(int ivar)
{
if (style[ivar] == EQUAL || style[ivar] == INTERNAL) return 1;
if (style[ivar] == EQUAL || style[ivar] == TIMER ||
style[ivar] == INTERNAL) return 1;
if (style[ivar] == PYTHON) {
int ifunc = python->variable_match(data[ivar][0],names[ivar],1);
if (ifunc < 0) return 0;
@ -925,7 +955,7 @@ char *Variable::retrieve(const char *name)
// then the Python class stores the result, query it via long_string()
char *strlong = python->long_string(ifunc);
if (strlong) str = strlong;
} else if (style[ivar] == INTERNAL) {
} else if (style[ivar] == TIMER || style[ivar] == INTERNAL) {
sprintf(data[ivar][0],"%.15g",dvalue[ivar]);
str = data[ivar][0];
} else if (style[ivar] == ATOM || style[ivar] == ATOMFILE ||
@ -938,7 +968,7 @@ char *Variable::retrieve(const char *name)
/* ----------------------------------------------------------------------
return result of equal-style variable evaluation
can be EQUAL or INTERNAL style or PYTHON numeric style
can be EQUAL or TIMER or INTERNAL style or PYTHON numeric style
for PYTHON, don't need to check python->variable_match() error return,
since caller will have already checked via equalstyle()
------------------------------------------------------------------------- */
@ -952,6 +982,7 @@ double Variable::compute_equal(int ivar)
double value = 0.0;
if (style[ivar] == EQUAL) value = evaluate(data[ivar][0],nullptr,ivar);
else if (style[ivar] == TIMER) value = dvalue[ivar];
else if (style[ivar] == INTERNAL) value = dvalue[ivar];
else if (style[ivar] == PYTHON) {
int ifunc = python->find(data[ivar][0]);

1
src/variable.h
View File

@ -71,6 +71,7 @@ class Variable : protected Pointers {
ATOM,
VECTOR,
PYTHON,
TIMER,
INTERNAL
};
static constexpr int VALUELENGTH = 64;