jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://github.com/lammps/lammps into gridcomm-tiled

This commit is contained in:
Stan Moore 2020-08-18 13:43:33 -06:00
parent 71e0646063 d63f3d87eb
commit e2923d2d8c
172 changed files with 9786 additions and 939 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/Documentation.cmake
View File

@ -29,8 +29,8 @@ if(BUILD_DOC)
OUTPUT requirements.txt
DEPENDS docenv
COMMAND ${CMAKE_COMMAND} -E copy ${LAMMPS_DOC_DIR}/utils/requirements.txt requirements.txt
COMMAND ${DOCENV_BINARY_DIR}/pip install -r requirements.txt --upgrade
COMMAND ${DOCENV_BINARY_DIR}/pip install --upgrade ${LAMMPS_DOC_DIR}/utils/converters
COMMAND ${DOCENV_BINARY_DIR}/pip install --use-feature=2020-resolver -r requirements.txt --upgrade
)
# download mathjax distribution and unpack to folder "mathjax"

2
doc/Makefile
View File

@ -194,7 +194,7 @@ $(VENV):
$(VIRTUALENV) -p $(PYTHON) $(VENV); \
. $(VENV)/bin/activate; \
pip install --upgrade pip; \
pip install -r requirements.txt; \
pip install --use-feature=2020-resolver -r requirements.txt; \
deactivate;\
)

60
doc/src/Build_development.rst
View File

@ -159,8 +159,8 @@ A test run is then a a collection multiple individual test runs each
with many comparisons to reference results based on template input
files, individual command settings, relative error margins, and
reference data stored in a YAML format file with ``.yaml``
suffix. Currently the programs ``pair_style``, ``bond_style``, and
``angle_style`` are implemented. They will compare forces, energies and
suffix. Currently the programs ``test_pair_style``, ``test_bond_style``, and
``test_angle_style`` are implemented. They will compare forces, energies and
(global) stress for all atoms after a ``run 0`` calculation and after a
few steps of MD with :doc:`fix nve <fix_nve>`, each in multiple variants
with different settings and also for multiple accelerated styles. If a
@ -172,7 +172,7 @@ Below is an example command and output:
.. parsed-literal::
[tests]$ pair_style mol-pair-lj_cut.yaml
[tests]$ test_pair_style mol-pair-lj_cut.yaml
[==========] Running 6 tests from 1 test suite.
[----------] Global test environment set-up.
[----------] 6 tests from PairStyle
@ -259,13 +259,61 @@ and working.
of mis-compiled code (or an undesired large loss of precision due
to significant reordering of operations and thus less error cancellation).
Unit tests for timestepping related fixes
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A substantial subset of :doc:`fix styles <fix>` are invoked regularly
during MD timestepping and manipulate per-atom properties like
positions, velocities, and forces. For those fix styles, testing can be
done in a very similar fashion as for force fields and thus there is a
test program `test_fix_timestep` that shares a lot of code, properties,
and command line flags with the force field style testers described in
the previous section.
This tester will set up a small molecular system run with verlet run
style for 4 MD steps, then write a binary restart and continue for
another 4 MD steps. At this point coordinates and velocities are
recorded and compared to reference data. Then the system is cleared,
restarted and running the second 4 MD steps again and the data is
compared to the same reference. That is followed by another restart
after which per atom type masses are replaced with per-atom masses and
the second 4 MD steps are repeated again and compared to the same
reference. Also global scalar and vector data of the fix is recorded
and compared. If the fix is a thermostat and thus the internal property
``t_target`` can be extracted, then this is compared to the reference
data. The tests are repeated with the respa run style.
If the fix has a multi-threaded version in the USER-OMP package, then
the entire set of tests is repeated for that version as well.
For this to work, some additional conditions have to be met by the
YAML format test inputs.
- The fix to be tested (and only this fix), should be listed in the
``prerequisites:`` section
- The fix to be tested must be specified in the ``post_commands:``
section with the fix-ID ``test``. This section may contain other
commands and other fixes (e.g. an instance of fix nve for testing
a thermostat or force manipulation fix)
- For fixes that can tally contributions to the global virial, the
line ``fix_modify test virial yes`` should be included in the
``post_commands:`` section of the test input.
- For thermostat fixes the target temperature should be ramped from
an arbitrary value (e.g. 50K) to a pre-defined target temperature
entered as ``${t_target}``.
- For fixes that have thermostatting support included, but do not
have it enabled in the input (e.g. fix rigid with default settings),
the ``post_commands:`` section should contain the line
``variable t_target delete`` to disable the target temperature ramp
check to avoid false positives.
Use custom linker for faster link times when ENABLE_TESTING is active
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
When compiling LAMMPS with enabled tests, most test executables will
need to be linked against the LAMMPS library. Since this can be a
large with many C++ objects when many packages are enabled, link times
can become very long on machines that use the GNU BFD linker (e.g.
need to be linked against the LAMMPS library. Since this can be a very
large library with many C++ objects when many packages are enabled, link
times can become very long on machines that use the GNU BFD linker (e.g.
Linux systems). Alternatives like the ``lld`` linker of the LLVM project
or the ``gold`` linker available with GNU binutils can speed up this step
substantially. CMake will by default test if any of the two can be

1
doc/src/Commands_bond.rst
View File

@ -46,6 +46,7 @@ OPT.
* :doc:`oxdna2/fene <bond_oxdna>`
* :doc:`oxrna2/fene <bond_oxdna>`
* :doc:`quartic (o) <bond_quartic>`
* :doc:`special <bond_special>`
* :doc:`table (o) <bond_table>`
.. _angle:

2
doc/src/bond_none.rst
View File

@ -13,7 +13,7 @@ Syntax
Examples
""""""""
.. code-blocK:: LAMMPS
.. code-block:: LAMMPS
bond_style none

106
doc/src/bond_special.rst Normal file
View File

@ -0,0 +1,106 @@
.. index:: bond_style special
bond_style special command
=================================
Syntax
""""""
.. code-block:: LAMMPS
bond_style special
Examples
""""""""
.. code-block:: LAMMPS
bond_style special
bond_coeff 0.5 0.5
Description
"""""""""""
The *special* bond style can be used to create conceptual bonds which
effectively impose weightings on the pairwise Lennard Jones and/or
Coulombic interactions between selected pairs of particles in the
system. The form of the pairwise interaction will be whatever is
computed by the :doc:`pair_style <pair_style>` command defined for the
system; this command defines the weightings for its two terms.
This command can thus be useful to apply weightings that cannot be
handled by the :doc:`special_bonds <special_bonds>` command, such as
on 1-5 or 1-6 interactions. Or it can be used to add pairwise forces
between one or more pairs of atoms that otherwise would not be include
in the :doc:`pair_style <pair_style>` computation.
The potential for this bond style has the form
.. math::
E = w_{LJ} E_{LJ} + w_{Coul} E_{Coul}
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
the data file or restart files read by the :doc:`read_data <read_data>`
or :doc:`read_restart <read_restart>` commands:
* :math:`w_{LJ}` weight (0.0 to 1.0) on pairwise Lennard-Jones interactions
* :math:`w_{Coul}` weight (0.0 to 1.0) on pairwise Coulombic interactions
----------
Normally this bond style should be used in conjunction with one (or
more) other bond styles which compute forces between atoms directly
bonded to each other in a molecule. This means the :doc:`bond_style
hybrid <bond_hybrid>` command should be used with bond_style special
as one of its sub-styles.
Note that the same as for any other bond style, pairs of bonded atoms
must be enumerated in the data file read by the :doc:`read_data
<read_data>` command. Thus if this command is used to weight all 1-5
interactions in the system, all the 1-5 pairs of atoms must be listed
in the "Bonds" section of the data file.
This bond style imposes strict requirements on settings made with the
:doc:`special_bonds <special_bonds>` command. These requirements
ensure that the new bonds created by this style do not create spurious
1-2, 1-3, or 1-4 interactions within the molecular topology.
Specifically 1-2 interactions must have weights of zero, 1-3
interactions must either have weights of unity or :doc:`special_bonds
angle yes <special_bonds>` must be used, and 1-4 interactions must
have weights of unity or :doc:`special_bonds dihedral yes <special_bonds>`
must be used.
If this command is used to create bonded interactions between
particles that are further apart than usual (e.g. 1-5 or 1-6
interactions), this style may require an increase in the communication
cutoff via the :doc:`comm_modify cutoff <comm_modify>` command. If
LAMMPS cannot find a partner atom in a bond, an error will be issued.
Restrictions
""""""""""""
This bond style can only be used if LAMMPS was built with the
USER-MISC package. See the :doc:`Build package <Build_package>` doc
page for more info.
This bond style requires the use of a :doc:`pair_style <pair_style>` which
computes a pairwise additive interaction and provides the ability to
compute interactions for individual pairs of atoms. Manybody potentials
are not compatible in general, but also some other pair styles are missing
the required functionality and thus will cause an error.
This command is not compatible with long-range Coulombic interactions. If a
`kspace_style <kspace_style>` is declared, an error will be issued.
Related commands
""""""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`special_bonds <special_bonds>`
**Default:** none

8
doc/src/compute_saed.rst
View File

@ -57,8 +57,8 @@ is computed from the structure factor F using the equations:
.. math::
I = & \frac{F^{*}F}{N} \\
F(\mathbf{k}) = & \sum_{j=1}^{N}f_j(\theta)exp(2\pi i \mathbf{k} \cdot \mathbf{r}_j)
I = & \frac{F^{*}F}{N} \\
F(\mathbf{k}) = & \sum_{j=1}^{N}f_j(\theta)exp(2\pi i \mathbf{k} \cdot \mathbf{r}_j)
Here, K is the location of the reciprocal lattice node, :math:`r_j` is the
position of each atom, :math:`f_j` are atomic scattering factors.
@ -116,8 +116,8 @@ The analytic approximation is computed using the formula
.. math::
f_j\left ( \frac{sin(\theta)}{\lambda} \right )=\sum_{i}^{5}
a_i exp\left ( -b_i \frac{sin^{2}(\theta)}{\lambda^{2}} \right )
f_j\left ( \frac{sin(\theta)}{\lambda} \right )=\sum_{i}^{5}
a_i exp\left ( -b_i \frac{sin^{2}(\theta)}{\lambda^{2}} \right )
Coefficients parameterized by :ref:`(Fox) <Fox>` are assigned for each
atom type designating the chemical symbol and charge of each atom

6
doc/src/fix_temp_berendsen.rst
View File

@ -119,7 +119,11 @@ thermal degrees of freedom, and the bias is added back in.
**Restart, fix_modify, output, run start/stop, minimize info:**
No information about this fix is written to :doc:`binary restart files <restart>`.
This fix writes the cumulative global energy change to
:doc:`binary restart files <restart>`. See the
:doc:`read_restart <read_restart>` command for info on how to
re-specify a fix in an input script that reads a restart file,
so that the fix continues in an uninterrupted fashion.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
fix. You can use it to assign a temperature :doc:`compute <compute>`

8
doc/src/fix_temp_csvr.rst
View File

@ -127,6 +127,14 @@ thermal degrees of freedom, and the bias is added back in.
**Restart, fix_modify, output, run start/stop, minimize info:**
These fixes write the cumulative global energy change and the
random number generator states to :doc:`binary restart files <restart>`.
See the :doc:`read_restart <read_restart>` command for info on how to
re-specify a fix in an input script that reads a restart file,
so that the selected fix continues in an uninterrupted fashion. The
random number generator state can only be restored when the number
of processors remains unchanged from what is recorded in the restart file.
No information about these fixes are written to :doc:`binary restart files <restart>`.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by these

30
doc/src/lattice.rst
View File

@ -260,25 +260,21 @@ of Aidan Thompson), with its 8 atom unit cell.
variable b equal $a*sqrt(3.0)
variable c equal $a*sqrt(8.0/3.0)
variable 1_3 equal 1.0/3.0
variable 2_3 equal 2.0/3.0
variable 1_6 equal 1.0/6.0
variable 5_6 equal 5.0/6.0
variable 1_12 equal 1.0/12.0
variable 5_12 equal 5.0/12.0
variable third equal 1.0/3.0
variable five6 equal 5.0/6.0
lattice custom 1.0 &
a1 $a 0.0 0.0 &
a2 0.0 $b 0.0 &
a3 0.0 0.0 $c &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.0 &
basis ${1_3} 0.0 0.5 &
basis ${5_6} 0.5 0.5 &
basis 0.0 0.0 0.625 &
basis 0.5 0.5 0.625 &
basis ${1_3} 0.0 0.125 &
basis ${5_6} 0.5 0.125
a1 $a 0.0 0.0 &
a2 0.0 $b 0.0 &
a3 0.0 0.0 $c &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.0 &
basis ${third} 0.0 0.5 &
basis ${five6} 0.5 0.5 &
basis 0.0 0.0 0.625 &
basis 0.5 0.5 0.625 &
basis ${third} 0.0 0.125 &
basis ${five6} 0.5 0.125
region myreg block 0 1 0 1 0 1
create_box 2 myreg

8
doc/src/pair_fep_soft.rst
View File

@ -161,11 +161,11 @@ Examples
pair_style coul/long/soft 1.0 10.0 9.5
pair_coeff * * 1.0
pair_coeff 1 1 1.0 9.5
pair_coeff 1 1 1.0
pair_style tip4p/long/soft 1 2 7 8 0.15 2.0 0.5 10.0 9.8
pair_coeff * * 1.0
pair_coeff 1 1 1.0 9.5
pair_coeff 1 1 1.0
pair_style morse/soft 4 0.9 10.0
pair_coeff * * 100.0 2.0 1.5 1.0
@ -284,7 +284,9 @@ core. Hence, if used by themselves, there will be no repulsion to keep two
oppositely charged particles from overlapping each other. In this case, if
:math:`\lambda = 1`, a singularity may occur. These sub-styles are suitable to
represent charges embedded in the Lennard-Jones radius of another site (for
example hydrogen atoms in several water models).
example hydrogen atoms in several water models). The :math:`\lambda` must
be defined for each pair, and *coul/cut/soft* can accept an optional cutoff as
the second coefficient.
.. note::

2
doc/src/pair_meam_spline.rst
View File

@ -16,7 +16,7 @@ Syntax
Examples
""""""""
.. code:: LAMMPS
.. code-block:: LAMMPS
pair_style meam/spline
pair_coeff * * Ti.meam.spline Ti

2
doc/src/pair_table.rst
View File

@ -60,7 +60,7 @@ For the *linear* style, the distance *R* is used to find the 2
surrounding table values from which an energy or force is computed by
linear interpolation.
For the *spline* style, a cubic spline coefficients are computed and
For the *spline* style, cubic spline coefficients are computed and
stored for each of the *N* values in the table, one set of splines for
energy, another for force. Note that these splines are different than
the ones used to pre-compute the *N* values. Those splines were fit

4
src/GRANULAR/fix_wall_gran.cpp
View File

@ -1541,6 +1541,7 @@ int FixWallGran::pack_restart(int i, double *buf)
if (!use_history) return 0;
int n = 0;
// pack buf[0] this way because other fixes unpack it
buf[n++] = size_history + 1;
for (int m = 0; m < size_history; m++)
buf[n++] = history_one[i][m];
@ -1558,7 +1559,8 @@ void FixWallGran::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;

4
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -479,6 +479,7 @@ int FixWallGranRegion::pack_restart(int i, double *buf)
for (m = 0; m < size_history; m++)
buf[n++] = history_many[i][iwall][m];
}
// pack buf[0] this way because other fixes unpack it
buf[0] = n;
return n;
}
@ -496,7 +497,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;

13
src/KOKKOS/kokkos_type.h
View File

@ -1058,7 +1058,7 @@ struct alignas(2*sizeof(real)) SNAComplex
{
real re,im;
KOKKOS_FORCEINLINE_FUNCTION SNAComplex() = default;
SNAComplex() = default;
KOKKOS_FORCEINLINE_FUNCTION SNAComplex(real re)
: re(re), im(static_cast<real>(0.)) { ; }
@ -1100,6 +1100,17 @@ KOKKOS_FORCEINLINE_FUNCTION SNAComplex<real> operator*(const real& r, const SNAC
typedef SNAComplex<SNAreal> SNAcomplex;
// Cayley-Klein pack
// Can guarantee it's aligned to 2 complex
struct alignas(32) CayleyKleinPack {
SNAcomplex a, b;
SNAcomplex da[3], db[3];
SNAreal sfac;
SNAreal dsfacu[3];
};
#if defined(KOKKOS_ENABLE_CXX11)
#undef ISFINITE

9
src/KOKKOS/pair_snap_kokkos.h
View File

@ -50,6 +50,7 @@ struct TagPairSNAPComputeFusedDeidrj{};
// CPU backend only
struct TagPairSNAPPreUiCPU{};
struct TagPairSNAPComputeUiCPU{};
struct TagPairSNAPTransformUiCPU{};
struct TagPairSNAPComputeZiCPU{};
struct TagPairSNAPBetaCPU{};
struct TagPairSNAPComputeBiCPU{};
@ -104,7 +105,7 @@ public:
void operator() (TagPairSNAPComputeUi,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeUi>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPTransformUi,const int iatom_mod, const int idxu, const int iatom_div) const;
void operator() (TagPairSNAPTransformUi,const int iatom_mod, const int j, const int iatom_div) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeZi,const int iatom_mod, const int idxz, const int iatom_div) const;
@ -134,15 +135,15 @@ public:
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeUiCPU,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeUiCPU>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPTransformUiCPU, const int j, const int iatom) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeZiCPU,const int& ii) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeBiCPU,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeBiCPU>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPZeroYiCPU,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPZeroYiCPU>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeYiCPU,const int& ii) const;

157
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -206,8 +206,6 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
EV_FLOAT ev;
int idxu_max = snaKK.idxu_max;
while (chunk_offset < inum) { // chunk up loop to prevent running out of memory
EV_FLOAT ev_tmp;
@ -246,6 +244,13 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
Kokkos::parallel_for("ComputeUiCPU",policy_ui_cpu,*this);
}
{
// Expand ulisttot -> ulisttot_full
// Zero out ylist
typename Kokkos::MDRangePolicy<DeviceType, Kokkos::IndexType<int>, Kokkos::Rank<2, Kokkos::Iterate::Left, Kokkos::Iterate::Left>, TagPairSNAPTransformUiCPU> policy_transform_ui_cpu({0,0},{twojmax+1,chunk_size});
Kokkos::parallel_for("TransformUiCPU",policy_transform_ui_cpu,*this);
}
//Compute bispectrum
if (quadraticflag || eflag) {
//ComputeZi
@ -261,20 +266,12 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
Kokkos::parallel_for("ComputeBiCPU",policy_bi_cpu,*this);
}
//ZeroYi,ComputeYi
//ComputeYi
{
int vector_length = vector_length_default;
int team_size = team_size_default;
//Compute beta = dE_i/dB_i for all i in list
typename Kokkos::RangePolicy<DeviceType,TagPairSNAPBetaCPU> policy_beta(0,chunk_size);
Kokkos::parallel_for("ComputeBetaCPU",policy_beta,*this);
//ZeroYi
check_team_size_for<TagPairSNAPZeroYiCPU>(chunk_size,team_size,vector_length);
typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPZeroYiCPU> policy_zero_yi(((idxu_max+team_size-1)/team_size)*chunk_size,team_size,vector_length);
Kokkos::parallel_for("ZeroYiCPU",policy_zero_yi,*this);
//ComputeYi
int idxz_max = snaKK.idxz_max;
typename Kokkos::RangePolicy<DeviceType,TagPairSNAPComputeYiCPU> policy_yi_cpu(0,chunk_size*idxz_max);
@ -294,6 +291,7 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
Kokkos::parallel_for("ComputeDeidrjCPU",policy_deidrj_cpu,*this);
}
} else { // GPU
#ifdef LMP_KOKKOS_GPU
@ -313,10 +311,10 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
int team_size = 4; // need to cap b/c of shared memory reqs
check_team_size_for<TagPairSNAPComputeUi>(chunk_size,team_size,vector_length);
// scratch size: 2 * team_size * (twojmax+1)^2, to cover all `m1`,`m2` values
// scratch size: 2 * team_size * (twojmax+1)^2, to cover all `m1`,`m2` values, div 2 for symmetry
// 2 is for double buffer
const int tile_size = (twojmax+1)*(twojmax+1);
const int tile_size = (twojmax+1)*(twojmax/2+1);
typedef Kokkos::View< SNAcomplex*,
Kokkos::DefaultExecutionSpace::scratch_memory_space,
Kokkos::MemoryTraits<Kokkos::Unmanaged> >
@ -329,7 +327,7 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
Kokkos::parallel_for("ComputeUi",policy_ui,*this);
//Transform data layout of ulisttot to AoSoA, zero ylist
typename Kokkos::MDRangePolicy<DeviceType, Kokkos::IndexType<int>, Kokkos::Rank<3, Kokkos::Iterate::Left, Kokkos::Iterate::Left>, TagPairSNAPTransformUi> policy_transform_ui({0,0,0},{32,idxu_max,(chunk_size + 32 - 1) / 32},{32,4,1});
typename Kokkos::MDRangePolicy<DeviceType, Kokkos::IndexType<int>, Kokkos::Rank<3, Kokkos::Iterate::Left, Kokkos::Iterate::Left>, TagPairSNAPTransformUi> policy_transform_ui({0,0,0},{32,twojmax+1,(chunk_size + 32 - 1) / 32},{32,4,1});
Kokkos::parallel_for("TransformUi",policy_transform_ui,*this);
}
@ -367,7 +365,8 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
Kokkos::parallel_for("ComputeYi",policy_compute_yi,*this);
//Transform data layout of ylist out of AoSoA
typename Kokkos::MDRangePolicy<DeviceType, Kokkos::IndexType<int>, Kokkos::Rank<3, Kokkos::Iterate::Left, Kokkos::Iterate::Left>, TagPairSNAPTransformYi> policy_transform_yi({0,0,0},{32,idxu_max,(chunk_size + 32 - 1) / 32},{32,4,1});
const int idxu_half_max = snaKK.idxu_half_max;
typename Kokkos::MDRangePolicy<DeviceType, Kokkos::IndexType<int>, Kokkos::Rank<3, Kokkos::Iterate::Left, Kokkos::Iterate::Left>, TagPairSNAPTransformYi> policy_transform_yi({0,0,0},{32,idxu_half_max,(chunk_size + 32 - 1) / 32},{32,4,1});
Kokkos::parallel_for("TransformYi",policy_transform_yi,*this);
}
@ -397,7 +396,7 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
}
}
#endif // KOKKOS_ENABLE_CUDA
#endif // LMP_KOKKOS_GPU
}
@ -608,12 +607,21 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeNeigh,const typen
if ( rsq < rnd_cutsq(itype,jtype) ) {
if (final) {
my_sna.rij(ii,offset,0) = dx;
my_sna.rij(ii,offset,1) = dy;
my_sna.rij(ii,offset,2) = dz;
#ifdef LMP_KOKKOS_GPU
if (std::is_same<DeviceType,Kokkos::Cuda>::value) {
my_sna.compute_cayley_klein(ii, offset, dx, dy, dz, (radi + d_radelem[elem_j])*rcutfac,
d_wjelem[elem_j]);
} else {
#endif
my_sna.rij(ii,offset,0) = dx;
my_sna.rij(ii,offset,1) = dy;
my_sna.rij(ii,offset,2) = dz;
my_sna.wj(ii,offset) = d_wjelem[elem_j];
my_sna.rcutij(ii,offset) = (radi + d_radelem[elem_j])*rcutfac;
#ifdef LMP_KOKKOS_GPU
}
#endif
my_sna.inside(ii,offset) = j;
my_sna.wj(ii,offset) = d_wjelem[elem_j];
my_sna.rcutij(ii,offset) = (radi + d_radelem[elem_j])*rcutfac;
if (chemflag)
my_sna.element(ii,offset) = elem_j;
else
@ -704,27 +712,54 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeUi,const typename
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPTransformUi,const int iatom_mod, const int idxu, const int iatom_div) const {
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPTransformUi,const int iatom_mod, const int j, const int iatom_div) const {
SNAKokkos<DeviceType> my_sna = snaKK;
const int iatom = iatom_mod + iatom_div * 32;
if (iatom >= chunk_size) return;
if (idxu >= my_sna.idxu_max) return;
if (j > twojmax) return;
int elem_count = chemflag ? nelements : 1;
for (int ielem = 0; ielem < elem_count; ielem++) {
const int jju_half = my_sna.idxu_half_block(j);
const int jju = my_sna.idxu_block(j);
const auto utot_re = my_sna.ulisttot_re(idxu, ielem, iatom);
const auto utot_im = my_sna.ulisttot_im(idxu, ielem, iatom);
for (int mb = 0; 2*mb <= j; mb++) {
for (int ma = 0; ma <= j; ma++) {
// Extract top half
my_sna.ulisttot_pack(iatom_mod, idxu, ielem, iatom_div) = { utot_re, utot_im };
const int idxu_shift = mb * (j + 1) + ma;
const int idxu_half = jju_half + idxu_shift;
const int idxu = jju + idxu_shift;
my_sna.ylist_pack_re(iatom_mod, idxu, ielem, iatom_div) = 0.;
my_sna.ylist_pack_im(iatom_mod, idxu, ielem, iatom_div) = 0.;
auto utot_re = my_sna.ulisttot_re(idxu_half, ielem, iatom);
auto utot_im = my_sna.ulisttot_im(idxu_half, ielem, iatom);
// Store
my_sna.ulisttot_pack(iatom_mod, idxu, ielem, iatom_div) = { utot_re, utot_im };
// Also zero yi
my_sna.ylist_pack_re(iatom_mod, idxu_half, ielem, iatom_div) = 0.;
my_sna.ylist_pack_im(iatom_mod, idxu_half, ielem, iatom_div) = 0.;
// Symmetric term
const int sign_factor = (((ma+mb)%2==0)?1:-1);
const int idxu_flip = jju + (j + 1 - mb) * (j + 1) - (ma + 1);
if (sign_factor == 1) {
utot_im = -utot_im;
} else {
utot_re = -utot_re;
}
my_sna.ulisttot_pack(iatom_mod, idxu_flip, ielem, iatom_div) = { utot_re, utot_im };
// No need to zero symmetrized ylist
}
}
}
}
template<class DeviceType>
@ -742,20 +777,20 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeYi,const int iato
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPTransformYi,const int iatom_mod, const int idxu, const int iatom_div) const {
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPTransformYi,const int iatom_mod, const int idxu_half, const int iatom_div) const {
SNAKokkos<DeviceType> my_sna = snaKK;
const int iatom = iatom_mod + iatom_div * 32;
if (iatom >= chunk_size) return;
if (idxu >= my_sna.idxu_max) return;
if (idxu_half >= my_sna.idxu_half_max) return;
int elem_count = chemflag ? nelements : 1;
for (int ielem = 0; ielem < elem_count; ielem++) {
const auto y_re = my_sna.ylist_pack_re(iatom_mod, idxu, ielem, iatom_div);
const auto y_im = my_sna.ylist_pack_im(iatom_mod, idxu, ielem, iatom_div);
const auto y_re = my_sna.ylist_pack_re(iatom_mod, idxu_half, ielem, iatom_div);
const auto y_im = my_sna.ylist_pack_im(iatom_mod, idxu_half, ielem, iatom_div);
my_sna.ylist(idxu, ielem, iatom) = { y_re, y_im };
my_sna.ylist(idxu_half, ielem, iatom) = { y_re, y_im };
}
}
@ -904,22 +939,52 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeUiCPU,const typen
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPZeroYiCPU,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPZeroYiCPU>::member_type& team) const {
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPTransformUiCPU, const int j, const int iatom) const {
SNAKokkos<DeviceType> my_sna = snaKK;
// Extract the quantum number
const int idx = team.team_rank() + team.team_size() * (team.league_rank() % ((my_sna.idxu_max+team.team_size()-1)/team.team_size()));
if (idx >= my_sna.idxu_max) return;
if (iatom >= chunk_size) return;
// Extract the atomic index
const int ii = team.league_rank() / ((my_sna.idxu_max+team.team_size()-1)/team.team_size());
if (ii >= chunk_size) return;
if (j > twojmax) return;
if (chemflag)
for(int ielem = 0; ielem < nelements; ielem++)
my_sna.zero_yi_cpu(idx,ii,ielem);
else
my_sna.zero_yi_cpu(idx,ii,0);
int elem_count = chemflag ? nelements : 1;
// De-symmetrize ulisttot
for (int ielem = 0; ielem < elem_count; ielem++) {
const int jju_half = my_sna.idxu_half_block(j);
const int jju = my_sna.idxu_block(j);
for (int mb = 0; 2*mb <= j; mb++) {
for (int ma = 0; ma <= j; ma++) {
// Extract top half
const int idxu_shift = mb * (j + 1) + ma;
const int idxu_half = jju_half + idxu_shift;
const int idxu = jju + idxu_shift;
// Load ulist
auto utot = my_sna.ulisttot(idxu_half, ielem, iatom);
// Store
my_sna.ulisttot_full(idxu, ielem, iatom) = utot;
// Zero Yi
my_sna.ylist(idxu_half, ielem, iatom) = {0., 0.};
// Symmetric term
const int sign_factor = (((ma+mb)%2==0)?1:-1);
const int idxu_flip = jju + (j + 1 - mb) * (j + 1) - (ma + 1);
if (sign_factor == 1) {
utot.im = -utot.im;
} else {
utot.re = -utot.re;
}
my_sna.ulisttot_full(idxu_flip, ielem, iatom) = utot;
}
}
}
}
template<class DeviceType>

21
src/KOKKOS/sna_kokkos.h
View File

@ -55,6 +55,8 @@ public:
typedef Kokkos::View<SNAcomplex**[3], DeviceType> t_sna_3c3;
typedef Kokkos::View<SNAcomplex*****, DeviceType> t_sna_5c;
typedef Kokkos::View<CayleyKleinPack**, DeviceType> t_sna_2ckp;
inline
SNAKokkos() {};
KOKKOS_INLINE_FUNCTION
@ -78,6 +80,9 @@ inline
// functions for bispectrum coefficients, GPU only
KOKKOS_INLINE_FUNCTION
void compute_cayley_klein(const int&, const int&, const double&, const double&,
const double&, const double&, const double&);
KOKKOS_INLINE_FUNCTION
void pre_ui(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team,const int&,const int&); // ForceSNAP
KOKKOS_INLINE_FUNCTION
void compute_ui(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int); // ForceSNAP
@ -97,8 +102,6 @@ inline
KOKKOS_INLINE_FUNCTION
void compute_zi_cpu(const int&); // ForceSNAP
KOKKOS_INLINE_FUNCTION
void zero_yi_cpu(const int&,const int&,const int&); // ForceSNAP
KOKKOS_INLINE_FUNCTION
void compute_yi_cpu(int,
const Kokkos::View<F_FLOAT**, DeviceType> &beta); // ForceSNAP
KOKKOS_INLINE_FUNCTION
@ -117,6 +120,8 @@ inline
double compute_sfac(double, double); // add_uarraytot, compute_duarray
KOKKOS_INLINE_FUNCTION
double compute_dsfac(double, double); // compute_duarray
KOKKOS_INLINE_FUNCTION
void compute_s_dsfac(const double, const double, double&, double&); // compute_cayley_klein
// efficient complex FMA
// efficient caxpy (i.e., y += a x)
@ -140,6 +145,9 @@ inline
//per sna class instance for OMP use
// Alternative to rij, wj, rcutij...
// just calculate everything up front
t_sna_2ckp cayleyklein;
// Per InFlight Particle
t_sna_3d rij;
@ -156,6 +164,7 @@ inline
t_sna_3d_ll blist;
t_sna_3c_ll ulisttot;
t_sna_3c_ll ulisttot_full; // un-folded ulisttot, cpu only
t_sna_3c_ll zlist;
t_sna_3c_ll ulist;
@ -173,7 +182,7 @@ inline
t_sna_4d_ll ylist_pack_re; // split real,
t_sna_4d_ll ylist_pack_im; // imag AoSoA layout
int idxcg_max, idxu_max, idxz_max, idxb_max;
int idxcg_max, idxu_max, idxu_half_max, idxu_cache_max, idxz_max, idxb_max;
// Chem snap counts
int nelements;
@ -188,7 +197,13 @@ private:
Kokkos::View<int*[10], DeviceType> idxz;
Kokkos::View<int*[3], DeviceType> idxb;
Kokkos::View<int***, DeviceType> idxcg_block;
public:
Kokkos::View<int*, DeviceType> idxu_block;
Kokkos::View<int*, DeviceType> idxu_half_block;
Kokkos::View<int*, DeviceType> idxu_cache_block;
private:
Kokkos::View<int***, DeviceType> idxz_block;
Kokkos::View<int***, DeviceType> idxb_block;

624
src/KOKKOS/sna_kokkos_impl.h
View File

@ -120,6 +120,36 @@ void SNAKokkos<DeviceType>::build_indexlist()
idxu_max = idxu_count;
Kokkos::deep_copy(idxu_block,h_idxu_block);
// index list for half uarray
idxu_half_block = Kokkos::View<int*, DeviceType>("SNAKokkos::idxu_half_block",jdim);
auto h_idxu_half_block = Kokkos::create_mirror_view(idxu_half_block);
int idxu_half_count = 0;
for(int j = 0; j <= twojmax; j++) {
h_idxu_half_block[j] = idxu_half_count;
for(int mb = 0; 2*mb <= j; mb++)
for(int ma = 0; ma <= j; ma++)
idxu_half_count++;
}
idxu_half_max = idxu_half_count;
Kokkos::deep_copy(idxu_half_block, h_idxu_half_block);
// index list for "cache" uarray
// this is the GPU scratch memory requirements
// applied the CPU structures
idxu_cache_block = Kokkos::View<int*, DeviceType>("SNAKokkos::idxu_cache_block",jdim);
auto h_idxu_cache_block = Kokkos::create_mirror_view(idxu_cache_block);
int idxu_cache_count = 0;
for(int j = 0; j <= twojmax; j++) {
h_idxu_cache_block[j] = idxu_cache_count;
for(int mb = 0; mb < ((j+3)/2); mb++)
for (int ma = 0; ma <= j; ma++)
idxu_cache_count++;
}
idxu_cache_max = idxu_cache_count;
Kokkos::deep_copy(idxu_cache_block, h_idxu_cache_block);
// index list for beta and B
int idxb_count = 0;
@ -201,15 +231,17 @@ void SNAKokkos<DeviceType>::build_indexlist()
h_idxz(idxz_count,8) = MIN(j1, (2 * mb - j + j2 + j1) / 2) - h_idxz(idxz_count,5) + 1;
// apply to z(j1,j2,j,ma,mb) to unique element of y(j)
const int jju = h_idxu_block[j] + (j+1)*mb + ma;
h_idxz(idxz_count,9) = jju;
// ylist is "compressed" via symmetry in its
// contraction with dulist
const int jju_half = h_idxu_half_block[j] + (j+1)*mb + ma;
h_idxz(idxz_count,9) = jju_half;
idxz_count++;
}
}
Kokkos::deep_copy(idxz,h_idxz);
Kokkos::deep_copy(idxz_block,h_idxz_block);
}
/* ---------------------------------------------------------------------- */
@ -230,44 +262,47 @@ void SNAKokkos<DeviceType>::grow_rij(int newnatom, int newnmax)
natom = newnatom;
nmax = newnmax;
rij = t_sna_3d("sna:rij",natom,nmax,3);
inside = t_sna_2i("sna:inside",natom,nmax);
wj = t_sna_2d("sna:wj",natom,nmax);
rcutij = t_sna_2d("sna:rcutij",natom,nmax);
element = t_sna_2i("sna:rcutij",natom,nmax);
dedr = t_sna_3d("sna:dedr",natom,nmax,3);
inside = t_sna_2i(Kokkos::ViewAllocateWithoutInitializing("sna:inside"),natom,nmax);
element = t_sna_2i(Kokkos::ViewAllocateWithoutInitializing("sna:rcutij"),natom,nmax);
dedr = t_sna_3d(Kokkos::ViewAllocateWithoutInitializing("sna:dedr"),natom,nmax,3);
#ifdef LMP_KOKKOS_GPU
if (std::is_same<DeviceType,Kokkos::Cuda>::value) {
// dummy allocation
ulisttot = t_sna_3c_ll("sna:ulisttot",1,1,1);
ulisttot_re = t_sna_3d_ll("sna:ulisttot_re",idxu_max,nelements,natom);
ulisttot_im = t_sna_3d_ll("sna:ulisttot_im",idxu_max,nelements,natom);
ulisttot_pack = t_sna_4c_ll("sna:ulisttot_pack",32,idxu_max,nelements,(natom+32-1)/32);
ulist = t_sna_3c_ll("sna:ulist",1,1,1);
zlist = t_sna_3c_ll("sna:zlist",1,1,1);
zlist_pack = t_sna_4c_ll("sna:zlist_pack",32,idxz_max,ndoubles,(natom+32-1)/32);
blist = t_sna_3d_ll("sna:blist",idxb_max,ntriples,natom);
blist_pack = t_sna_4d_ll("sna:blist_pack",32,idxb_max,ntriples,(natom+32-1)/32);
ylist = t_sna_3c_ll("sna:ylist",idxu_max,nelements,natom);
ylist_pack_re = t_sna_4d_ll("sna:ylist_pack_re",32,idxu_max,nelements,(natom+32-1)/32);
ylist_pack_im = t_sna_4d_ll("sna:ylist_pack_im",32,idxu_max,nelements,(natom+32-1)/32);
dulist = t_sna_4c3_ll("sna:dulist",1,1,1);
cayleyklein = t_sna_2ckp(Kokkos::ViewAllocateWithoutInitializing("sna:cayleyklein"), natom, nmax);
ulisttot = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot"),1,1,1); // dummy allocation
ulisttot_full = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot"),1,1,1);
ulisttot_re = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_re"),idxu_half_max,nelements,natom);
ulisttot_im = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_im"),idxu_half_max,nelements,natom);
ulisttot_pack = t_sna_4c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_pack"),32,idxu_max,nelements,(natom+32-1)/32);
ulist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulist"),1,1,1);
zlist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:zlist"),1,1,1);
zlist_pack = t_sna_4c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:zlist_pack"),32,idxz_max,ndoubles,(natom+32-1)/32);
blist = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:blist"),idxb_max,ntriples,natom);
blist_pack = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:blist_pack"),32,idxb_max,ntriples,(natom+32-1)/32);
ylist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist"),idxu_half_max,nelements,natom);
ylist_pack_re = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist_pack_re"),32,idxu_half_max,nelements,(natom+32-1)/32);
ylist_pack_im = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist_pack_im"),32,idxu_half_max,nelements,(natom+32-1)/32);
dulist = t_sna_4c3_ll(Kokkos::ViewAllocateWithoutInitializing("sna:dulist"),1,1,1);
} else {
#endif
ulisttot = t_sna_3c_ll("sna:ulisttot",idxu_max,nelements,natom);
ulisttot_re = t_sna_3d_ll("sna:ulisttot_re",1,1,1);
ulisttot_im = t_sna_3d_ll("sna:ulisttot_im",1,1,1);
ulisttot_pack = t_sna_4c_ll("sna:ulisttot_pack",1,1,1,1);
ulist = t_sna_3c_ll("sna:ulist",idxu_max,natom,nmax);
zlist = t_sna_3c_ll("sna:zlist",idxz_max,ndoubles,natom);
zlist_pack = t_sna_4c_ll("sna:zlist_pack",1,1,1,1);
blist = t_sna_3d_ll("sna:blist",idxb_max,ntriples,natom);
blist_pack = t_sna_4d_ll("sna:blist_pack",1,1,1,1);
ylist = t_sna_3c_ll("sna:ylist",idxu_max,nelements,natom);
ylist_pack_re = t_sna_4d_ll("sna:ylist_pack_re",1,1,1,1);
ylist_pack_im = t_sna_4d_ll("sna:ylist_pack_im",1,1,1,1);
dulist = t_sna_4c3_ll("sna:dulist",idxu_max,natom,nmax);
rij = t_sna_3d(Kokkos::ViewAllocateWithoutInitializing("sna:rij"),natom,nmax,3);
wj = t_sna_2d(Kokkos::ViewAllocateWithoutInitializing("sna:wj"),natom,nmax);
rcutij = t_sna_2d(Kokkos::ViewAllocateWithoutInitializing("sna:rcutij"),natom,nmax);
ulisttot = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot"),idxu_half_max,nelements,natom);
ulisttot_full = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_full"),idxu_max,nelements,natom);
ulisttot_re = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_re"),1,1,1);
ulisttot_im = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_im"),1,1,1);
ulisttot_pack = t_sna_4c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulisttot_pack"),1,1,1,1);
ulist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ulist"),idxu_cache_max,natom,nmax);
zlist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:zlist"),idxz_max,ndoubles,natom);
zlist_pack = t_sna_4c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:zlist_pack"),1,1,1,1);
blist = t_sna_3d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:blist"),idxb_max,ntriples,natom);
blist_pack = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:blist_pack"),1,1,1,1);
ylist = t_sna_3c_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist"),idxu_half_max,nelements,natom);
ylist_pack_re = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist_pack_re"),1,1,1,1);
ylist_pack_im = t_sna_4d_ll(Kokkos::ViewAllocateWithoutInitializing("sna:ylist_pack_im"),1,1,1,1);
dulist = t_sna_4c3_ll(Kokkos::ViewAllocateWithoutInitializing("sna:dulist"),idxu_cache_max,natom,nmax);
#ifdef LMP_KOKKOS_GPU
}
@ -278,20 +313,111 @@ void SNAKokkos<DeviceType>::grow_rij(int newnatom, int newnmax)
* GPU routines
* ----------------------------------------------------------------------*/
/* ----------------------------------------------------------------------
Precompute the Cayley-Klein parameters and the derivatives thereof.
This routine better exploits parallelism than the GPU ComputeUi and
ComputeFusedDeidrj, which are one warp per atom-neighbor pair.
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void SNAKokkos<DeviceType>::compute_cayley_klein(const int& iatom, const int& jnbor, const double& x, const double& y,
const double& z, const double& rcut, const double& wj_local)
{
//const double x = rij(iatom,jnbor,0);
//const double y = rij(iatom,jnbor,1);
//const double z = rij(iatom,jnbor,2);
const double rsq = x * x + y * y + z * z;
const double r = sqrt(rsq);
//const double rcut = rcutij(iatom, jnbor);
const double rscale0 = rfac0 * MY_PI / (rcut - rmin0);
const double theta0 = (r - rmin0) * rscale0;
double sn, cs;
sincos(theta0, &sn, &cs);
const double z0 = r * cs / sn;
const double dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq;
//const double wj_local = wj(iatom, jnbor);
double sfac, dsfac;
compute_s_dsfac(r, rcut, sfac, dsfac);
sfac *= wj_local;
dsfac *= wj_local;
const double rinv = 1.0 / r;
const double ux = x * rinv;
const double uy = y * rinv;
const double uz = z * rinv;
const double r0inv = 1.0 / sqrt(r * r + z0 * z0);
const SNAcomplex a = { z0 * r0inv, -z * r0inv };
const SNAcomplex b = { r0inv * y, -r0inv * x };
const double dr0invdr = -r0inv * r0inv * r0inv * (r + z0 * dz0dr);
const double dr0invx = dr0invdr * ux;
const double dr0invy = dr0invdr * uy;
const double dr0invz = dr0invdr * uz;
const double dz0x = dz0dr * ux;
const double dz0y = dz0dr * uy;
const double dz0z = dz0dr * uz;
const SNAcomplex dax = { dz0x * r0inv + z0 * dr0invx, -z * dr0invx };
const SNAcomplex day = { dz0y * r0inv + z0 * dr0invy, -z * dr0invy };
const SNAcomplex daz = { dz0z * r0inv + z0 * dr0invz, -z * dr0invz - r0inv };
const SNAcomplex dbx = { y * dr0invx, -x * dr0invx - r0inv };
const SNAcomplex dby = { y * dr0invy + r0inv, -x * dr0invy };
const SNAcomplex dbz = { y * dr0invz, -x * dr0invz };
const double dsfacux = dsfac * ux;
const double dsfacuy = dsfac * uy;
const double dsfacuz = dsfac * uz;
CayleyKleinPack ckp{};
ckp.a = a;
ckp.b = b;
ckp.da[0] = dax;
ckp.db[0] = dbx;
ckp.da[1] = day;
ckp.db[1] = dby;
ckp.da[2] = daz;
ckp.db[2] = dbz;
ckp.sfac = sfac;
ckp.dsfacu[0] = dsfacux;
ckp.dsfacu[1] = dsfacuy;
ckp.dsfacu[2] = dsfacuz;
// Yes, this breaks the standard mantra of using SoA
// instead of AoS, but it's net fine because of the
// one warp per atom/neighbor pair for the recursive
// polynomials. There's good L1 reuse, anyway.
cayleyklein(iatom, jnbor) = ckp;
}
/* ----------------------------------------------------------------------
Initialize ulisttot with self-energy terms.
Ulisttot uses a "half" data layout which takes
advantage of the symmetry of the Wigner U matrices.
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void SNAKokkos<DeviceType>::pre_ui(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int& iatom, const int& ielem)
{
for (int jelem = 0; jelem < nelements; jelem++) {
for (int j = 0; j <= twojmax; j++) {
const int jju = idxu_block(j);
const int jju_half = idxu_half_block(j);
// Only diagonal elements get initialized
// for (int m = 0; m < (j+1)*(j+1); m++)
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, (j+1)*(j+1)),
// Top half only: gets symmetrized by TransformUi
// for (int m = 0; m < (j+1)*(j/2+1); m++)
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, (j+1)*(j/2+1)),
[&] (const int m) {
const int jjup = jju + m;
const int jjup = jju_half + m;
// if m is on the "diagonal", initialize it with the self energy.
// Otherwise zero it out
@ -323,40 +449,27 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
// utot(j,ma,mb) += u(r0;j,ma,mb) for all j,ma,mb
// get shared memory offset
const int max_m_tile = (twojmax+1)*(twojmax+1);
const int max_m_tile = (twojmax+1)*(twojmax/2+1);
const int team_rank = team.team_rank();
const int scratch_shift = team_rank * max_m_tile;
// double buffer
// shared memory double buffer
// Each level stores a (j+1) x ((j+3)/2) tile in preparation
// for the computation on the next level. This is the "cached"
// data layout, which also gets used on the CPU.
SNAcomplex* buf1 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
SNAcomplex* buf2 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
const double x = rij(iatom,jnbor,0);
const double y = rij(iatom,jnbor,1);
const double z = rij(iatom,jnbor,2);
// Each warp is accessing the same pack: take advantage of
// L1 reuse
const CayleyKleinPack& ckp = cayleyklein(iatom, jnbor);
const SNAcomplex a = ckp.a;
const SNAcomplex b = ckp.b;
const double sfac = ckp.sfac;
const double wj_local = wj(iatom, jnbor);
const double rcut = rcutij(iatom, jnbor);
const int jelem = element(iatom, jnbor);
const double rsq = x * x + y * y + z * z;
const double r = sqrt(rsq);
const double theta0 = (r - rmin0) * rfac0 * MY_PI / (rcutij(iatom,jnbor) - rmin0);
// theta0 = (r - rmin0) * rscale0;
const double cs = cos(theta0);
const double sn = sin(theta0);
const double z0 = r * cs / sn; // r / tan(theta0)
// Compute cutoff function
const double sfac = compute_sfac(r, rcut) * wj_local;
// compute Cayley-Klein parameters for unit quaternion,
// pack into complex number
const double r0inv = 1.0 / sqrt(r * r + z0 * z0);
const SNAcomplex a = { r0inv * z0, -r0inv * z };
const SNAcomplex b = { r0inv * y, -r0inv * x };
// VMK Section 4.8.2
// All writes go to global memory and shared memory
@ -367,7 +480,8 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
});
for (int j = 1; j <= twojmax; j++) {
const int jju = idxu_block[j];
const int jju = idxu_half_block[j];
// fill in left side of matrix layer from previous layer
@ -377,8 +491,7 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
// for (int mb = 0; 2*mb <= j; mb++)
const int n_mb = j/2+1;
// the last (j / 2) can be avoided due to symmetry
const int total_iters = n_ma * n_mb - (j % 2 == 0 ? (j / 2) : 0);
const int total_iters = n_ma * n_mb;
//for (int m = 0; m < total_iters; m++) {
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, total_iters),
@ -423,25 +536,22 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j,mb-j] = (-1)^(ma-mb)*Conj([u[ma,mb))
// if j is even (-> physical j integer), last element maps to self, skip
//if (!(m == total_iters - 1 && j % 2 == 0)) {
if (m < total_iters - 1 || j % 2 == 1) {
const int sign_factor = (((ma+mb)%2==0)?1:-1);
const int jju_shared_flip = (j+1-mb)*(j+1)-(ma+1);
const int jjup_flip = jju + jju_shared_flip; // jju+(j+1-mb)*(j+1)-(ma+1);
const int sign_factor = (((ma+mb)%2==0)?1:-1);
const int jju_shared_flip = (j+1-mb)*(j+1)-(ma+1);
if (sign_factor == 1) {
u_accum.im = -u_accum.im;
} else {
u_accum.re = -u_accum.re;
}
buf2[jju_shared_flip] = u_accum;
// split re, im to get fully coalesced atomic add
Kokkos::atomic_add(&(ulisttot_re(jjup_flip,jelem,iatom)), sfac * u_accum.re);
Kokkos::atomic_add(&(ulisttot_im(jjup_flip,jelem,iatom)), sfac * u_accum.im);
if (sign_factor == 1) {
u_accum.im = -u_accum.im;
} else {
u_accum.re = -u_accum.re;
}
if (j%2==1 && mb+1==n_mb) {
buf2[jju_shared_flip] = u_accum;
}
// symmetric part of ulisttot is generated in TransformUi
});
// In CUDA backend,
// ThreadVectorRange has a __syncwarp (appropriately masked for
@ -454,8 +564,11 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
}
}
/* ----------------------------------------------------------------------
compute Zi by summing over products of Ui, GPU version
compute Zi by summing over products of Ui,
AoSoA data layout to take advantage of coalescing, avoiding warp
divergence. GPU version
------------------------------------------------------------------------- */
template<class DeviceType>
@ -525,6 +638,8 @@ void SNAKokkos<DeviceType>::compute_zi(const int& iatom_mod, const int& jjz, con
/* ----------------------------------------------------------------------
compute Bi by summing conj(Ui)*Zi
AoSoA data layout to take advantage of coalescing, avoiding warp
divergence.
------------------------------------------------------------------------- */
template<class DeviceType>
@ -615,8 +730,9 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
/* ----------------------------------------------------------------------
compute Yi from Ui without storing Zi, looping over zlist indices,
GPU version
compute Yi from Ui without storing Zi, looping over zlist indices.
AoSoA data layout to take advantage of coalescing, avoiding warp
divergence. GPU version.
------------------------------------------------------------------------- */
template<class DeviceType>
@ -635,7 +751,7 @@ void SNAKokkos<DeviceType>::compute_yi(int iatom_mod, int jjz, int iatom_div,
const int mb2max = idxz(jjz, 6);
const int na = idxz(jjz, 7);
const int nb = idxz(jjz, 8);
const int jju = idxz(jjz, 9);
const int jju_half = idxz(jjz, 9);
const double *cgblock = cglist.data() + idxcg_block(j1,j2,j);
//int mb = (2 * (mb1min+mb2max) - j1 - j2 + j) / 2;
@ -677,8 +793,8 @@ void SNAKokkos<DeviceType>::compute_yi(int iatom_mod, int jjz, int iatom_div,
} // end loop over ib
if (bnorm_flag) {
ztmp_r /= j + 1;
ztmp_i /= j + 1;
ztmp_r /= (j + 1);
ztmp_i /= (j + 1);
}
// apply to z(j1,j2,j,ma,mb) to unique element of y(j)
@ -710,8 +826,8 @@ void SNAKokkos<DeviceType>::compute_yi(int iatom_mod, int jjz, int iatom_div,
betaj *= (j1 + 1) / (j + 1.0);
Kokkos::atomic_add(&(ylist_pack_re(iatom_mod, jju, elem3, iatom_div)), betaj*ztmp_r);
Kokkos::atomic_add(&(ylist_pack_im(iatom_mod, jju, elem3, iatom_div)), betaj*ztmp_i);
Kokkos::atomic_add(&(ylist_pack_re(iatom_mod, jju_half, elem3, iatom_div)), betaj*ztmp_r);
Kokkos::atomic_add(&(ylist_pack_im(iatom_mod, jju_half, elem3, iatom_div)), betaj*ztmp_i);
} // end loop over elem3
} // end loop over elem2
} // end loop over elem1
@ -719,7 +835,7 @@ void SNAKokkos<DeviceType>::compute_yi(int iatom_mod, int jjz, int iatom_div,
/* ----------------------------------------------------------------------
Fused calculation of the derivative of Ui w.r.t. atom j
and of dEidRj. GPU only.
and accumulation into dEidRj. GPU only.
------------------------------------------------------------------------- */
template<class DeviceType>
@ -732,6 +848,9 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
const int scratch_shift = team_rank * max_m_tile;
const int jelem = element(iatom, jnbor);
// See notes on data layouts for shared memory caching
// in `compute_ui`.
// double buffer for ulist
SNAcomplex* ulist_buf1 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
SNAcomplex* ulist_buf2 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
@ -740,50 +859,20 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
SNAcomplex* dulist_buf1 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
SNAcomplex* dulist_buf2 = (SNAcomplex*)team.team_shmem( ).get_shmem(team.team_size()*max_m_tile*sizeof(SNAcomplex), 0) + scratch_shift;
const double x = rij(iatom,jnbor,0);
const double y = rij(iatom,jnbor,1);
const double z = rij(iatom,jnbor,2);
const double rsq = x * x + y * y + z * z;
const double r = sqrt(rsq);
const double rcut = rcutij(iatom, jnbor);
const double rscale0 = rfac0 * MY_PI / (rcut - rmin0);
const double theta0 = (r - rmin0) * rscale0;
const double cs = cos(theta0);
const double sn = sin(theta0);
const double z0 = r * cs / sn;
const double dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq;
const CayleyKleinPack& ckp = cayleyklein(iatom, jnbor);
const double wj_local = wj(iatom, jnbor);
const double sfac = wj_local * compute_sfac(r, rcut);
const double dsfac = wj_local * compute_dsfac(r, rcut);
const double rinv = 1.0 / r;
// extract a single unit vector
const double u = (dir == 0 ? x * rinv : dir == 1 ? y * rinv : z * rinv);
// Compute Cayley-Klein parameters for unit quaternion
const double r0inv = 1.0 / sqrt(r * r + z0 * z0);
const SNAcomplex a = { r0inv * z0, -r0inv * z };
const SNAcomplex b = { r0inv * y, -r0inv * x };
const double dr0invdr = -r0inv * r0inv * r0inv * (r + z0 * dz0dr);
const double dr0inv = dr0invdr * u;
const double dz0 = dz0dr * u;
const SNAcomplex da = { dz0 * r0inv + z0 * dr0inv,
- z * dr0inv + (dir == 2 ? - r0inv : 0.) };
const SNAcomplex db = { y * dr0inv + (dir==1?r0inv:0.),
-x * dr0inv + (dir==0?-r0inv:0.) };
const SNAcomplex a = ckp.a;
const SNAcomplex b = ckp.b;
const SNAcomplex da = ckp.da[dir];
const SNAcomplex db = ckp.db[dir];
const double sfac = ckp.sfac;
const double dsfacu = ckp.dsfacu[dir]; // dsfac * u
// Accumulate the full contribution to dedr on the fly
const double du_prod = dsfac * u; // chain rule
const SNAcomplex y_local = ylist(0, jelem, iatom);
// Symmetry factor of 0.5 b/c 0 element is on diagonal for even j==0
double dedr_full_sum = 0.5 * du_prod * y_local.re;
double dedr_full_sum = 0.5 * dsfacu * y_local.re;
// single has a warp barrier at the end
Kokkos::single(Kokkos::PerThread(team), [=]() {
@ -793,8 +882,7 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
});
for (int j = 1; j <= twojmax; j++) {
int jju = idxu_block[j];
int jjup = idxu_block[j-1];
int jju_half = idxu_half_block[j];
// flatten the loop over ma,mb
@ -812,11 +900,10 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
[&] (const int m, double& sum_tmp) {
// ma fast, mb slow
// Equivalent to `int ma = m % n_ma; int mb = m / n_ma;` IF everything's positive.
const int mb = m / n_ma;
const int ma = m - mb * n_ma;
const int jju_index = jju+m;
const int jju_half_index = jju_half+m;
// Load y_local, apply the symmetry scaling factor
// The "secret" of the shared memory optimization is it eliminates
@ -824,7 +911,7 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
// shared memory and otherwise always writing, making the kernel
// ultimately compute bound. We take advantage of that by adding
// some reads back in.
auto y_local = ylist(jju_index, jelem, iatom);
auto y_local = ylist(jju_half_index, jelem, iatom);
if (j % 2 == 0 && 2*mb == j) {
if (ma == mb) { y_local = 0.5*y_local; }
else if (ma > mb) { y_local = { 0., 0. }; } // can probably avoid this outright
@ -870,7 +957,7 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
dulist_buf2[jju_shared_idx] = du_accum;
// Directly accumulate deidrj into sum_tmp
const SNAcomplex du_prod = ((dsfac * u)*u_accum) + (sfac*du_accum);
const SNAcomplex du_prod = (dsfacu * u_accum) + (sfac * du_accum);
sum_tmp += du_prod.re * y_local.re + du_prod.im * y_local.im;
// copy left side to right side with inversion symmetry VMK 4.4(2)
@ -923,8 +1010,9 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
* ----------------------------------------------------------------------*/
/* ----------------------------------------------------------------------
* compute Ui by summing over neighbors j
* ------------------------------------------------------------------------- */
Ulisttot uses a "half" data layout which takes
advantage of the symmetry of the Wigner U matrices.
* ------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
@ -932,11 +1020,11 @@ void SNAKokkos<DeviceType>::pre_ui_cpu(const typename Kokkos::TeamPolicy<DeviceT
{
for (int jelem = 0; jelem < nelements; jelem++) {
for (int j = 0; j <= twojmax; j++) {
const int jju = idxu_block(j);
const int jju = idxu_half_block(j);
// Only diagonal elements get initialized
// for (int m = 0; m < (j+1)*(j+1); m++)
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, (j+1)*(j+1)),
// for (int m = 0; m < (j+1)*(j/2+1); m++)
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, (j+1)*(j/2+1)),
[&] (const int m) {
const int jjup = jju + m;
@ -956,6 +1044,8 @@ void SNAKokkos<DeviceType>::pre_ui_cpu(const typename Kokkos::TeamPolicy<DeviceT
/* ----------------------------------------------------------------------
compute Ui by summing over bispectrum components. CPU only.
See comments above compute_uarray_cpu and add_uarraytot for
data layout comments.
------------------------------------------------------------------------- */
template<class DeviceType>
@ -1026,10 +1116,10 @@ void SNAKokkos<DeviceType>::compute_zi_cpu(const int& iter)
int ma2 = ma2max;
int icga = ma1min * (j2 + 1) + ma2max;
for(int ia = 0; ia < na; ia++) {
suma1_r += cgblock[icga] * (ulisttot(jju1+ma1, elem1, iatom).re * ulisttot(jju2+ma2, elem2, iatom).re -
ulisttot(jju1+ma1, elem1, iatom).im * ulisttot(jju2+ma2, elem2, iatom).im);
suma1_i += cgblock[icga] * (ulisttot(jju1+ma1, elem1, iatom).re * ulisttot(jju2+ma2, elem2, iatom).im +
ulisttot(jju1+ma1, elem1, iatom).im * ulisttot(jju2+ma2, elem2, iatom).re);
suma1_r += cgblock[icga] * (ulisttot_full(jju1+ma1, elem1, iatom).re * ulisttot_full(jju2+ma2, elem2, iatom).re -
ulisttot_full(jju1+ma1, elem1, iatom).im * ulisttot_full(jju2+ma2, elem2, iatom).im);
suma1_i += cgblock[icga] * (ulisttot_full(jju1+ma1, elem1, iatom).re * ulisttot_full(jju2+ma2, elem2, iatom).im +
ulisttot_full(jju1+ma1, elem1, iatom).im * ulisttot_full(jju2+ma2, elem2, iatom).re);
ma1++;
ma2--;
icga += j2;
@ -1098,8 +1188,8 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
const int jjz_index = jjz + mb * (j + 1) + ma;
if (2*mb == j) return;
sum +=
ulisttot(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im;
ulisttot_full(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot_full(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im;
},sumzu_temp); // end loop over ma, mb
sumzu += sumzu_temp;
@ -1113,8 +1203,8 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
const int jju_index = jju+(mb-1)*(j+1)+(j+1)+ma;
const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma;
sum +=
ulisttot(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im;
ulisttot_full(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot_full(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im;
},sumzu_temp); // end loop over ma
sumzu += sumzu_temp;
@ -1122,8 +1212,8 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
const int jju_index = jju+(mb-1)*(j+1)+(j+1)+ma;
const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma;
sumzu += 0.5*
(ulisttot(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im);
(ulisttot_full(jju_index, elem3, iatom).re * zlist(jjz_index, jalloy, iatom).re +
ulisttot_full(jju_index, elem3, iatom).im * zlist(jjz_index, jalloy, iatom).im);
} // end if jeven
Kokkos::single(Kokkos::PerThread(team), [&] () {
@ -1154,17 +1244,6 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
}
/* ----------------------------------------------------------------------
compute Yi from Ui without storing Zi, looping over zlist indices
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void SNAKokkos<DeviceType>::zero_yi_cpu(const int& idx, const int& iatom, const int& ielem)
{
ylist(idx,ielem,iatom) = {0.0, 0.0};
}
/* ----------------------------------------------------------------------
compute Yi from Ui without storing Zi, looping over zlist indices,
CPU version
@ -1188,7 +1267,7 @@ void SNAKokkos<DeviceType>::compute_yi_cpu(int iter,
const int mb2max = idxz(jjz, 6);
const int na = idxz(jjz, 7);
const int nb = idxz(jjz, 8);
const int jju = idxz(jjz, 9);
const int jju_half = idxz(jjz, 9);
const double *cgblock = cglist.data() + idxcg_block(j1,j2,j);
//int mb = (2 * (mb1min+mb2max) - j1 - j2 + j) / 2;
@ -1214,10 +1293,10 @@ void SNAKokkos<DeviceType>::compute_yi_cpu(int iter,
int icga = ma1min*(j2+1) + ma2max;
for (int ia = 0; ia < na; ia++) {
suma1_r += cgblock[icga] * (ulisttot(jju1+ma1, elem1, iatom).re * ulisttot(jju2+ma2, elem2, iatom).re -
ulisttot(jju1+ma1, elem1, iatom).im * ulisttot(jju2+ma2, elem2, iatom).im);
suma1_i += cgblock[icga] * (ulisttot(jju1+ma1, elem1, iatom).re * ulisttot(jju2+ma2, elem2, iatom).im +
ulisttot(jju1+ma1, elem1, iatom).im * ulisttot(jju2+ma2, elem2, iatom).re);
suma1_r += cgblock[icga] * (ulisttot_full(jju1+ma1, elem1, iatom).re * ulisttot_full(jju2+ma2, elem2, iatom).re -
ulisttot_full(jju1+ma1, elem1, iatom).im * ulisttot_full(jju2+ma2, elem2, iatom).im);
suma1_i += cgblock[icga] * (ulisttot_full(jju1+ma1, elem1, iatom).re * ulisttot_full(jju2+ma2, elem2, iatom).im +
ulisttot_full(jju1+ma1, elem1, iatom).im * ulisttot_full(jju2+ma2, elem2, iatom).re);
ma1++;
ma2--;
icga += j2;
@ -1264,8 +1343,8 @@ void SNAKokkos<DeviceType>::compute_yi_cpu(int iter,
if (!bnorm_flag && j1 > j)
betaj *= (j1 + 1) / (j + 1.0);
Kokkos::atomic_add(&(ylist(jju, elem3, iatom).re), betaj*ztmp_r);
Kokkos::atomic_add(&(ylist(jju, elem3, iatom).im), betaj*ztmp_i);
Kokkos::atomic_add(&(ylist(jju_half, elem3, iatom).re), betaj*ztmp_r);
Kokkos::atomic_add(&(ylist(jju_half, elem3, iatom).im), betaj*ztmp_i);
} // end loop over elem3
} // end loop over elem2
} // end loop over elem1
@ -1274,6 +1353,8 @@ void SNAKokkos<DeviceType>::compute_yi_cpu(int iter,
/* ----------------------------------------------------------------------
calculate derivative of Ui w.r.t. atom j
see comments above compute_duarray_cpu for comments on the
data layout
------------------------------------------------------------------------- */
template<class DeviceType>
@ -1301,6 +1382,10 @@ void SNAKokkos<DeviceType>::compute_duidrj_cpu(const typename Kokkos::TeamPolicy
/* ----------------------------------------------------------------------
compute dEidRj, CPU path only.
dulist takes advantage of a `cached` data layout, similar to the
shared memory layout for the GPU routines, which is efficient for
compressing the calculation in compute_duarray_cpu. That said,
dulist only uses the "half" data layout part of that structure.
------------------------------------------------------------------------- */
@ -1314,14 +1399,18 @@ void SNAKokkos<DeviceType>::compute_deidrj_cpu(const typename Kokkos::TeamPolicy
//for(int j = 0; j <= twojmax; j++) {
Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(team,twojmax+1),
[&] (const int& j, t_scalar3<double>& sum_tmp) {
int jju = idxu_block[j];
int jju_half = idxu_half_block[j];
int jju_cache = idxu_cache_block[j];
for(int mb = 0; 2*mb < j; mb++)
for(int ma = 0; ma <= j; ma++) {
sum_tmp.x += dulist(jju,iatom,jnbor,0).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,0).im * ylist(jju,jelem,iatom).im;
sum_tmp.y += dulist(jju,iatom,jnbor,1).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,1).im * ylist(jju,jelem,iatom).im;
sum_tmp.z += dulist(jju,iatom,jnbor,2).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,2).im * ylist(jju,jelem,iatom).im;
jju++;
sum_tmp.x += dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im;
sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im;
sum_tmp.z += dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im;
jju_half++; jju_cache++;
} //end loop over ma mb
// For j even, handle middle column
@ -1330,16 +1419,22 @@ void SNAKokkos<DeviceType>::compute_deidrj_cpu(const typename Kokkos::TeamPolicy
int mb = j/2;
for(int ma = 0; ma < mb; ma++) {
sum_tmp.x += dulist(jju,iatom,jnbor,0).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,0).im * ylist(jju,jelem,iatom).im;
sum_tmp.y += dulist(jju,iatom,jnbor,1).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,1).im * ylist(jju,jelem,iatom).im;
sum_tmp.z += dulist(jju,iatom,jnbor,2).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,2).im * ylist(jju,jelem,iatom).im;
jju++;
sum_tmp.x += dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im;
sum_tmp.y += dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im;
sum_tmp.z += dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im;
jju_half++; jju_cache++;
}
//int ma = mb;
sum_tmp.x += (dulist(jju,iatom,jnbor,0).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,0).im * ylist(jju,jelem,iatom).im)*0.5;
sum_tmp.y += (dulist(jju,iatom,jnbor,1).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,1).im * ylist(jju,jelem,iatom).im)*0.5;
sum_tmp.z += (dulist(jju,iatom,jnbor,2).re * ylist(jju,jelem,iatom).re + dulist(jju,iatom,jnbor,2).im * ylist(jju,jelem,iatom).im)*0.5;
sum_tmp.x += (dulist(jju_cache,iatom,jnbor,0).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,0).im * ylist(jju_half,jelem,iatom).im)*0.5;
sum_tmp.y += (dulist(jju_cache,iatom,jnbor,1).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,1).im * ylist(jju_half,jelem,iatom).im)*0.5;
sum_tmp.z += (dulist(jju_cache,iatom,jnbor,2).re * ylist(jju_half,jelem,iatom).re +
dulist(jju_cache,iatom,jnbor,2).im * ylist(jju_half,jelem,iatom).im)*0.5;
} // end if jeven
},final_sum); // end loop over j
@ -1355,6 +1450,10 @@ void SNAKokkos<DeviceType>::compute_deidrj_cpu(const typename Kokkos::TeamPolicy
/* ----------------------------------------------------------------------
add Wigner U-functions for one neighbor to the total
ulist is in a "cached" data layout, which is a compressed layout
which still keeps the recursive calculation simple. On the other hand
`ulisttot` uses a "half" data layout, which fully takes advantage
of the symmetry of the Wigner U matrices.
------------------------------------------------------------------------- */
template<class DeviceType>
@ -1364,15 +1463,27 @@ void SNAKokkos<DeviceType>::add_uarraytot(const typename Kokkos::TeamPolicy<Devi
{
const double sfac = compute_sfac(r, rcut) * wj;
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,ulisttot.extent(0)),
[&] (const int& i) {
Kokkos::atomic_add(&(ulisttot(i,jelem,iatom).re), sfac * ulist(i,iatom,jnbor).re);
Kokkos::atomic_add(&(ulisttot(i,jelem,iatom).im), sfac * ulist(i,iatom,jnbor).im);
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,twojmax+1),
[&] (const int& j) {
int jju_half = idxu_half_block[j]; // index into ulisttot
int jju_cache = idxu_cache_block[j]; // index into ulist
int count = 0;
for (int mb = 0; 2*mb <= j; mb++) {
for (int ma = 0; ma <= j; ma++) {
Kokkos::atomic_add(&(ulisttot(jju_half+count, jelem, iatom).re), sfac * ulist(jju_cache+count, iatom, jnbor).re);
Kokkos::atomic_add(&(ulisttot(jju_half+count, jelem, iatom).im), sfac * ulist(jju_cache+count, iatom, jnbor).im);
count++;
}
}
});
}
/* ----------------------------------------------------------------------
compute Wigner U-functions for one neighbor
compute Wigner U-functions for one neighbor.
`ulisttot` uses a "cached" data layout, matching the amount of
information stored between layers via scratch memory on the GPU path
------------------------------------------------------------------------- */
template<class DeviceType>
@ -1399,8 +1510,8 @@ void SNAKokkos<DeviceType>::compute_uarray_cpu(const typename Kokkos::TeamPolicy
ulist(0,iatom,jnbor).im = 0.0;
for (int j = 1; j <= twojmax; j++) {
int jju = idxu_block[j];
int jjup = idxu_block[j-1];
const int jju = idxu_cache_block[j];
const int jjup = idxu_cache_block[j-1];
// fill in left side of matrix layer from previous layer
@ -1434,37 +1545,37 @@ void SNAKokkos<DeviceType>::compute_uarray_cpu(const typename Kokkos::TeamPolicy
(b_r * ulist(jjup_index,iatom,jnbor).im -
b_i * ulist(jjup_index,iatom,jnbor).re);
}
});
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j,mb-j] = (-1)^(ma-mb)*Conj([u[ma,mb))
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j,mb-j] = (-1)^(ma-mb)*Conj([u[ma,mb))
jju = idxu_block[j];
jjup = jju+(j+1)*(j+1)-1;
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,(j+2)/2),
[&] (const int& mb) {
// for (int mb = 0; 2*mb <= j; mb++) {
int mbpar = (mb)%2==0?1:-1;
int mapar = mbpar;
for (int ma = 0; ma <= j; ma++) {
const int jju_index = jju+mb*(j+1)+ma;
const int jjup_index = jjup-mb*(j+1)-ma;
if (mapar == 1) {
ulist(jjup_index,iatom,jnbor).re = ulist(jju_index,iatom,jnbor).re;
ulist(jjup_index,iatom,jnbor).im = -ulist(jju_index,iatom,jnbor).im;
} else {
ulist(jjup_index,iatom,jnbor).re = -ulist(jju_index,iatom,jnbor).re;
ulist(jjup_index,iatom,jnbor).im = ulist(jju_index,iatom,jnbor).im;
// Only need to add one symmetrized row for convenience
// Symmetry gets "unfolded" in accumulating ulisttot
if (j%2==1 && mb==(j/2)) {
const int mbpar = (mb)%2==0?1:-1;
int mapar = mbpar;
for (int ma = 0; ma <= j; ma++) {
const int jju_index = jju + mb*(j+1) + ma;
const int jjup_index = jju + (j+1-mb)*(j+1)-(ma+1);
if (mapar == 1) {
ulist(jjup_index,iatom,jnbor).re = ulist(jju_index,iatom,jnbor).re;
ulist(jjup_index,iatom,jnbor).im = -ulist(jju_index,iatom,jnbor).im;
} else {
ulist(jjup_index,iatom,jnbor).re = -ulist(jju_index,iatom,jnbor).re;
ulist(jjup_index,iatom,jnbor).im = ulist(jju_index,iatom,jnbor).im;
}
mapar = -mapar;
}
mapar = -mapar;
}
});
}
}
/* ----------------------------------------------------------------------
compute derivatives of Wigner U-functions for one neighbor
see comments in compute_uarray_cpu()
Uses same cached data layout of ulist
------------------------------------------------------------------------- */
template<class DeviceType>
@ -1524,8 +1635,8 @@ double r0inv;
dulist(0,iatom,jnbor,2).im = 0.0;
for (int j = 1; j <= twojmax; j++) {
int jju = idxu_block[j];
int jjup = idxu_block[j-1];
int jju = idxu_cache_block[j];
int jjup = idxu_cache_block[j-1];
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,(j+2)/2),
[&] (const int& mb) {
//for (int mb = 0; 2*mb <= j; mb++) {
@ -1568,33 +1679,32 @@ double r0inv;
b_i * dulist(jjup_index,iatom,jnbor,k).re);
}
}
});
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])
// Only need to add one symmetrized row for convenience
// Symmetry gets "unfolded" during the dedr accumulation
jju = idxu_block[j];
jjup = jju+(j+1)*(j+1)-1;
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,(j+2)/2),
[&] (const int& mb) {
// for (int mb = 0; 2*mb <= j; mb++) {
int mbpar = (mb)%2==0?1:-1;
int mapar = mbpar;
for (int ma = 0; ma <= j; ma++) {
const int jju_index = jju+mb*(j+1)+ma;
const int jjup_index = jjup-mb*(j+1)-ma;
if (mapar == 1) {
for (int k = 0; k < 3; k++) {
dulist(jjup_index,iatom,jnbor,k).re = dulist(jju_index,iatom,jnbor,k).re;
dulist(jjup_index,iatom,jnbor,k).im = -dulist(jju_index,iatom,jnbor,k).im;
}
} else {
for (int k = 0; k < 3; k++) {
dulist(jjup_index,iatom,jnbor,k).re = -dulist(jju_index,iatom,jnbor,k).re;
dulist(jjup_index,iatom,jnbor,k).im = dulist(jju_index,iatom,jnbor,k).im;
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])
if (j%2==1 && mb==(j/2)) {
const int mbpar = (mb)%2==0?1:-1;
int mapar = mbpar;
for (int ma = 0; ma <= j; ma++) {
const int jju_index = jju+mb*(j+1)+ma;
const int jjup_index = jju+(mb+2)*(j+1)-(ma+1);
if (mapar == 1) {
for (int k = 0; k < 3; k++) {
dulist(jjup_index,iatom,jnbor,k).re = dulist(jju_index,iatom,jnbor,k).re;
dulist(jjup_index,iatom,jnbor,k).im = -dulist(jju_index,iatom,jnbor,k).im;
}
} else {
for (int k = 0; k < 3; k++) {
dulist(jjup_index,iatom,jnbor,k).re = -dulist(jju_index,iatom,jnbor,k).re;
dulist(jjup_index,iatom,jnbor,k).im = dulist(jju_index,iatom,jnbor,k).im;
}
}
mapar = -mapar;
}
mapar = -mapar;
}
});
}
@ -1605,8 +1715,11 @@ double r0inv;
sfac *= wj;
dsfac *= wj;
// Even though we fill out a full "cached" data layout above,
// we only need the "half" data for the accumulation into dedr.
// Thus we skip updating any unnecessary data.
for (int j = 0; j <= twojmax; j++) {
int jju = idxu_block[j];
int jju = idxu_cache_block[j];
for (int mb = 0; 2*mb <= j; mb++)
for (int ma = 0; ma <= j; ma++) {
dulist(jju,iatom,jnbor,0).re = dsfac * ulist(jju,iatom,jnbor).re * ux +
@ -1980,6 +2093,24 @@ double SNAKokkos<DeviceType>::compute_dsfac(double r, double rcut)
return 0.0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void SNAKokkos<DeviceType>::compute_s_dsfac(const double r, const double rcut, double& sfac, double& dsfac) {
if (switch_flag == 0) { sfac = 0.; dsfac = 0.; }
else if (switch_flag == 1) {
if (r <= rmin0) { sfac = 1.0; dsfac = 0.0; }
else if (r > rcut) { sfac = 0.; dsfac = 0.; }
else {
const double rcutfac = MY_PI / (rcut - rmin0);
double sn, cs;
sincos((r - rmin0) * rcutfac, &sn, &cs);
sfac = 0.5 * (cs + 1.0);
dsfac = -0.5 * sn * rcutfac;
}
} else { sfac = 0.; dsfac = 0.; }
}
/* ---------------------------------------------------------------------- */
// efficient complex FMA (i.e., y += a x)
@ -2031,30 +2162,31 @@ double SNAKokkos<DeviceType>::memory_usage()
#ifdef LMP_KOKKOS_GPU
if (std::is_same<DeviceType,Kokkos::Cuda>::value) {
int natom_pad = ((natom + 32 - 1) / 32) * 32; // for AoSoA layouts
auto natom_pad = (natom+32-1)/32;
bytes += natom * idxu_max * nelements * sizeof(double); // ulisttot_re
bytes += natom * idxu_max * nelements * sizeof(double); // ulisttot_im
bytes += natom * idxu_half_max * nelements * sizeof(double); // ulisttot_re
bytes += natom * idxu_half_max * nelements * sizeof(double); // ulisttot_im
bytes += natom_pad * idxu_max * nelements * sizeof(double) * 2; // ulisttot_pack
bytes += natom_pad * idxz_max * ndoubles * sizeof(double) * 2; // zlist_pack
bytes += natom_pad * idxb_max * ntriples * sizeof(double); // blist_pack
bytes += natom_pad * idxu_max * nelements * sizeof(double); // ylist_pack_re
bytes += natom_pad * idxu_max * nelements * sizeof(double); // ylist_pack_im
bytes += natom * idxu_max * nelements * sizeof(double) * 2; // ylist
bytes += natom_pad * idxu_half_max * nelements * sizeof(double); // ylist_pack_re
bytes += natom_pad * idxu_half_max * nelements * sizeof(double); // ylist_pack_im
bytes += natom * idxu_half_max * nelements * sizeof(double) * 2; // ylist
} else {
#endif
bytes += natom * nmax * idxu_max * sizeof(double) * 2; // ulist
bytes += natom * idxu_max * nelements * sizeof(double) * 2; // ulisttot
bytes += natom * nmax * idxu_cache_max * sizeof(double) * 2; // ulist
bytes += natom * idxu_half_max * nelements * sizeof(double) * 2; // ulisttot
bytes += natom * idxz_max * ndoubles * sizeof(double) * 2; // zlist
bytes += natom * idxb_max * ntriples * sizeof(double); // blist
bytes += natom * idxu_max * nelements * sizeof(double) * 2; // ylist
bytes += natom * idxu_half_max * nelements * sizeof(double) * 2; // ylist
bytes += natom * nmax * idxu_max * 3 * sizeof(double) * 2; // dulist
bytes += natom * nmax * idxu_cache_max * 3 * sizeof(double) * 2; // dulist
#ifdef LMP_KOKKOS_GPU
}
#endif
@ -2063,6 +2195,8 @@ double SNAKokkos<DeviceType>::memory_usage()
bytes += jdim * jdim * jdim * sizeof(int); // idxcg_block
bytes += jdim * sizeof(int); // idxu_block
bytes += jdim * sizeof(int); // idxu_half_block
bytes += jdim * sizeof(int); // idxu_cache_block
bytes += jdim * jdim * jdim * sizeof(int); // idxz_block
bytes += jdim * jdim * jdim * sizeof(int); // idxb_block

4
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

@ -689,7 +689,7 @@ void PairLJCharmmCoulLong::init_style()
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
@ -718,7 +718,7 @@ void PairLJCharmmCoulLong::init_style()
// set & error check interior rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
cut_in_off = cut_respa[0];

9
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -545,10 +545,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
}
} else forcecoul = 0.0;
if (rsq <= cut_in_off_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else if (rsq <= cut_in_on_sq) {
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
}
@ -659,7 +656,7 @@ void PairLJCutCoulLong::init_style()
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
@ -676,7 +673,7 @@ void PairLJCutCoulLong::init_style()
// set rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = NULL;

4
src/MISC/fix_gld.cpp
View File

@ -547,6 +547,7 @@ int FixGLD::unpack_exchange(int nlocal, double *buf)
int FixGLD::pack_restart(int i, double *buf)
{
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 3*prony_terms + 1;
for (int k = 0; k < 3*prony_terms; k=k+3)
{
@ -566,7 +567,8 @@ void FixGLD::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to the nth set of extended variables
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i< nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;

4
src/MISC/fix_ttm.cpp
View File

@ -656,6 +656,7 @@ void FixTTM::restart(char *buf)
int FixTTM::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = flangevin[i][0];
buf[2] = flangevin[i][1];
@ -672,7 +673,8 @@ void FixTTM::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;

6
src/MOLECULE/fix_cmap.cpp
View File

@ -1041,7 +1041,7 @@ void FixCMAP::read_data_header(char *line)
sscanf(line,BIGINT_FORMAT,&ncmap);
} else error->all(FLERR,"Invalid read data header line for fix cmap");
// didn't set in constructor b/c this fix could be defined
// didn't set in constructor because this fix could be defined
// before newton command
newton_bond = force->newton_bond;
@ -1291,6 +1291,7 @@ int FixCMAP::pack_restart(int i, double *buf)
buf[n++] = ubuf(crossterm_atom4[i][m]).d;
buf[n++] = ubuf(crossterm_atom5[i][m]).d;
}
// pack buf[0] this way because other fixes unpack it
buf[0] = n;
return n;
@ -1305,7 +1306,8 @@ void FixCMAP::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int n = 0;
for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]);

2
src/PERI/fix_peri_neigh.cpp
View File

@ -579,6 +579,7 @@ void FixPeriNeigh::restart(char *buf)
int FixPeriNeigh::pack_restart(int i, double *buf)
{
int m = 0;
// pack buf[0] this way b/c other fixes unpack it
if (isVES) buf[m++] = 4*npartner[i] + 4;
else if (isEPS) buf[m++] = 3*npartner[i] + 5;
else buf[m++] = 2*npartner[i] + 4;
@ -608,6 +609,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way b/c other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

33
src/USER-CGSDK/pair_lj_sdk_coul_long.cpp
View File

@ -593,17 +593,32 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
if (rsq < cut_ljsq[itype][jtype]) {
const int ljt = lj_type[itype][jtype];
const double ljpow1 = lj_pow1[ljt];
const double ljpow2 = lj_pow2[ljt];
const double ljpref = lj_prefact[ljt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv
- lj2[itype][jtype]);
forcelj = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)
- ljpow2*pow(ratio,ljpow2))/rsq;
philj = factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))
- offset[itype][jtype]);
philj = r4inv*(lj3[itype][jtype]*r4inv*r4inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
forcelj = r6inv*(lj1[itype][jtype]*r3inv
- lj2[itype][jtype]);
philj = r6inv*(lj3[itype][jtype]*r3inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(lj1[itype][jtype]*r6inv
- lj2[itype][jtype]);
philj = r6inv*(lj3[itype][jtype]*r6inv
- lj4[itype][jtype]) - offset[itype][jtype];
}
forcelj *= factor_lj;
philj *= factor_lj;
}
fforce = (forcecoul + forcelj) * r2inv;

33
src/USER-CGSDK/pair_lj_sdk_coul_msm.cpp
View File

@ -257,17 +257,32 @@ double PairLJSDKCoulMSM::single(int i, int j, int itype, int jtype,
if (rsq < cut_ljsq[itype][jtype]) {
const int ljt = lj_type[itype][jtype];
const double ljpow1 = lj_pow1[ljt];
const double ljpow2 = lj_pow2[ljt];
const double ljpref = lj_prefact[ljt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv
- lj2[itype][jtype]);
forcelj = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)
- ljpow2*pow(ratio,ljpow2))/rsq;
philj = factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))
- offset[itype][jtype]);
philj = r4inv*(lj3[itype][jtype]*r4inv*r4inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
forcelj = r6inv*(lj1[itype][jtype]*r3inv
- lj2[itype][jtype]);
philj = r6inv*(lj3[itype][jtype]*r3inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(lj1[itype][jtype]*r6inv
- lj2[itype][jtype]);
philj = r6inv*(lj3[itype][jtype]*r6inv
- lj4[itype][jtype]) - offset[itype][jtype];
}
forcelj *= factor_lj;
philj *= factor_lj;
}
fforce = (forcecoul + forcelj) * r2inv;

42
src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp
View File

@ -691,38 +691,20 @@ void PairLJCharmmCoulLongSoft::init_style()
// request regular or rRESPA neighbor lists
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
if (respa == 0) irequest = neighbor->request(this,instance_me);
else if (respa == 1) {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
} else {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 2;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respamiddle = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
}
irequest = neighbor->request(this,instance_me);
} else irequest = neighbor->request(this,instance_me);
if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1;
neighbor->requests[irequest]->respainner = 1;
}
if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
// require cut_lj_inner < cut_lj
@ -739,7 +721,7 @@ void PairLJCharmmCoulLongSoft::init_style()
// set & error check interior rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
if (MIN(cut_lj,cut_coul) < cut_respa[3])

38
src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
View File

@ -634,44 +634,26 @@ void PairLJCutCoulLongSoft::init_style()
// request regular or rRESPA neighbor lists
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
int respa = 0;
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
if (respa == 0) irequest = neighbor->request(this,instance_me);
else if (respa == 1) {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
} else {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 2;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respamiddle = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
}
irequest = neighbor->request(this,instance_me);
} else irequest = neighbor->request(this,instance_me);
if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1;
neighbor->requests[irequest]->respainner = 1;
}
if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
cut_coulsq = cut_coul * cut_coul;
// set rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = NULL;

39
src/USER-FEP/pair_lj_cut_soft.cpp
View File

@ -512,46 +512,27 @@ void PairLJCutSoft::init_style()
// request regular or rRESPA neighbor lists
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
int respa = 0;
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
if (respa == 0) irequest = neighbor->request(this,instance_me);
else if (respa == 1) {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
} else {
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 2;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respamiddle = 1;
irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
}
irequest = neighbor->request(this,instance_me);
} else irequest = neighbor->request(this,instance_me);
if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1;
neighbor->requests[irequest]->respainner = 1;
}
if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
// set rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = NULL;
}
/* ----------------------------------------------------------------------

31
src/USER-FEP/pair_morse_soft.cpp
View File

@ -333,6 +333,37 @@ void PairMorseSoft::read_restart(FILE *fp)
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairMorseSoft::write_restart_settings(FILE *fp)
{
fwrite(&nlambda,sizeof(double),1,fp);
fwrite(&shift_range,sizeof(double),1,fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairMorseSoft::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
utils::sfread(FLERR,&nlambda,sizeof(double),1,fp,NULL,error);
utils::sfread(FLERR,&shift_range,sizeof(double),1,fp,NULL,error);
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,NULL,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,NULL,error);
}
MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
MPI_Bcast(&shift_range,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file

2
src/USER-FEP/pair_morse_soft.h
View File

@ -38,6 +38,8 @@ class PairMorseSoft : public PairMorse {
virtual double init_one(int, int);
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);

1
src/USER-MISC/README
View File

@ -26,6 +26,7 @@ angle_style dipole, Mario Orsi, orsimario at gmail.com, 10 Jan 12
angle_style quartic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
bond_style special, David Nicholson, davidanich at gmail.com, 31 Jan 2020
compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017

222
src/USER-MISC/bond_special.cpp Normal file
View File

@ -0,0 +1,222 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
Contributing Author: David Nicholson (MIT)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "bond_special.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondSpecial::BondSpecial(LAMMPS *lmp) : Bond(lmp) {}
/* ---------------------------------------------------------------------- */
BondSpecial::~BondSpecial()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(factor_lj);
memory->destroy(factor_coul);
}
}
/* ---------------------------------------------------------------------- */
void BondSpecial::init_style()
{
if (force->pair == NULL) error->all(FLERR,"No pair style defined");
else if ((force->pair->single_enable == 0) || force->pair->manybody_flag)
error->all(FLERR,"Pair style does not support bond style special");
if (force->special_lj[1] != 0.0 || force->special_coul[1] != 0)
error->all(FLERR,"Invalid 1-2 setting for bond style special.");
if (force->special_angle != 1 && (force->special_lj[2] != 1.0 ||
force->special_coul[2] != 1.0))
error->all(FLERR,"Invalid 1-3 setting for bond style special.");
if (force->special_dihedral != 1 && (force->special_lj[3] != 1.0 ||
force->special_coul[3] != 1.0))
error->all(FLERR,"Invalid 1-4 setting for bond style special.");
if (force->kspace != NULL)
error->all(FLERR,"Bond style special is not compatible with long range "
"Coulombic interactions");
}
/* ---------------------------------------------------------------------- */
void BondSpecial::compute(int eflag, int vflag)
{
int i1,i2,i1type,i2type,n,type;
double delx,dely,delz,ebond,fbond;
double rsq;
ev_init(eflag,vflag);
double **x = atom->x;
double **f = atom->f;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
for (n = 0; n < nbondlist; n++) {
i1 = bondlist[n][0];
i2 = bondlist[n][1];
type = bondlist[n][2];
i1type = atom->type[i1];
i2type = atom->type[i2];
delx = x[i1][0] - x[i2][0];
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
rsq = delx*delx + dely*dely + delz*delz;
ebond = force->pair->single(i1,i2,i1type,i2type,rsq,factor_lj[type],
factor_coul[type],fbond);
// apply force to each of 2 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fbond;
f[i1][1] += dely*fbond;
f[i1][2] += delz*fbond;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx*fbond;
f[i2][1] -= dely*fbond;
f[i2][2] -= delz*fbond;
}
if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
}
}
/* ---------------------------------------------------------------------- */
void BondSpecial::allocate()
{
allocated = 1;
int n = atom->nbondtypes;
memory->create(factor_lj,n+1,"bond:factor_lj");
memory->create(factor_coul,n+1,"bond:factor_coul");
memory->create(setflag,n+1,"bond:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void BondSpecial::coeff(int narg, char **arg)
{
if (narg != 3) error->all(FLERR,"Incorrect args for bond coefficients");
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi);
double factor_lj_one = force->numeric(FLERR,arg[1]);
double factor_coul_one = force->numeric(FLERR,arg[2]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
factor_lj[i] = factor_lj_one;
factor_coul[i] = factor_coul_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
}
/* ----------------------------------------------------------------------
return an equilbrium bond length
------------------------------------------------------------------------- */
double BondSpecial::equilibrium_distance(int i)
{
return 0.;
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void BondSpecial::write_restart(FILE *fp)
{
fwrite(&factor_lj[1],sizeof(double),atom->nbondtypes,fp);
fwrite(&factor_coul[1],sizeof(double),atom->nbondtypes,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void BondSpecial::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
fread(&factor_lj[1],sizeof(double),atom->nbondtypes,fp);
fread(&factor_coul[1],sizeof(double),atom->nbondtypes,fp);
}
MPI_Bcast(&factor_lj[1],atom->nbondtypes,MPI_DOUBLE,0,world);
MPI_Bcast(&factor_coul[1],atom->nbondtypes,MPI_DOUBLE,0,world);
for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void BondSpecial::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++)
fprintf(fp,"%d %g %g\n",i,factor_lj[i],factor_coul[i]);
}
/* ---------------------------------------------------------------------- */
double BondSpecial::single(int type, double rsq, int i, int j,
double &fforce)
{
int itype = atom->type[i];
int jtype = atom->type[j];
double ebond;
ebond = force->pair->single(i,j,itype,jtype,rsq,factor_lj[type],
factor_coul[type],fforce);
return ebond;
}

78
src/USER-MISC/bond_special.h Normal file
View File

@ -0,0 +1,78 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
Contributing Author: David Nicholson (MIT)
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
BondStyle(special,BondSpecial)
#else
#ifndef LMP_BOND_SPECIAL_H
#define LMP_BOND_SPECIAL_H
#include <cstdio>
#include "bond.h"
namespace LAMMPS_NS {
class BondSpecial: public Bond {
public:
BondSpecial(class LAMMPS *);
virtual ~BondSpecial();
void init_style();
void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_data(FILE *);
double single(int, double, int, int, double &);
protected:
double *factor_lj,*factor_coul;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
E: Invalid 1-2 setting for bond style special.
Bond style special must be used with zero factors for 1-2 special bonds.
E: Invalid 1-3 setting for bond style special.
Bond style special must be used with 1.0 factors for 1-3 special bonds or the
angle keyword set to yes.
E: Invalid 1-4 setting for bond style special.
Bond style special must be used with 1.0 factors for 1-4 special bonds or the
dihedral keyword set to yes.
E: Bond style special is not compatible with long range Coulombic interactions.
Self-explanatory.
*/

4
src/USER-MISC/fix_gle.cpp
View File

@ -822,6 +822,7 @@ int FixGLE::unpack_exchange(int nlocal, double *buf)
int FixGLE::pack_restart(int i, double *buf)
{
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 3*ns + 1;
for (int k = 0; k < 3*ns; k=k+3)
{
@ -841,7 +842,8 @@ void FixGLE::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to the nth set of extended variables
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i< nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;

2
src/USER-MISC/fix_pimd.cpp
View File

@ -794,6 +794,7 @@ int FixPIMD::pack_restart(int i, double *buf)
{
int offset=0;
int pos = i * 3;
// pack buf[0] this way because other fixes unpack it
buf[offset++] = size_peratom_cols+1;
memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1;
@ -811,6 +812,7 @@ void FixPIMD::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i=0; i<nth; i++) m += static_cast<int> (extra[nlocal][m]);

9
src/USER-MISC/fix_srp.cpp
View File

@ -114,7 +114,7 @@ void FixSRP::init()
error->all(FLERR,"Illegal bond particle type");
// this fix must come before any fix which migrates atoms in its pre_exchange()
// b/c this fix's pre_exchange() creates per-atom data structure
// because this fix's pre_exchange() creates per-atom data structure
// that data must be current for atom migration to carry it along
for (int i = 0; i < modify->nfix; i++) {
@ -555,6 +555,7 @@ void FixSRP::post_run()
int FixSRP::pack_restart(int i, double *buf)
{
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 3;
buf[m++] = array[i][0];
buf[m++] = array[i][1];
@ -569,10 +570,12 @@ void FixSRP::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++){
m += extra[nlocal][m];
m += static_cast<int> (extra[nlocal][m]);
}
m++;

2
src/USER-MISC/fix_ti_spring.cpp
View File

@ -315,6 +315,7 @@ int FixTISpring::unpack_exchange(int nlocal, double *buf)
int FixTISpring::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1];
@ -331,6 +332,7 @@ void FixTISpring::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

2
src/USER-MISC/fix_ttm_mod.cpp
View File

@ -980,6 +980,7 @@ void FixTTMMod::restart(char *buf)
int FixTTMMod::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = flangevin[i][0];
buf[2] = flangevin[i][1];
@ -996,6 +997,7 @@ void FixTTMMod::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

2
src/USER-MISC/pair_coul_diel.cpp
View File

@ -281,6 +281,8 @@ void PairCoulDiel::read_restart(FILE *fp)
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&rme[i][j],sizeof(double),1,fp,NULL,error);

51
src/USER-MISC/pair_coul_shield.cpp
View File

@ -115,8 +115,8 @@ void PairCoulShield::compute(int eflag, int vflag)
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
Tap = calc_Tap(r,cut[itype][jtype]);
dTap = calc_dTap(r,cut[itype][jtype]);
} else {Tap = 1.0; dTap = 0.0;}
forcecoul = qqrd2e*qtmp*q[j]*r*depsdr;
@ -297,6 +297,8 @@ void PairCoulShield::read_restart(FILE *fp)
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&sigmae[i][j],sizeof(double),1,fp,NULL,error);
@ -338,34 +340,41 @@ void PairCoulShield::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */
double PairCoulShield::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r, rarg,Vc,fvc,forcecoul,phishieldec;
double r3,th,epsr,depsdr,Tap,dTap;
double *q = atom->q;
double qqrd2e = force->qqrd2e;
r = sqrt(rsq);
r3 = rsq*r;
rarg = 1.0/sigmae[itype][jtype];
th = r3 + MathSpecial::cube(rarg);
epsr = 1.0/pow(th,0.333333333333333333);
depsdr = epsr*epsr;
depsdr *= depsdr;
Vc = qqrd2e*q[i]*q[j]*epsr;
// only computed between different layers as indicated by different molecule ids.
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
} else {Tap = 1.0; dTap = 0.0;}
if (atom->molecule[i] == atom->molecule[j]) {
fforce = 0.0;
return 0.0;
}
forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;
fvc = forcecoul*Tap - Vc*dTap/r;
fforce = factor_coul*fvc;
r = sqrt(rsq);
r3 = rsq*r;
rarg = 1.0/sigmae[itype][jtype];
th = r3 + MathSpecial::cube(rarg);
epsr = 1.0/pow(th,0.333333333333333333);
depsdr = epsr*epsr;
depsdr *= depsdr;
Vc = qqrd2e*q[i]*q[j]*epsr;
if (tap_flag) phishieldec = factor_coul*Vc*Tap;
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,cut[itype][jtype]);
dTap = calc_dTap(r,cut[itype][jtype]);
} else {Tap = 1.0; dTap = 0.0;}
forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;
fvc = forcecoul*Tap - Vc*dTap/r;
fforce = factor_coul*fvc;
if (tap_flag) phishieldec = Vc*Tap;
else phishieldec = Vc - offset[itype][jtype];
return factor_coul*phishieldec;
}

4
src/USER-MISC/pair_e3b.cpp
View File

@ -382,8 +382,6 @@ void PairE3B::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
typeO=force->inumeric(FLERR,arg[0]);
if (typeO<1 || typeO>atom->ntypes)
error->all(FLERR,"Invalid Otype: out of bounds");
}
/* ----------------------------------------------------------------------
@ -477,6 +475,8 @@ void PairE3B::init_style()
error->all(FLERR,"Pair style E3B requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style E3B requires newton pair on");
if (typeO<1 || typeO>atom->ntypes)
error->all(FLERR,"Invalid Otype: out of bounds");
// need a half neighbor list
neighbor->request(this,instance_me);

11
src/USER-MISC/pair_edip.cpp
View File

@ -51,7 +51,13 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairEDIP::PairEDIP(LAMMPS *lmp) : Pair(lmp)
PairEDIP::PairEDIP(LAMMPS *lmp) :
Pair(lmp), preInvR_ij(NULL), preExp3B_ij(NULL), preExp3BDerived_ij(NULL),
preExp2B_ij(NULL), preExp2BDerived_ij(NULL), prePow2B_ij(NULL),
preForceCoord(NULL), cutoffFunction(NULL), cutoffFunctionDerived(NULL),
pow2B(NULL), exp2B(NULL), exp3B(NULL), qFunctionGrid(NULL),
expMinusBetaZeta_iZeta_iGrid(NULL), tauFunctionGrid(NULL),
tauFunctionDerivedGrid(NULL)
{
single_enable = 0;
restartinfo = 0;
@ -503,6 +509,8 @@ void PairEDIP::allocateGrids(void)
double maxArgumentExpMinusBetaZeta_iZeta_i;
double const leftLimitToZero = -DBL_MIN * 1000.0;
deallocateGrids();
// tauFunctionGrid
maxArgumentTauFunctionGrid = leadDimInteractionList;
@ -561,6 +569,7 @@ void PairEDIP::allocatePreLoops(void)
{
int nthreads = comm->nthreads;
deallocatePreLoops();
memory->create(preInvR_ij,nthreads*leadDimInteractionList,"edip:preInvR_ij");
memory->create(preExp3B_ij,nthreads*leadDimInteractionList,"edip:preExp3B_ij");
memory->create(preExp3BDerived_ij,nthreads*leadDimInteractionList,

10
src/USER-MISC/pair_lennard_mdf.cpp
View File

@ -226,8 +226,8 @@ void PairLennardMDF::coeff(int narg, char **arg)
cut_inner_one = force->numeric(FLERR,arg[4]);
cut_one = force->numeric(FLERR,arg[5]);
}
if (cut_inner_global <= 0.0 || cut_inner_global > cut_global)
error->all(FLERR,"Illegal pair_style command");
if (cut_inner_one <= 0.0 || cut_inner_one > cut_one)
error->all(FLERR,"Illegal pair_coeff command");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
@ -328,6 +328,8 @@ void PairLennardMDF::read_restart(FILE *fp)
void PairLennardMDF::write_restart_settings(FILE *fp)
{
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&cut_inner_global,sizeof(double),1,fp);
}
/* ----------------------------------------------------------------------
@ -339,8 +341,12 @@ void PairLennardMDF::read_restart_settings(FILE *fp)
int me = comm->me;
if (me == 0) {
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,NULL,error);
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,NULL,error);
utils::sfread(FLERR,&cut_inner_global,sizeof(double),1,fp,NULL,error);
}
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
}
/* ----------------------------------------------------------------------

48
src/USER-MISC/pair_lj_expand_coul_long.cpp
View File

@ -53,6 +53,7 @@ PairLJExpandCoulLong::PairLJExpandCoulLong(LAMMPS *lmp) : Pair(lmp)
writedata = 1;
ftable = NULL;
qdist = 0.0;
cut_respa = NULL;
}
/* ---------------------------------------------------------------------- */
@ -167,7 +168,7 @@ void PairLJExpandCoulLong::compute(int eflag, int vflag)
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj = factor_lj*forcelj/rshift/r;
forcelj *= factor_lj/rshift/r;
} else forcelj = 0.0;
fpair = forcecoul*r2inv + forcelj;
@ -276,7 +277,7 @@ void PairLJExpandCoulLong::compute_inner()
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj = factor_lj*forcelj/rshift/r;
forcelj *= factor_lj/rshift/r;
} else forcelj = 0.0;
fpair = forcecoul*r2inv + forcelj;
@ -371,7 +372,7 @@ void PairLJExpandCoulLong::compute_middle()
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj = factor_lj*forcelj/rshift/r;
forcelj *= factor_lj/rshift/r;
} else forcelj = 0.0;
fpair = forcecoul*r2inv + forcelj;
@ -506,7 +507,7 @@ void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj = factor_lj*forcelj/rshift/r;
forcelj *= factor_lj/rshift/r;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
@ -524,6 +525,7 @@ void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
if (evflag) ecoul = evdwl = 0.0;
if (eflag) {
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
@ -538,14 +540,18 @@ void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
ecoul -= (1.0-factor_coul)*prefactor;
}
}
} else ecoul = 0.0;
}
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
rshift = r - shift[itype][jtype];
rshiftsq = rshift*rshift;
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
}
}
if (vflag) {
@ -564,15 +570,18 @@ void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
}
} else forcecoul = 0.0;
if (rsq <= cut_in_off_sq) {
r6inv = r2inv*r2inv*r2inv;
if (rsq < cut_ljsq[itype][jtype]) {
r = sqrt(rsq);
rshift = r - shift[itype][jtype];
rshiftsq = rshift*rshift;
rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else if (rsq <= cut_in_on_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
}
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}
forcelj *= factor_lj/rshift/r;
} else forcelj = 0.0;
fpair = forcecoul*r2inv + forcelj;
} else fpair = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);
@ -946,11 +955,16 @@ double PairLJExpandCoulLong::single(int i, int j, int itype, int jtype,
} else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
double rshift = r - shift[itype][jtype];
double rshiftsq = rshift*rshift;
double rshift2inv = 1.0/rshiftsq;
r6inv = rshift2inv*rshift2inv*rshift2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj *= factor_lj/rshift/r;
} else forcelj = 0.0;
fforce = (forcecoul + factor_lj*forcelj) * r2inv;
fforce = forcecoul*r2inv + forcelj;
double eng = 0.0;
if (rsq < cut_coulsq) {

40
src/USER-MISC/pair_meam_spline.cpp
View File

@ -57,6 +57,7 @@ PairMEAMSpline::PairMEAMSpline(LAMMPS *lmp) : Pair(lmp)
nelements = 0;
elements = NULL;
map = NULL;
Uprime_values = NULL;
nmax = 0;
@ -65,6 +66,14 @@ PairMEAMSpline::PairMEAMSpline(LAMMPS *lmp) : Pair(lmp)
comm_forward = 1;
comm_reverse = 0;
phis = NULL;
Us = NULL;
rhos = NULL;
fs = NULL;
gs = NULL;
zero_atom_energies = NULL;
}
/* ---------------------------------------------------------------------- */
@ -332,6 +341,8 @@ void PairMEAMSpline::allocate()
allocated = 1;
int n = nelements;
memory->destroy(setflag);
memory->destroy(cutsq);
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
@ -339,15 +350,23 @@ void PairMEAMSpline::allocate()
//Change the functional form
//f_ij->f_i
//g_i(cos\theta_ijk)->g_jk(cos\theta_ijk)
delete[] phis;
delete[] Us;
delete[] rhos;
delete[] fs;
delete[] gs;
phis = new SplineFunction[nmultichoose2];
Us = new SplineFunction[n];
rhos = new SplineFunction[n];
fs = new SplineFunction[n];
gs = new SplineFunction[nmultichoose2];
delete[] zero_atom_energies;
zero_atom_energies = new double[n];
delete[] map;
map = new int[n+1];
for (int i=0; i <= n; ++i) map[i] = -1;
}
/* ----------------------------------------------------------------------
@ -400,7 +419,7 @@ void PairMEAMSpline::coeff(int narg, char **arg)
if (strcmp(arg[i],elements[j]) == 0)
break;
if (j < nelements) map[i-2] = j;
else error->all(FLERR,"No matching element in EAM potential file");
else error->all(FLERR,"No matching element in meam/spline potential file");
}
}
// clear setflag since coeff() called once with I,J = * *
@ -419,8 +438,17 @@ void PairMEAMSpline::coeff(int narg, char **arg)
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
// check that each element is mapped to exactly one atom type
for (int i = 0; i < nelements; i++) {
count = 0;
for (int j = 1; j <= n; j++)
if (map[j] == i) count++;
if (count != 1)
error->all(FLERR,"Pair style meam/spline requires one atom type per element");
}
}
#define MAXLINE 1024
@ -460,6 +488,9 @@ void PairMEAMSpline::read_file(const char* filename)
if (nelements < 1)
error->one(FLERR, "Invalid number of atomic species on"
" meam/spline line in potential file");
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
elements = new char*[nelements];
for (int i=0; i<nelements; ++i) {
ptr = strtok(NULL," \t\n\r\f");
@ -480,7 +511,10 @@ void PairMEAMSpline::read_file(const char* filename)
}
nmultichoose2 = ((nelements+1)*nelements)/2;
// allocate!!
if (nelements != atom->ntypes)
error->all(FLERR,"Pair style meam/spline requires one atom type per element");
allocate();
// Parse spline functions.

11
src/USER-MISC/pair_morse_smooth_linear.cpp
View File

@ -201,14 +201,6 @@ void PairMorseSmoothLinear::coeff(int narg, char **arg)
alpha[i][j] = alpha_one;
r0[i][j] = r0_one;
cut[i][j] = cut_one;
morse1[i][j] = 2.0*d0[i][j]*alpha[i][j];
double alpha_dr = -alpha[i][j] * (cut[i][j] - r0[i][j]);
offset[i][j] = d0[i][j] * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr));
der_at_cutoff[i][j] = -2.0*alpha[i][j]*d0[i][j] * (exp(2.0*alpha_dr) - exp(alpha_dr));
setflag[i][j] = 1;
count++;
}
@ -227,6 +219,9 @@ double PairMorseSmoothLinear::init_one(int i, int j)
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
morse1[i][j] = 2.0*d0[i][j]*alpha[i][j];
double alpha_dr = -alpha[i][j] * (cut[i][j] - r0[i][j]);
offset[i][j] = d0[i][j] * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr));
der_at_cutoff[i][j] = -2.0*alpha[i][j]*d0[i][j] * (exp(2.0*alpha_dr) - exp(alpha_dr));
d0[j][i] = d0[i][j];
alpha[j][i] = alpha[i][j];

100
src/USER-OMP/ewald_omp.cpp
View File

@ -36,7 +36,7 @@ using namespace MathConst;
EwaldOMP::EwaldOMP(LAMMPS *lmp) : Ewald(lmp), ThrOMP(lmp, THR_KSPACE)
{
triclinic_support = 0;
triclinic_support = 1;
suffix_flag |= Suffix::OMP;
}
@ -79,7 +79,11 @@ void EwaldOMP::compute(int eflag, int vflag)
// partial structure factors on each processor
// total structure factor by summing over procs
eik_dot_r();
if (triclinic == 0)
eik_dot_r();
else
eik_dot_r_triclinic();
MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world);
@ -421,3 +425,95 @@ void EwaldOMP::eik_dot_r()
} // end of parallel region
}
/* ---------------------------------------------------------------------- */
void EwaldOMP::eik_dot_r_triclinic()
{
const double * const * const x = atom->x;
const double * const q = atom->q;
const int nlocal = atom->nlocal;
const int nthreads = comm->nthreads;
#if defined(_OPENMP)
#pragma omp parallel LMP_DEFAULT_NONE
#endif
{
int i,ifrom,ito,k,l,m,n,ic,tid;
double cstr1,sstr1;
double sqk,clpm,slpm;
double unitk_lamda[3];
loop_setup_thr(ifrom,ito,tid,nlocal,nthreads);
double max_kvecs[3];
max_kvecs[0] = kxmax;
max_kvecs[1] = kymax;
max_kvecs[2] = kzmax;
// (k,0,0), (0,l,0), (0,0,m)
for (ic = 0; ic < 3; ic++) {
unitk_lamda[0] = 0.0;
unitk_lamda[1] = 0.0;
unitk_lamda[2] = 0.0;
unitk_lamda[ic] = 2.0*MY_PI;
x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
sqk = unitk_lamda[ic]*unitk_lamda[ic];
if (sqk <= gsqmx) {
for (i = ifrom; i < ito; i++) {
cs[0][ic][i] = 1.0;
sn[0][ic][i] = 0.0;
cs[1][ic][i] = cos(unitk_lamda[0]*x[i][0] + unitk_lamda[1]*x[i][1] + unitk_lamda[2]*x[i][2]);
sn[1][ic][i] = sin(unitk_lamda[0]*x[i][0] + unitk_lamda[1]*x[i][1] + unitk_lamda[2]*x[i][2]);
cs[-1][ic][i] = cs[1][ic][i];
sn[-1][ic][i] = -sn[1][ic][i];
}
}
}
for (ic = 0; ic < 3; ic++) {
for (m = 2; m <= max_kvecs[ic]; m++) {
unitk_lamda[0] = 0.0;
unitk_lamda[1] = 0.0;
unitk_lamda[2] = 0.0;
unitk_lamda[ic] = 2.0*MY_PI*m;
x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
sqk = unitk_lamda[ic]*unitk_lamda[ic];
for (i = ifrom; i < ito; i++) {
cs[m][ic][i] = cs[m-1][ic][i]*cs[1][ic][i] -
sn[m-1][ic][i]*sn[1][ic][i];
sn[m][ic][i] = sn[m-1][ic][i]*cs[1][ic][i] +
cs[m-1][ic][i]*sn[1][ic][i];
cs[-m][ic][i] = cs[m][ic][i];
sn[-m][ic][i] = -sn[m][ic][i];
}
}
}
double * const sfacrl_thr = sfacrl + tid*kmax3d;
double * const sfacim_thr = sfacim + tid*kmax3d;
for (n = 0; n < kcount; n++) {
k = kxvecs[n];
l = kyvecs[n];
m = kzvecs[n];
cstr1 = 0.0;
sstr1 = 0.0;
for (i = ifrom; i < ito; i++) {
clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i];
slpm = sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i];
cstr1 += q[i]*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr1 += q[i]*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
}
sfacrl_thr[n] = cstr1;
sfacim_thr[n] = sstr1;
}
sync_threads();
data_reduce_thr(sfacrl,kmax3d,nthreads,1,tid);
data_reduce_thr(sfacim,kmax3d,nthreads,1,tid);
} // end of parallel region
}

1
src/USER-OMP/ewald_omp.h
View File

@ -34,6 +34,7 @@ namespace LAMMPS_NS {
protected:
virtual void eik_dot_r();
virtual void eik_dot_r_triclinic();
};
}

20
src/USER-OMP/pair_edip_omp.cpp
View File

@ -331,14 +331,14 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr)
f_ij[1] = forceMod2B * directorCos_ij_y;
f_ij[2] = forceMod2B * directorCos_ij_z;
f[j].x -= f_ij[0];
f[j].y -= f_ij[1];
f[j].z -= f_ij[2];
f[i].x += f_ij[0];
f[i].y += f_ij[1];
f[i].z += f_ij[2];
f[j].x -= f_ij[0];
f[j].y -= f_ij[1];
f[j].z -= f_ij[2];
// potential energy
evdwl = (exp2B_ij * potential2B_factor);
@ -460,13 +460,13 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr)
f_ij[1] = forceModCoord_ij * dr_ij[1];
f_ij[2] = forceModCoord_ij * dr_ij[2];
f[j].x -= f_ij[0];
f[j].y -= f_ij[1];
f[j].z -= f_ij[2];
f[i].x -= f_ij[0];
f[i].y -= f_ij[1];
f[i].z -= f_ij[2];
f[i].x += f_ij[0];
f[i].y += f_ij[1];
f[i].z += f_ij[2];
f[j].x += f_ij[0];
f[j].y += f_ij[1];
f[j].z += f_ij[2];
// potential energy

8
src/USER-OMP/pair_morse_smooth_linear_omp.cpp
View File

@ -35,8 +35,6 @@ PairMorseSmoothLinearOMP::PairMorseSmoothLinearOMP(LAMMPS *lmp) :
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairMorseSmoothLinearOMP::compute(int eflag, int vflag)
@ -76,8 +74,6 @@ void PairMorseSmoothLinearOMP::compute(int eflag, int vflag)
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairMorseSmoothLinearOMP::eval(int iifrom, int iito, ThrData * const thr)
{
@ -129,7 +125,7 @@ void PairMorseSmoothLinearOMP::eval(int iifrom, int iito, ThrData * const thr)
dexp = exp(-alpha[itype][jtype] * dr);
fpartial = morse1[itype][jtype] * (dexp*dexp - dexp) / r;
fpair = factor_lj * ( fpartial - der_at_cutoff[itype][jtype] / r);
fpair = factor_lj * ( fpartial + der_at_cutoff[itype][jtype] / r);
fxtmp += delx*fpair;
fytmp += dely*fpair;
@ -143,7 +139,7 @@ void PairMorseSmoothLinearOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EFLAG) {
evdwl = d0[itype][jtype] * (dexp*dexp - 2.0*dexp) -
offset[itype][jtype];
evdwl += ( r - cut[itype][jtype] ) * der_at_cutoff[itype][jtype];
evdwl -= ( r - cut[itype][jtype] ) * der_at_cutoff[itype][jtype];
evdwl *= factor_lj;
}

2
src/USER-SDPD/fix_meso_move.cpp
View File

@ -948,6 +948,7 @@ int FixMesoMove::unpack_exchange (int nlocal, double *buf) {
------------------------------------------------------------------------- */
int FixMesoMove::pack_restart (int i, double *buf) {
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1];
@ -963,6 +964,7 @@ void FixMesoMove::unpack_restart (int nlocal, int nth) {
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

2
src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp
View File

@ -450,6 +450,7 @@ int FixSMD_TLSPH_ReferenceConfiguration::unpack_exchange(int nlocal, double *buf
int FixSMD_TLSPH_ReferenceConfiguration::pack_restart(int i, double *buf) {
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 4 * npartner[i] + 2;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
@ -470,6 +471,7 @@ void FixSMD_TLSPH_ReferenceConfiguration::unpack_restart(int /*nlocal*/, int /*n
// ipage = NULL if being called from granular pair style init()
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
// double **extra = atom->extra;
//

14
src/create_atoms.cpp
View File

@ -71,6 +71,20 @@ void CreateAtoms::command(int narg, char **arg)
error->all(FLERR,"Cannot create_atoms after "
"reading restart file with per-atom info");
// check for compatible lattice
int latsty = domain->lattice->style;
if (domain->dimension == 2) {
if (latsty == Lattice::SC || latsty == Lattice::BCC
|| latsty == Lattice::FCC || latsty == Lattice::HCP
|| latsty == Lattice::DIAMOND)
error->all(FLERR,"Lattice style incompatible with simulation dimension");
} else {
if (latsty == Lattice::SQ ||latsty == Lattice::SQ2
|| latsty == Lattice::HEX)
error->all(FLERR,"Lattice style incompatible with simulation dimension");
}
// parse arguments
if (narg < 2) error->all(FLERR,"Illegal create_atoms command");

2
src/fix_move.cpp
View File

@ -1193,6 +1193,7 @@ int FixMove::pack_restart(int i, double *buf)
buf[n++] = qoriginal[i][2];
buf[n++] = qoriginal[i][3];
}
// pack buf[0] this way because other fixes unpack it
buf[0] = n;
return n;
}
@ -1206,6 +1207,7 @@ void FixMove::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

18
src/fix_neigh_history.cpp
View File

@ -149,7 +149,7 @@ void FixNeighHistory::init()
error->all(FLERR,"Neighbor history requires atoms have IDs");
// this fix must come before any fix which migrates atoms in its pre_exchange()
// b/c this fix's pre_exchange() creates per-atom data structure
// because this fix's pre_exchange() creates per-atom data structure
// that data must be current for atom migration to carry it along
for (int i = 0; i < modify->nfix; i++) {
@ -215,7 +215,7 @@ void FixNeighHistory::setup_post_neighbor()
called during run before atom exchanges, including for restart files
called at end of run via post_run()
do not call during setup of run (setup_pre_exchange)
b/c there is no guarantee of a current NDS (even on continued run)
because there is no guarantee of a current NDS (even on continued run)
if run command does a 2nd run with pre = no, then no neigh list
will be built, but old neigh list will still have the info
onesided and newton on and newton off versions
@ -242,7 +242,7 @@ void FixNeighHistory::pre_exchange_onesided()
double *allvalues,*onevalues;
// NOTE: all operations until very end are with nlocal_neigh <= current nlocal
// b/c previous neigh list was built with nlocal_neigh
// because previous neigh list was built with nlocal_neigh
// nlocal can be larger if other fixes added atoms at this pre_exchange()
// clear two paged data structures
@ -335,7 +335,7 @@ void FixNeighHistory::pre_exchange_newton()
// NOTE: all operations until very end are with
// nlocal_neigh <= current nlocal and nall_neigh
// b/c previous neigh list was built with nlocal_neigh & nghost_neigh
// because previous neigh list was built with nlocal_neigh & nghost_neigh
// nlocal can be larger if other fixes added atoms at this pre_exchange()
// clear two paged data structures
@ -430,7 +430,7 @@ void FixNeighHistory::pre_exchange_newton()
// perform reverse comm to augment
// owned atom partner/valuepartner with ghost info
// use variable variant b/c size of packed data can be arbitrarily large
// use variable variant because size of packed data can be arbitrarily large
// if many touching neighbors for large particle
commflag = PERPARTNER;
@ -463,7 +463,7 @@ void FixNeighHistory::pre_exchange_no_newton()
double *allvalues,*onevalues,*jvalues;
// NOTE: all operations until very end are with nlocal_neigh <= current nlocal
// b/c previous neigh list was built with nlocal_neigh
// because previous neigh list was built with nlocal_neigh
// nlocal can be larger if other fixes added atoms at this pre_exchange()
// clear two paged data structures
@ -708,10 +708,10 @@ void FixNeighHistory::grow_arrays(int nmax)
void FixNeighHistory::copy_arrays(int i, int j, int /*delflag*/)
{
// just copy pointers for partner and valuepartner
// b/c can't overwrite chunk allocation inside ipage_atom,dpage_atom
// because can't overwrite chunk allocation inside ipage_atom,dpage_atom
// incoming atoms in unpack_exchange just grab new chunks
// so are orphaning chunks for migrating atoms
// OK, b/c will reset ipage_atom,dpage_atom on next reneighboring
// OK, because will reset ipage_atom,dpage_atom on next reneighboring
npartner[j] = npartner[i];
partner[j] = partner[i];
@ -853,6 +853,7 @@ int FixNeighHistory::pack_restart(int i, double *buf)
memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes);
m += dnum;
}
// pack buf[0] this way because other fixes unpack it
buf[0] = m;
return m;
}
@ -868,6 +869,7 @@ void FixNeighHistory::unpack_restart(int nlocal, int nth)
if (ipage_atom == NULL) allocate_pages();
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
double **extra = atom->extra;

4
src/fix_property_atom.cpp
View File

@ -192,7 +192,7 @@ int FixPropertyAtom::setmask()
void FixPropertyAtom::init()
{
// error if atom style has changed since fix was defined
// don't allow this b/c user could change to style that defines molecule,q
// don't allow this because user could change to style that defines molecule,q
if (strcmp(astyle,atom->atom_style) != 0)
error->all(FLERR,"Atom style was redefined after using fix property/atom");
@ -579,6 +579,7 @@ int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
int FixPropertyAtom::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = nvalue+1;
int m = 1;
@ -602,6 +603,7 @@ void FixPropertyAtom::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

13
src/fix_setforce.cpp
View File

@ -297,8 +297,14 @@ void FixSetForce::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
// set force to desired value on requested level, 0.0 on other levels
if (ilevel == ilevel_respa) post_force(vflag);
else {
if (ilevel == 0) foriginal_saved[0] = foriginal_saved[1] = foriginal_saved[2] = 0.0;
if (ilevel == ilevel_respa) {
post_force(vflag);
foriginal[0] += foriginal_saved[0];
foriginal[1] += foriginal_saved[1];
foriginal[2] += foriginal_saved[2];
} else {
Region *region = NULL;
if (iregion >= 0) {
region = domain->regions[iregion];
@ -313,6 +319,9 @@ void FixSetForce::post_force_respa(int vflag, int ilevel, int /*iloop*/)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
foriginal_saved[0] += f[i][0];
foriginal_saved[1] += f[i][1];
foriginal_saved[2] += f[i][2];
if (xstyle) f[i][0] = 0.0;
if (ystyle) f[i][1] = 0.0;
if (zstyle) f[i][2] = 0.0;

2
src/fix_setforce.h
View File

@ -45,7 +45,7 @@ class FixSetForce : public Fix {
char *xstr,*ystr,*zstr;
char *idregion;
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[3],foriginal_all[3];
double foriginal[3],foriginal_all[3],foriginal_saved[3];
int force_flag;
int nlevels_respa,ilevel_respa;

2
src/fix_spring_self.cpp
View File

@ -264,6 +264,7 @@ int FixSpringSelf::unpack_exchange(int nlocal, double *buf)
int FixSpringSelf::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1];
@ -280,6 +281,7 @@ void FixSpringSelf::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

6
src/fix_store.cpp
View File

@ -280,7 +280,8 @@ int FixStore::pack_restart(int i, double *buf)
return 1;
}
buf[0] = ubuf(nvalues+1).d;
// pack buf[0] this way because other fixes unpack it
buf[0] = nvalues+1;
if (vecflag) buf[1] = vstore[i];
else
for (int m = 0; m < nvalues; m++)
@ -299,9 +300,10 @@ void FixStore::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += (int) ubuf(extra[nlocal][m]).i;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
if (vecflag) vstore[nlocal] = extra[nlocal][m];

2
src/fix_store_state.cpp
View File

@ -611,6 +611,7 @@ int FixStoreState::unpack_exchange(int nlocal, double *buf)
int FixStoreState::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = nvalues+1;
for (int m = 0; m < nvalues; m++) buf[m+1] = values[i][m];
return nvalues+1;
@ -625,6 +626,7 @@ void FixStoreState::unpack_restart(int nlocal, int nth)
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

29
src/fix_temp_berendsen.cpp
View File

@ -43,6 +43,7 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
// Berendsen thermostat should be applied every step
restart_global = 1;
dynamic_group_allow = 1;
nevery = 1;
scalar_flag = 1;
@ -241,6 +242,34 @@ double FixTempBerendsen::compute_scalar()
return energy;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTempBerendsen::write_restart(FILE *fp)
{
int n = 0;
double list[1];
list[n++] = energy;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTempBerendsen::restart(char *buf)
{
double *list = (double *) buf;
energy = list[0];
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */

2
src/fix_temp_berendsen.h
View File

@ -34,6 +34,8 @@ class FixTempBerendsen : public Fix {
int modify_param(int, char **);
void reset_target(double);
double compute_scalar();
void write_restart(FILE *);
void restart(char *buf);
virtual void *extract(const char *, int &);
private:

41
src/fix_temp_csld.cpp
View File

@ -49,6 +49,7 @@ FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) :
// CSLD thermostat should be applied every step
restart_global = 1;
nevery = 1;
scalar_flag = 1;
global_freq = nevery;
@ -297,6 +298,46 @@ double FixTempCSLD::compute_scalar()
return energy;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTempCSLD::write_restart(FILE *fp)
{
int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy
double *list;
if (comm->me == 0) list = new double[nsize];
list[0] = energy;
list[1] = comm->nprocs;
double state[103];
random->get_state(state);
MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world);
if (comm->me == 0) {
int size = nsize * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),nsize,fp);
}
if (comm->me == 0) delete[] list;
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTempCSLD::restart(char *buf)
{
double *list = (double *) buf;
energy = list[0];
int nprocs = (int) list[1];
if (nprocs != comm->nprocs) {
if (comm->me == 0)
error->warning(FLERR,"Different number of procs. Cannot restore RNG state.");
} else random->set_state(list+2+comm->me*103);
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */

2
src/fix_temp_csld.h
View File

@ -34,6 +34,8 @@ class FixTempCSLD : public Fix {
int modify_param(int, char **);
void reset_target(double);
virtual double compute_scalar();
void write_restart(FILE *);
void restart(char *buf);
virtual void *extract(const char *, int &);
private:

40
src/fix_temp_csvr.cpp
View File

@ -126,6 +126,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
// CSVR thermostat should be applied every step
restart_global = 1;
nevery = 1;
scalar_flag = 1;
global_freq = nevery;
@ -331,6 +332,45 @@ double FixTempCSVR::compute_scalar()
return energy;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTempCSVR::write_restart(FILE *fp)
{
int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy
double *list;
if (comm->me == 0) list = new double[nsize];
list[0] = energy;
list[1] = comm->nprocs;
double state[103];
random->get_state(state);
MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world);
if (comm->me == 0) {
int size = nsize * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),nsize,fp);
}
if (comm->me == 0) delete[] list;
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTempCSVR::restart(char *buf)
{
double *list = (double *) buf;
energy = list[0];
int nprocs = (int) list[1];
if (nprocs != comm->nprocs) {
if (comm->me == 0)
error->warning(FLERR,"Different number of procs. Cannot restore RNG state.");
} else random->set_state(list+2+comm->me*103);
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */

2
src/fix_temp_csvr.h
View File

@ -34,6 +34,8 @@ class FixTempCSVR : public Fix {
int modify_param(int, char **);
void reset_target(double);
virtual double compute_scalar();
void write_restart(FILE *);
void restart(char *buf);
virtual void *extract(const char *, int &);
private:

7
src/lammps.cpp
View File

@ -854,9 +854,6 @@ void LAMMPS::destroy()
delete neighbor;
neighbor = NULL;
delete comm;
comm = NULL;
delete force;
force = NULL;
@ -870,6 +867,10 @@ void LAMMPS::destroy()
// since they delete fixes
modify = NULL;
delete comm; // comm must come after modify
// since fix destructors may access comm
comm = NULL;
delete domain; // domain must come after modify
// since fix destructors access domain
domain = NULL;

2
src/lattice.cpp
View File

@ -28,8 +28,6 @@ using namespace LAMMPS_NS;
#define BIG 1.0e30
enum{NONE,SC,BCC,FCC,HCP,DIAMOND,SQ,SQ2,HEX,CUSTOM};
/* ---------------------------------------------------------------------- */
Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)

2
src/lattice.h
View File

@ -20,6 +20,8 @@ namespace LAMMPS_NS {
class Lattice : protected Pointers {
public:
enum{NONE,SC,BCC,FCC,HCP,DIAMOND,SQ,SQ2,HEX,CUSTOM};
int style; // NONE,SC,FCC,etc
double xlattice,ylattice,zlattice; // lattice scale factors in 3 dims
double a1[3],a2[3],a3[3]; // edge vectors of unit cell

3
src/pair.cpp
View File

@ -86,6 +86,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
manybody_flag = 0;
offset_flag = 0;
mix_flag = GEOMETRIC;
mixed_flag = 1;
tail_flag = 0;
etail = ptail = etail_ij = ptail_ij = 0.0;
ncoultablebits = 12;
@ -251,10 +252,12 @@ void Pair::init()
cutforce = 0.0;
etail = ptail = 0.0;
mixed_flag = 1;
double cut;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if ((i != j) && setflag[i][j]) mixed_flag = 0;
cut = init_one(i,j);
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutforce = MAX(cutforce,cut);

1
src/pair.h
View File

@ -105,6 +105,7 @@ class Pair : protected Pointers {
int allocated; // 0/1 = whether arrays are allocated
// public so external driver can check
int compute_flag; // 0 if skip compute()
int mixed_flag; // 1 if all itype != jtype coeffs are from mixing
enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options

2
src/pair_coul_streitz.cpp
View File

@ -249,7 +249,7 @@ void PairCoulStreitz::read_file(char *file)
FILE *fp;
if (comm->me == 0) {
fp = fopen(file,"r");
fp = force->open_potential(file);
if (fp == NULL)
error->one(FLERR,fmt::format("Cannot open coul/streitz potential "
"file {}",file));

4
src/pair_lj_cut.cpp
View File

@ -484,7 +484,7 @@ void PairLJCut::init_style()
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
@ -499,7 +499,7 @@ void PairLJCut::init_style()
// set rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
if (utils::strmatch(update->integrate_style,"^respa") &&
((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = NULL;

14
src/pair_table.cpp
View File

@ -458,25 +458,25 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
}
if (ferror)
error->warning(FLERR,fmt::format("{} of {} force values in table are "
error->warning(FLERR,fmt::format("{} of {} force values in table {} are "
"inconsistent with -dE/dr.\n Should "
"only be flagged at inflection points",
ferror,tb->ninput));
ferror,tb->ninput,keyword));
// warn if re-computed distance values differ from file values
if (rerror)
error->warning(FLERR,fmt::format("{} of {} distance values in table with "
"relative error\n over {} to "
error->warning(FLERR,fmt::format("{} of {} distance values in table {} "
"with relative error\n over {} to "
"re-computed values",
rerror,tb->ninput,EPSILONR));
rerror,tb->ninput,EPSILONR,keyword));
// warn if data was read incompletely, e.g. columns were missing
if (cerror)
error->warning(FLERR,fmt::format("{} of {} lines in table were "
error->warning(FLERR,fmt::format("{} of {} lines in table {} were "
"incomplete\n or could not be parsed "
"completely",cerror,tb->ninput));
"completely",cerror,tb->ninput,keyword));
}
/* ----------------------------------------------------------------------

36
src/random_mars.cpp
View File

@ -16,6 +16,7 @@
#include "random_mars.h"
#include <cmath>
#include <cstring>
#include "error.h"
#include "math_const.h"
@ -36,6 +37,7 @@ RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp),
save = 0;
u = new double[97+1];
memset(u,0,98*sizeof(double));
ij = (seed-1)/30082;
kl = (seed-1) - 30082*ij;
@ -181,7 +183,7 @@ double RanMars::besselexp(double theta, double alpha, double cp)
void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
{
int mode,index,oldindex,newvalue,nflip,which,niter;
int active[2],first[2],last[2];
int active[2],first[2];
int newactive[2],newfirst[2],newlast[2];
bigint nmark,nflipall;
bigint activeall[2],bsum[3],bsumall[3];
@ -191,8 +193,6 @@ void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
active[1] = 0;
first[0] = 0;
first[1] = -1;
last[0] = nmine-1;
last[1] = -1;
bigint bnmine = nmine;
bigint bnall;
@ -266,8 +266,6 @@ void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
active[1] = newactive[1];
first[0] = newfirst[0];
first[1] = newfirst[1];
last[0] = newlast[0];
last[1] = newlast[1];
}
// update nmark and activeall
@ -291,3 +289,31 @@ void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
// niter,nmark,nactiveall,thresh,nflipall);
}
}
/* ----------------------------------------------------------------------
store state in buffer
------------------------------------------------------------------------- */
void RanMars::get_state(double *state)
{
for (int i=0; i < 98; ++i) state[i] = u[i];
state[98] = i97;
state[99] = j97;
state[100]= c;
state[101]= cd;
state[102]= cm;
}
/* ----------------------------------------------------------------------
restore state from buffer
------------------------------------------------------------------------- */
void RanMars::set_state(double *state)
{
for (int i=0; i < 98; ++i) u[i] = state[i];
i97 = state[98];
j97 = state[99];
c = state[100];
cd = state[101];
cm = state[102];
}

2
src/random_mars.h
View File

@ -28,6 +28,8 @@ class RanMars : protected Pointers {
double rayleigh(double sigma);
double besselexp(double theta, double alpha, double cp);
void select_subset(bigint, int, int *, int *);
void get_state(double *);
void set_state(double *);
private:
int save;

6
src/velocity.cpp
View File

@ -67,7 +67,7 @@ void Velocity::command(int narg, char **arg)
// check if velocities of atoms in rigid bodies are updated
if (modify->check_rigid_group_overlap(groupbit))
if (modify->check_rigid_group_overlap(groupbit) && (comm->me == 0))
error->warning(FLERR,"Changing velocities of atoms in rigid bodies. "
"This has no effect unless rigid bodies are rebuild");
@ -200,7 +200,7 @@ void Velocity::create(double t_desired, int seed)
// initialize temperature computation(s)
// warn if groups don't match
if (igroup != temperature->igroup && comm->me == 0)
if ((igroup != temperature->igroup) && (comm->me == 0))
error->warning(FLERR,"Mismatch between velocity and compute groups");
temperature->init();
temperature->setup();
@ -596,7 +596,7 @@ void Velocity::scale(int /*narg*/, char **arg)
// initialize temperature computation
// warn if groups don't match
if (igroup != temperature->igroup && comm->me == 0)
if ((igroup != temperature->igroup) && (comm->me == 0))
error->warning(FLERR,"Mismatch between velocity and compute groups");
temperature->init();
temperature->setup();

3
src/write_data.cpp
View File

@ -303,6 +303,9 @@ void WriteData::force_fields()
{
if (force->pair && force->pair->writedata) {
if (pairflag == II) {
if ((comm->me == 0) && (force->pair->mixed_flag == 0))
error->warning(FLERR,"Not all mixed pair coeffs generated from mixing. "
"Use write_data with 'pair ij' option to store all pair coeffs.");
fmt::print(fp,"\nPair Coeffs # {}\n\n", force->pair_style);
force->pair->write_data(fp);
} else if (pairflag == IJ) {

2
unittest/CMakeLists.txt
View File

@ -1,9 +1,9 @@
include(GTest)
add_subdirectory(force-styles)
add_subdirectory(utils)
add_subdirectory(formats)
add_subdirectory(commands)
add_subdirectory(force-styles)
find_package(ClangFormat 8.0)

4
unittest/commands/CMakeLists.txt
View File

@ -3,6 +3,10 @@ add_executable(test_simple_commands test_simple_commands.cpp)
target_link_libraries(test_simple_commands PRIVATE lammps GTest::GMock GTest::GTest)
add_test(NAME SimpleCommands COMMAND test_simple_commands WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
add_executable(test_lattice_region test_lattice_region.cpp)
target_link_libraries(test_lattice_region PRIVATE lammps GTest::GMock GTest::GTest)
add_test(NAME LatticeRegion COMMAND test_lattice_region WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
add_executable(test_kim_commands test_kim_commands.cpp)
target_link_libraries(test_kim_commands PRIVATE lammps GTest::GMock GTest::GTest)
add_test(NAME KimCommands COMMAND test_kim_commands WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})

703
unittest/commands/test_lattice_region.cpp Normal file
View File

@ -0,0 +1,703 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, 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 "atom.h"
#include "domain.h"
#include "fmt/format.h"
#include "info.h"
#include "input.h"
#include "lammps.h"
#include "lattice.h"
#include "region.h"
#include "utils.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include <cstdio>
#include <cstring>
#include <fstream>
#include <iostream>
#include <mpi.h>
// whether to print verbose output (i.e. not capturing LAMMPS screen output).
bool verbose = false;
#if defined(OMPI_MAJOR_VERSION)
const bool have_openmpi = true;
#else
const bool have_openmpi = false;
#endif
using LAMMPS_NS::utils::split_words;
namespace LAMMPS_NS {
using ::testing::ExitedWithCode;
using ::testing::MatchesRegex;
using ::testing::StrEq;
#define TEST_FAILURE(errmsg, ...) \
if (Info::has_exceptions()) { \
::testing::internal::CaptureStdout(); \
ASSERT_ANY_THROW({__VA_ARGS__}); \
auto mesg = ::testing::internal::GetCapturedStdout(); \
ASSERT_THAT(mesg, MatchesRegex(errmsg)); \
} else { \
if (!have_openmpi) { \
::testing::internal::CaptureStdout(); \
ASSERT_DEATH({__VA_ARGS__}, ""); \
auto mesg = ::testing::internal::GetCapturedStdout(); \
ASSERT_THAT(mesg, MatchesRegex(errmsg)); \
} \
}
class LatticeRegionTest : public ::testing::Test {
protected:
LAMMPS *lmp;
void SetUp() override
{
const char *args[] = {"LatticeRegionTest", "-log", "none", "-echo", "screen", "-nocite"};
char **argv = (char **)args;
int argc = sizeof(args) / sizeof(char *);
if (!verbose) ::testing::internal::CaptureStdout();
lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
lmp->input->one("units metal");
if (!verbose) ::testing::internal::GetCapturedStdout();
}
void TearDown() override
{
if (!verbose) ::testing::internal::CaptureStdout();
delete lmp;
if (!verbose) ::testing::internal::GetCapturedStdout();
}
};
TEST_F(LatticeRegionTest, lattice_none)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice none 2.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::NONE);
ASSERT_EQ(lattice->xlattice, 2.0);
ASSERT_EQ(lattice->ylattice, 2.0);
ASSERT_EQ(lattice->zlattice, 2.0);
ASSERT_EQ(lattice->nbasis, 0);
ASSERT_EQ(lattice->basis, nullptr);
TEST_FAILURE(".*ERROR: Illegal lattice command.*", lmp->input->one("lattice"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*", lmp->input->one("lattice xxx"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*", lmp->input->one("lattice none 1.0 origin"););
TEST_FAILURE(".*ERROR: Expected floating point.*", lmp->input->one("lattice none xxx"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("units lj");
lmp->input->one("lattice none 1.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lattice->xlattice, 1.0);
ASSERT_EQ(lattice->ylattice, 1.0);
ASSERT_EQ(lattice->zlattice, 1.0);
}
TEST_F(LatticeRegionTest, lattice_sc)
{
::testing::internal::CaptureStdout();
lmp->input->one("lattice sc 1.0 spacing 1.5 2.0 3.0");
auto output = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << output;
ASSERT_THAT(output, MatchesRegex(".*Lattice spacing in x,y,z = 1.50* 2.0* 3.0*.*"));
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->xlattice, 1.5);
ASSERT_EQ(lattice->ylattice, 2.0);
ASSERT_EQ(lattice->zlattice, 3.0);
::testing::internal::CaptureStdout();
lmp->input->one("lattice sc 2.0");
output = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << output;
ASSERT_THAT(output, MatchesRegex(".*Lattice spacing in x,y,z = 2.0* 2.0* 2.0*.*"));
lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::SC);
ASSERT_EQ(lattice->xlattice, 2.0);
ASSERT_EQ(lattice->ylattice, 2.0);
ASSERT_EQ(lattice->zlattice, 2.0);
ASSERT_EQ(lattice->nbasis, 1);
ASSERT_NE(lattice->basis, nullptr);
ASSERT_EQ(lattice->a1[0], 1.0);
ASSERT_EQ(lattice->a1[1], 0.0);
ASSERT_EQ(lattice->a1[2], 0.0);
ASSERT_EQ(lattice->a2[0], 0.0);
ASSERT_EQ(lattice->a2[1], 1.0);
ASSERT_EQ(lattice->a2[2], 0.0);
ASSERT_EQ(lattice->a3[0], 0.0);
ASSERT_EQ(lattice->a3[1], 0.0);
ASSERT_EQ(lattice->a3[2], 1.0);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice sc 1.0 origin 1.0 1.0 1.0"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice sc 1.0 origin 1.0"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice sc 1.0 origin 0.0 0.0 0.0 xxx"););
TEST_FAILURE(".*ERROR: Expected floating point.*",
lmp->input->one("lattice sc 1.0 origin xxx 1.0 1.0"););
TEST_FAILURE(".*ERROR: Lattice orient vectors are not orthogonal.*",
lmp->input->one("lattice sc 1.0 orient x 2 2 0"););
TEST_FAILURE(".*ERROR: Lattice orient vectors are not right-handed.*",
lmp->input->one("lattice sc 1.0 orient y 0 -1 0"););
TEST_FAILURE(".*ERROR: Lattice spacings are invalid.*",
lmp->input->one("lattice sc 1.0 spacing 0.0 1.0 1.0"););
TEST_FAILURE(".*ERROR: Lattice spacings are invalid.*",
lmp->input->one("lattice sc 1.0 spacing 1.0 -0.1 1.0"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("units lj");
lmp->input->one("lattice sc 2.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_DOUBLE_EQ(lattice->xlattice, pow(0.5, 1.0 / 3.0));
ASSERT_DOUBLE_EQ(lattice->ylattice, pow(0.5, 1.0 / 3.0));
ASSERT_DOUBLE_EQ(lattice->zlattice, pow(0.5, 1.0 / 3.0));
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice sc 1.0"););
}
TEST_F(LatticeRegionTest, lattice_bcc)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice bcc 4.2 orient x 1 1 0 orient y -1 1 0");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::BCC);
ASSERT_DOUBLE_EQ(lattice->xlattice, sqrt(2.0) * 4.2);
ASSERT_DOUBLE_EQ(lattice->ylattice, sqrt(2.0) * 4.2);
ASSERT_DOUBLE_EQ(lattice->zlattice, 4.2);
ASSERT_EQ(lattice->nbasis, 2);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.5);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.5);
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice bcc 1.0"););
}
TEST_F(LatticeRegionTest, lattice_fcc)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 3.5 origin 0.5 0.5 0.5");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::FCC);
ASSERT_DOUBLE_EQ(lattice->xlattice, 3.5);
ASSERT_DOUBLE_EQ(lattice->ylattice, 3.5);
ASSERT_DOUBLE_EQ(lattice->zlattice, 3.5);
ASSERT_EQ(lattice->nbasis, 4);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.5);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.0);
ASSERT_EQ(lattice->basis[2][0], 0.5);
ASSERT_EQ(lattice->basis[2][1], 0.0);
ASSERT_EQ(lattice->basis[2][2], 0.5);
ASSERT_EQ(lattice->basis[3][0], 0.0);
ASSERT_EQ(lattice->basis[3][1], 0.5);
ASSERT_EQ(lattice->basis[3][2], 0.5);
TEST_FAILURE(".*ERROR: Invalid option in lattice command for non-custom style.*",
lmp->input->one("lattice fcc 1.0 basis 0.0 0.0 0.0"););
TEST_FAILURE(".*ERROR: Invalid option in lattice command for non-custom style.*",
lmp->input->one("lattice fcc 1.0 a1 0.0 1.0 0.0"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice fcc 1.0 orient w 1 0 0"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice fcc 1.0"););
}
TEST_F(LatticeRegionTest, lattice_hcp)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice hcp 3.0 orient z 0 0 1");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::HCP);
ASSERT_DOUBLE_EQ(lattice->xlattice, 3.0);
ASSERT_DOUBLE_EQ(lattice->ylattice, 3.0 * sqrt(3.0));
ASSERT_DOUBLE_EQ(lattice->zlattice, 2.0 * sqrt(6.0));
ASSERT_EQ(lattice->nbasis, 4);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.5);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.0);
ASSERT_EQ(lattice->basis[2][0], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[2][1], 5.0 / 6.0);
ASSERT_EQ(lattice->basis[2][2], 0.5);
ASSERT_EQ(lattice->basis[3][0], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[3][1], 1.0 / 3.0);
ASSERT_EQ(lattice->basis[3][2], 0.5);
ASSERT_EQ(lattice->a1[0], 1.0);
ASSERT_EQ(lattice->a1[1], 0.0);
ASSERT_EQ(lattice->a1[2], 0.0);
ASSERT_EQ(lattice->a2[0], 0.0);
ASSERT_DOUBLE_EQ(lattice->a2[1], sqrt(3.0));
ASSERT_EQ(lattice->a2[2], 0.0);
ASSERT_EQ(lattice->a3[0], 0.0);
ASSERT_EQ(lattice->a3[1], 0.0);
ASSERT_DOUBLE_EQ(lattice->a3[2], sqrt(8.0 / 3.0));
TEST_FAILURE(".*ERROR: Invalid option in lattice command for non-custom style.*",
lmp->input->one("lattice hcp 1.0 a2 0.0 1.0 0.0"););
TEST_FAILURE(".*ERROR: Invalid option in lattice command for non-custom style.*",
lmp->input->one("lattice hcp 1.0 a3 0.0 1.0 0.0"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice hcp 1.0"););
}
TEST_F(LatticeRegionTest, lattice_diamond)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice diamond 4.1 orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::DIAMOND);
ASSERT_DOUBLE_EQ(lattice->xlattice, 6.6952719636073539);
ASSERT_DOUBLE_EQ(lattice->ylattice, 5.7982756057296889);
ASSERT_DOUBLE_EQ(lattice->zlattice, 7.1014083110323973);
ASSERT_EQ(lattice->nbasis, 8);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.0);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.5);
ASSERT_EQ(lattice->basis[2][0], 0.5);
ASSERT_EQ(lattice->basis[2][1], 0.0);
ASSERT_EQ(lattice->basis[2][2], 0.5);
ASSERT_EQ(lattice->basis[3][0], 0.5);
ASSERT_EQ(lattice->basis[3][1], 0.5);
ASSERT_EQ(lattice->basis[3][2], 0.0);
ASSERT_EQ(lattice->basis[4][0], 0.25);
ASSERT_EQ(lattice->basis[4][1], 0.25);
ASSERT_EQ(lattice->basis[4][2], 0.25);
ASSERT_EQ(lattice->basis[5][0], 0.25);
ASSERT_EQ(lattice->basis[5][1], 0.75);
ASSERT_EQ(lattice->basis[5][2], 0.75);
ASSERT_EQ(lattice->basis[6][0], 0.75);
ASSERT_EQ(lattice->basis[6][1], 0.25);
ASSERT_EQ(lattice->basis[6][2], 0.75);
ASSERT_EQ(lattice->basis[7][0], 0.75);
ASSERT_EQ(lattice->basis[7][1], 0.75);
ASSERT_EQ(lattice->basis[7][2], 0.25);
ASSERT_EQ(lattice->a1[0], 1.0);
ASSERT_EQ(lattice->a1[1], 0.0);
ASSERT_EQ(lattice->a1[2], 0.0);
ASSERT_EQ(lattice->a2[0], 0.0);
ASSERT_EQ(lattice->a2[1], 1.0);
ASSERT_EQ(lattice->a2[2], 0.0);
ASSERT_EQ(lattice->a3[0], 0.0);
ASSERT_EQ(lattice->a3[1], 0.0);
ASSERT_EQ(lattice->a3[2], 1.0);
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice diamond 1.0"););
}
TEST_F(LatticeRegionTest, lattice_sq)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
lmp->input->one("lattice sq 3.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::SQ);
ASSERT_DOUBLE_EQ(lattice->xlattice, 3.0);
ASSERT_DOUBLE_EQ(lattice->ylattice, 3.0);
ASSERT_DOUBLE_EQ(lattice->zlattice, 3.0);
ASSERT_EQ(lattice->nbasis, 1);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
TEST_FAILURE(
".*ERROR: Lattice settings are not compatible with 2d simulation.*",
lmp->input->one("lattice sq 1.0 orient x 1 1 2 orient y -1 1 0 orient z -1 -1 1"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 3");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice sq 1.0"););
}
TEST_F(LatticeRegionTest, lattice_sq2)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
lmp->input->one("lattice sq2 2.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::SQ2);
ASSERT_DOUBLE_EQ(lattice->xlattice, 2.0);
ASSERT_DOUBLE_EQ(lattice->ylattice, 2.0);
ASSERT_DOUBLE_EQ(lattice->zlattice, 2.0);
ASSERT_EQ(lattice->nbasis, 2);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.5);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.0);
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 3");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice sq2 1.0"););
}
TEST_F(LatticeRegionTest, lattice_hex)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
lmp->input->one("lattice hex 2.0");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::HEX);
ASSERT_DOUBLE_EQ(lattice->xlattice, 2.0);
ASSERT_DOUBLE_EQ(lattice->ylattice, 3.4641016151377544);
ASSERT_DOUBLE_EQ(lattice->zlattice, 2.0);
ASSERT_EQ(lattice->nbasis, 2);
ASSERT_EQ(lattice->basis[0][0], 0.0);
ASSERT_EQ(lattice->basis[0][1], 0.0);
ASSERT_EQ(lattice->basis[0][2], 0.0);
ASSERT_EQ(lattice->basis[1][0], 0.5);
ASSERT_EQ(lattice->basis[1][1], 0.5);
ASSERT_EQ(lattice->basis[1][2], 0.0);
ASSERT_EQ(lattice->a1[0], 1.0);
ASSERT_EQ(lattice->a1[1], 0.0);
ASSERT_EQ(lattice->a1[2], 0.0);
ASSERT_EQ(lattice->a2[0], 0.0);
ASSERT_DOUBLE_EQ(lattice->a2[1], sqrt(3.0));
ASSERT_EQ(lattice->a2[2], 0.0);
ASSERT_EQ(lattice->a3[0], 0.0);
ASSERT_EQ(lattice->a3[1], 0.0);
ASSERT_EQ(lattice->a3[2], 1.0);
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 3");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Lattice style incompatible with simulation dimension.*",
lmp->input->one("lattice hex 1.0"););
}
TEST_F(LatticeRegionTest, lattice_custom)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("variable a equal 4.34");
lmp->input->one("variable b equal $a*sqrt(3.0)");
lmp->input->one("variable c equal $a*sqrt(8.0/3.0)");
lmp->input->one("variable t equal 1.0/3.0");
lmp->input->one("variable f equal 5.0/6.0");
lmp->input->one("lattice custom 1.0 "
"a1 $a 0.0 0.0 "
"a2 0.0 $b 0.0 "
"a3 0.0 0.0 $c "
"basis 0.0 0.0 0.0 "
"basis 0.5 0.5 0.0 "
"basis $t 0.0 0.5 "
"basis $f 0.5 0.5 "
"basis 0.0 0.0 0.625 "
"basis 0.5 0.5 0.625 "
"basis $t 0.0 0.125 "
"basis $f 0.5 0.125 ");
if (!verbose) ::testing::internal::GetCapturedStdout();
auto lattice = lmp->domain->lattice;
ASSERT_EQ(lattice->style, Lattice::CUSTOM);
ASSERT_DOUBLE_EQ(lattice->xlattice, 4.34);
ASSERT_DOUBLE_EQ(lattice->ylattice, 4.34 * sqrt(3.0));
ASSERT_DOUBLE_EQ(lattice->zlattice, 4.34 * sqrt(8.0 / 3.0));
ASSERT_EQ(lattice->nbasis, 8);
ASSERT_DOUBLE_EQ(lattice->basis[0][0], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[0][1], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[0][2], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[1][0], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[1][1], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[1][2], 0.0);
ASSERT_NEAR(lattice->basis[2][0], 1.0 / 3.0, 1.0e-14);
ASSERT_DOUBLE_EQ(lattice->basis[2][1], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[2][2], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[3][0], 5.0 / 6.0);
ASSERT_DOUBLE_EQ(lattice->basis[3][1], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[3][2], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[4][0], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[4][1], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[4][2], 0.625);
ASSERT_DOUBLE_EQ(lattice->basis[5][0], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[5][1], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[5][2], 0.625);
ASSERT_NEAR(lattice->basis[6][0], 1.0 / 3.0, 1.0e-14);
ASSERT_DOUBLE_EQ(lattice->basis[6][1], 0.0);
ASSERT_DOUBLE_EQ(lattice->basis[6][2], 0.125);
ASSERT_DOUBLE_EQ(lattice->basis[7][0], 5.0 / 6.0);
ASSERT_DOUBLE_EQ(lattice->basis[7][1], 0.5);
ASSERT_DOUBLE_EQ(lattice->basis[7][2], 0.125);
ASSERT_DOUBLE_EQ(lattice->a1[0], 4.34);
ASSERT_DOUBLE_EQ(lattice->a1[1], 0.0);
ASSERT_DOUBLE_EQ(lattice->a1[2], 0.0);
ASSERT_DOUBLE_EQ(lattice->a2[0], 0.0);
ASSERT_DOUBLE_EQ(lattice->a2[1], 4.34 * sqrt(3.0));
ASSERT_DOUBLE_EQ(lattice->a2[2], 0.0);
ASSERT_DOUBLE_EQ(lattice->a3[0], 0.0);
ASSERT_DOUBLE_EQ(lattice->a3[1], 0.0);
ASSERT_DOUBLE_EQ(lattice->a3[2], 4.34 * sqrt(8.0 / 3.0));
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice custom 1.0 basis -0.1 0 0"););
TEST_FAILURE(".*ERROR: Illegal lattice command.*",
lmp->input->one("lattice custom 1.0 basis 0.0 1.0 0"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("dimension 2");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: No basis atoms in lattice.*", lmp->input->one("lattice custom 1.0"););
TEST_FAILURE(".*ERROR: Lattice settings are not compatible with 2d simulation.*",
lmp->input->one("lattice custom 1.0 origin 0.5 0.5 0.5 basis 0.0 0.0 0.0"););
TEST_FAILURE(".*ERROR: Lattice settings are not compatible with 2d simulation.*",
lmp->input->one("lattice custom 1.0 a1 1.0 1.0 1.0 basis 0.0 0.0 0.0"););
TEST_FAILURE(".*ERROR: Lattice settings are not compatible with 2d simulation.*",
lmp->input->one("lattice custom 1.0 a2 1.0 1.0 1.0 basis 0.0 0.0 0.0"););
TEST_FAILURE(".*ERROR: Lattice settings are not compatible with 2d simulation.*",
lmp->input->one("lattice custom 1.0 a3 1.0 1.0 1.0 basis 0.0 0.0 0.0"););
}
TEST_F(LatticeRegionTest, region_fail)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice none 2.0");
lmp->input->one("region box block 0 1 0 1 0 1");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Create_atoms command before simulation box is defined.*",
lmp->input->one("create_atoms 1 box"););
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("create_box 1 box");
if (!verbose) ::testing::internal::GetCapturedStdout();
TEST_FAILURE(".*ERROR: Cannot create atoms with undefined lattice.*",
lmp->input->one("create_atoms 1 box"););
}
TEST_F(LatticeRegionTest, region_block_lattice)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice sc 1.5");
lmp->input->one("region box block 0 2 0 2 0 2 units lattice");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 8);
ASSERT_DOUBLE_EQ(x[0][0], 0.0);
ASSERT_DOUBLE_EQ(x[0][1], 0.0);
ASSERT_DOUBLE_EQ(x[0][2], 0.0);
ASSERT_DOUBLE_EQ(x[1][0], 1.5);
ASSERT_DOUBLE_EQ(x[1][1], 0.0);
ASSERT_DOUBLE_EQ(x[1][2], 0.0);
ASSERT_DOUBLE_EQ(x[2][0], 0.0);
ASSERT_DOUBLE_EQ(x[2][1], 1.5);
ASSERT_DOUBLE_EQ(x[2][2], 0.0);
ASSERT_DOUBLE_EQ(x[3][0], 1.5);
ASSERT_DOUBLE_EQ(x[3][1], 1.5);
ASSERT_DOUBLE_EQ(x[3][2], 0.0);
}
TEST_F(LatticeRegionTest, region_block_box)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice sc 1.5 origin 0.75 0.75 0.75");
lmp->input->one("region box block 0 2 0 2 0 2 units box");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 1);
ASSERT_DOUBLE_EQ(x[0][0], 1.125);
ASSERT_DOUBLE_EQ(x[0][1], 1.125);
ASSERT_DOUBLE_EQ(x[0][2], 1.125);
}
TEST_F(LatticeRegionTest, region_cone)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region box cone x 1.0 1.0 0.5 2.1 0.0 2.0");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region box");
lmp->input->one("write_dump all atom init.lammpstrj");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 42);
}
TEST_F(LatticeRegionTest, region_cylinder)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region box cylinder z 1.0 1.0 2.1 0.0 2.0 ");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 114);
}
TEST_F(LatticeRegionTest, region_prism)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice bcc 2.5 origin 0.75 0.75 0.75");
lmp->input->one("region box prism 0 2 0 2 0 2 0.5 0.0 0.0");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 1);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 16);
}
TEST_F(LatticeRegionTest, region_sphere)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region box sphere 1.0 1.0 1.0 1.1");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 14);
}
TEST_F(LatticeRegionTest, region_union)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region part1 sphere 2.0 1.0 1.0 1.1");
lmp->input->one("region part2 block 0.0 2.0 0.0 2.0 0.0 2.0");
lmp->input->one("region box union 2 part1 part2");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 67);
}
TEST_F(LatticeRegionTest, region_intersect)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region part1 sphere 2.0 1.0 1.0 1.8");
lmp->input->one("region part2 block 0.0 2.0 0.0 2.0 0.0 2.0");
lmp->input->one("region box intersect 2 part1 part2");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region box");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 21);
}
TEST_F(LatticeRegionTest, region_plane)
{
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("lattice fcc 2.5 origin 0.5 0.5 0.5");
lmp->input->one("region box block 0.0 2.0 0.0 2.0 0.0 2.0");
lmp->input->one("region part1 plane 0.5 1.0 0.0 0.75 0.0 0.0");
lmp->input->one("region part2 plane 1.5 1.0 0.0 0.75 0.0 0.0 side out");
lmp->input->one("region atoms intersect 2 part1 part2");
lmp->input->one("create_box 1 box");
lmp->input->one("create_atoms 1 region atoms");
lmp->input->one("write_dump all atom init.lammpstrj");
if (!verbose) ::testing::internal::GetCapturedStdout();
ASSERT_EQ(lmp->domain->triclinic, 0);
auto x = lmp->atom->x;
ASSERT_EQ(lmp->atom->natoms, 16);
}
} // namespace LAMMPS_NS
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
::testing::InitGoogleMock(&argc, argv);
if (have_openmpi && !LAMMPS_NS::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 = 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;
}

2
unittest/commands/test_reset_ids.cpp
View File

@ -660,6 +660,8 @@ TEST_F(ResetIDsTest, TopologyData)
TEST_F(ResetIDsTest, DeathTests)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*", lmp->input->one("reset_mol_ids"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*",
lmp->input->one("reset_mol_ids all offset 1 1"););

13
unittest/force-styles/CMakeLists.txt
View File

@ -114,3 +114,16 @@ foreach(TEST ${KSPACE_TESTS})
add_test(NAME ${TNAME} COMMAND test_pair_style ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
endforeach()
# tester for timestepping fixes
add_executable(test_fix_timestep test_fix_timestep.cpp)
target_link_libraries(test_fix_timestep PRIVATE lammps style_tests)
# tests for timestep related fixes (time integration, thermostat, force manipulation, constraints/restraints)
file(GLOB FIX_TIMESTEP_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/fix-timestep-*.yaml)
foreach(TEST ${FIX_TIMESTEP_TESTS})
string(REGEX REPLACE "^.*fix-timestep-(.*)\.yaml" "FixTimestep:\\1" TNAME ${TEST})
add_test(NAME ${TNAME} COMMAND test_fix_timestep ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
endforeach()

9
unittest/force-styles/check_tests.py
View File

@ -55,6 +55,7 @@ kspace = {}
pair = {}
style_pattern = re.compile("(.+)Style\((.+),(.+)\)")
upper = re.compile("[A-Z]+")
gpu = re.compile("(.+)/gpu$")
intel = re.compile("(.+)/intel$")
kokkos = re.compile("(.+)/kk$")
@ -85,7 +86,7 @@ for header in headers:
for m in matches:
# skip over internal styles w/o explicit documentation
style = m[1]
if style.isupper():
if upper.match(style):
continue
# detect, process, and flag suffix styles:
@ -131,6 +132,8 @@ for header in headers:
register_style(kspace,style,info)
elif m[0] == 'Pair':
register_style(pair,style,info)
elif m[0] == 'Fix':
register_style(fix,style,info)
@ -179,6 +182,8 @@ counter += check_tests('improper',improper,'improper-*.yaml',
'.*improper_style:\s*((\S+).*)?')
counter += check_tests('kspace',kspace,'kspace-*.yaml',
'.*kspace_style\s*((\S+).*)?')
counter += check_tests('fix',fix,'fix-*.yaml',
' fix\s+((\S+)\s*)?')
total = len(pair)+len(bond)+len(angle)+len(dihedral)+len(improper)+len(kspace)
total = len(pair)+len(bond)+len(angle)+len(dihedral)+len(improper)+len(kspace)+len(fix)
print(f"\nTotal tests missing: {counter} of {total}")

1
unittest/force-styles/error_stats.h
View File

@ -29,6 +29,7 @@ public:
double dev() const;
double max() const { return maxerr; }
double idx() const { return maxidx; }
bool has_data() const { return num > 0; }
private:
double sum, sumsq, maxerr;

13
unittest/force-styles/test_config.h
View File

@ -58,15 +58,21 @@ public:
double run_coul;
stress_t init_stress;
stress_t run_stress;
double global_scalar;
std::vector<double> global_vector;
std::vector<coord_t> init_forces;
std::vector<coord_t> run_forces;
std::vector<coord_t> run_pos;
std::vector<coord_t> restart_pos;
std::vector<coord_t> run_vel;
std::vector<coord_t> restart_vel;
TestConfig() :
lammps_version(""), date_generated(""), basename(""), epsilon(1.0e-14), input_file(""),
pair_style("zero"), bond_style("zero"), angle_style("zero"), dihedral_style("zero"),
improper_style("zero"), kspace_style("none"), natoms(0), init_energy(0), run_energy(0),
init_vdwl(0), run_vdwl(0), init_coul(0), run_coul(0), init_stress({0, 0, 0, 0, 0, 0}),
run_stress({0, 0, 0, 0, 0, 0})
run_stress({0, 0, 0, 0, 0, 0}), global_scalar(0)
{
prerequisites.clear();
pre_commands.clear();
@ -79,6 +85,11 @@ public:
extract.clear();
init_forces.clear();
run_forces.clear();
run_pos.clear();
restart_pos.clear();
run_vel.clear();
restart_vel.clear();
global_vector.clear();
}
virtual ~TestConfig(){};

53
unittest/force-styles/test_config_reader.cpp
View File

@ -40,6 +40,8 @@ TestConfigReader::TestConfigReader(TestConfig &config) : YamlReader(), config(co
consumers["run_stress"] = &TestConfigReader::run_stress;
consumers["init_forces"] = &TestConfigReader::init_forces;
consumers["run_forces"] = &TestConfigReader::run_forces;
consumers["run_pos"] = &TestConfigReader::run_pos;
consumers["run_vel"] = &TestConfigReader::run_vel;
consumers["pair_style"] = &TestConfigReader::pair_style;
consumers["pair_coeff"] = &TestConfigReader::pair_coeff;
@ -48,6 +50,9 @@ TestConfigReader::TestConfigReader(TestConfig &config) : YamlReader(), config(co
consumers["run_vdwl"] = &TestConfigReader::run_vdwl;
consumers["run_coul"] = &TestConfigReader::run_coul;
consumers["global_scalar"] = &TestConfigReader::global_scalar;
consumers["global_vector"] = &TestConfigReader::global_vector;
consumers["bond_style"] = &TestConfigReader::bond_style;
consumers["bond_coeff"] = &TestConfigReader::bond_coeff;
consumers["angle_style"] = &TestConfigReader::angle_style;
@ -176,6 +181,36 @@ void TestConfigReader::run_forces(const yaml_event_t &event)
}
}
void TestConfigReader::run_pos(const yaml_event_t &event)
{
config.run_pos.clear();
config.run_pos.resize(config.natoms + 1);
std::stringstream data((char *)event.data.scalar.value);
std::string line;
while (std::getline(data, line, '\n')) {
int tag;
coord_t xyz;
sscanf(line.c_str(), "%d %lg %lg %lg", &tag, &xyz.x, &xyz.y, &xyz.z);
config.run_pos[tag] = xyz;
}
}
void TestConfigReader::run_vel(const yaml_event_t &event)
{
config.run_vel.clear();
config.run_vel.resize(config.natoms + 1);
std::stringstream data((char *)event.data.scalar.value);
std::string line;
while (std::getline(data, line, '\n')) {
int tag;
coord_t xyz;
sscanf(line.c_str(), "%d %lg %lg %lg", &tag, &xyz.x, &xyz.y, &xyz.z);
config.run_vel[tag] = xyz;
}
}
void TestConfigReader::pair_style(const yaml_event_t &event)
{
config.pair_style = (char *)event.data.scalar.value;
@ -267,3 +302,21 @@ void TestConfigReader::run_energy(const yaml_event_t &event)
{
config.run_energy = atof((char *)event.data.scalar.value);
}
void TestConfigReader::global_scalar(const yaml_event_t &event)
{
config.global_scalar = atof((char *)event.data.scalar.value);
}
void TestConfigReader::global_vector(const yaml_event_t &event)
{
std::stringstream data((char *)event.data.scalar.value);
config.global_vector.clear();
double value;
std::size_t num;
data >> num;
for (std::size_t i = 0; i < num; ++i) {
data >> value;
config.global_vector.push_back(value);
}
}

4
unittest/force-styles/test_config_reader.h
View File

@ -36,6 +36,8 @@ public:
void run_stress(const yaml_event_t &event);
void init_forces(const yaml_event_t &event);
void run_forces(const yaml_event_t &event);
void run_pos(const yaml_event_t &event);
void run_vel(const yaml_event_t &event);
void pair_style(const yaml_event_t &event);
void pair_coeff(const yaml_event_t &event);
void bond_style(const yaml_event_t &event);
@ -49,6 +51,8 @@ public:
void run_coul(const yaml_event_t &event);
void init_energy(const yaml_event_t &event);
void run_energy(const yaml_event_t &event);
void global_scalar(const yaml_event_t &event);
void global_vector(const yaml_event_t &event);
};
#endif

2
unittest/force-styles/test_error_stats.cpp
View File

@ -29,12 +29,14 @@ TEST(ErrorStats, test)
ASSERT_EQ(out.str(), "Average: 5.800e-01 StdDev: 7.305e-01 MaxErr: 2.000e+00 @ item: 3.0");
stats.reset();
ASSERT_EQ(stats.has_data(), false);
ASSERT_DOUBLE_EQ(stats.avg(), 0.0);
ASSERT_DOUBLE_EQ(stats.dev(), 0.0);
ASSERT_DOUBLE_EQ(stats.max(), 0.0);
ASSERT_EQ(stats.idx(), -1);
stats.add(1.0);
ASSERT_EQ(stats.has_data(), true);
ASSERT_DOUBLE_EQ(stats.avg(), 1.0);
ASSERT_DOUBLE_EQ(stats.dev(), 0.0);
ASSERT_DOUBLE_EQ(stats.max(), 1.0);

1166
unittest/force-styles/test_fix_timestep.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

5
unittest/force-styles/test_main.cpp
View File

@ -102,13 +102,16 @@ int main(int argc, char **argv)
if (strcmp(argv[iarg], "-g") == 0) {
if (iarg + 1 < argc) {
generate_yaml_file(argv[iarg + 1], test_config);
MPI_Finalize();
return 0;
} else {
usage(std::cerr, argv[0]);
MPI_Finalize();
return 1;
}
} else if (strcmp(argv[iarg], "-u") == 0) {
generate_yaml_file(argv[1], test_config);
MPI_Finalize();
return 0;
} else if (strcmp(argv[iarg], "-d") == 0) {
if (iarg + 1 < argc) {
@ -116,6 +119,7 @@ int main(int argc, char **argv)
iarg += 2;
} else {
usage(std::cerr, argv[0]);
MPI_Finalize();
return 1;
}
} else if (strcmp(argv[iarg], "-s") == 0) {
@ -127,6 +131,7 @@ int main(int argc, char **argv)
} else {
std::cerr << "unknown option: " << argv[iarg] << "\n\n";
usage(std::cerr, argv[0]);
MPI_Finalize();
return 1;
}
}

12
unittest/force-styles/test_pair_style.cpp
View File

@ -130,8 +130,9 @@ LAMMPS *init_lammps(int argc, char **argv, const TestConfig &cfg, const bool new
}
command("run 0 post no");
command("variable write_data_pair index ii");
command("write_restart " + cfg.basename + ".restart");
command("write_data " + cfg.basename + ".data");
command("write_data " + cfg.basename + ".data pair ${write_data_pair}");
command("write_coeff " + cfg.basename + "-coeffs.in");
return lmp;
@ -1101,12 +1102,12 @@ TEST(PairStyle, single)
// The single function in EAM is different from what we assume
// here, therefore we have to skip testing those pair styles.
// Pair style colloid is also not compatible with this single tester
// Pair styles colloid and yukawa/colloid are also not compatible with this single tester
if ((test_config.pair_style.substr(0, 7) == "colloid") ||
(test_config.pair_style.substr(0, 14) == "yukawa/colloid") ||
(test_config.pair_style.substr(0, 3) == "eam") ||
((test_config.pair_style.substr(0, 6) == "hybrid") &&
(test_config.pair_style.find("eam") != std::string::npos))
) {
(test_config.pair_style.find("eam") != std::string::npos))) {
if (!verbose) ::testing::internal::CaptureStdout();
cleanup_lammps(lmp, test_config);
if (!verbose) ::testing::internal::GetCapturedStdout();
@ -1161,11 +1162,14 @@ TEST(PairStyle, single)
}
// create (only) two atoms
command("mass * 1.0");
command("create_atoms 1 single 0.0 -0.75 0.4 units box");
command("create_atoms 2 single 1.5 0.25 -0.1 units box");
command("set atom 1 charge -0.5");
command("set atom 2 charge 0.5");
command("set atom 1 mol 1");
command("set atom 2 mol 2");
command("special_bonds lj/coul 1.0 1.0 1.0");
if (molecular) {

121
unittest/force-styles/tests/atomic-pair-colloid.yaml
View File

@ -1,12 +1,15 @@
---
lammps_version: 21 Jul 2020
date_generated: Mon Aug 3 16:06:24 202
date_generated: Fri Aug 7 11:57:40 202
epsilon: 5e-14
prerequisites: ! |
pair colloid
pre_commands: ! |
variable units index lj
variable atom_style index atomic
post_commands: ! |
mass 1 9.0
mass 2 1.0
neighbor 1.0 bin
comm_style brick
comm_modify mode single
@ -18,68 +21,68 @@ pair_coeff: ! |
2 2 10.0 1.0 0.0 0.0 2.5
extract: ! ""
natoms: 27
init_vdwl: -0.250265293649414
init_vdwl: -0.250320882681727
init_coul: 0
init_stress: ! |2-
1.6693554149342897e+00 1.1752557173260264e+01 1.0084916384631697e+00 -5.5164068587184527e+00 3.4448938909043297e-01 2.8524984331506240e+00
1.6687682615057471e+00 1.1747891073356403e+01 1.0083450567161805e+00 -5.5144333681331359e+00 3.4548099128807436e-01 2.8508471154795298e+00
init_forces: ! |2
1 7.2698902368314028e-02 7.6228314794484830e-02 1.0726717051410481e-01
2 -4.1487528743785314e-02 -9.8021436494101673e-03 1.5709289685835706e-01
3 -6.4514052467860251e-03 -1.0256365723614392e-03 -6.5497018575823654e-04
4 -1.5879716299045554e-02 -1.3313626328419375e-02 -1.1960220488237115e-02
5 -3.6475105049685382e-02 -2.2992583864092243e-02 1.6166097759476125e-02
6 5.4769180740248721e-04 -2.4426732691233789e-02 -3.6444902428178980e-02
7 1.2006858444817004e-03 8.6089326153587114e-04 4.8033026817420977e-05
8 6.2576206692511894e-03 1.1942031910109577e-03 8.5598267089392085e-04
9 1.8304752931477772e+00 -4.3380838069003760e+00 -1.2763141904695217e+00
10 1.0456768773213965e-01 -2.8797785965171453e-02 -1.0768387599117879e-02
11 7.4457152913473112e-04 5.3043528596594623e-05 1.1848179655097060e-04
12 1.8275344545112628e-04 4.1948901318387846e-04 9.8529270880004928e-05
13 -1.1259391902077470e-01 1.6396787556459839e-02 -9.0455786475091124e-04
14 -2.4411189619188785e-04 -5.8284002050384879e-05 1.8046434825690733e-05
15 -2.9421019290851403e-03 -7.4831176382853144e-04 -1.2326238660052020e-03
16 2.2461127490925330e+00 -1.1352527056039741e+00 2.3679980564815390e+00
17 -1.8842133286566962e-02 1.7091738246445273e-02 -1.4247568288438442e-02
18 1.0304762560041125e-01 7.8062023784307467e-03 1.9257495383127633e-02
19 -4.3351162132476757e-02 -2.2995928402219718e-02 -1.1392806036292676e-02
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.2501672350519291e-02 6.8582908647593869e-03 -6.8677246013076093e-02
22 -2.2945818607502835e+00 1.1603382884260454e+00 -2.3594708301087310e+00
23 2.5606115583434265e-02 -1.0270997331622232e-01 -1.1451933998746673e-01
24 -8.0648597970214809e-04 1.1943689027481706e-04 -9.8441554909815322e-04
25 6.1663239781204479e-04 -2.9236784591333024e-04 3.5387752678593650e-03
26 -1.7197958284241512e+00 4.3707144027711022e+00 1.2977194913045051e+00
27 -1.4110864281012736e-01 4.2418795982941855e-02 -6.2606997884263993e-02
run_vdwl: 0.0367090517931739
1 4.2502214537152215e-02 6.8580119035896456e-03 -6.8677100496229701e-02
2 -2.4409971231138453e-04 -5.8282143696384295e-05 1.8044935805594862e-05
3 -1.5879676810899343e-02 -1.3312213962585376e-02 -1.1959790230130810e-02
4 5.4778343227770802e-04 -2.4424984321614867e-02 -3.6445490162239011e-02
5 7.4458982140296120e-04 5.3043116325023629e-05 1.1848540218007466e-04
6 2.2466945989092872e+00 -1.1355562227620619e+00 2.3685617789648648e+00
7 -2.9420817632825609e-03 -7.4830177547552694e-04 -1.2326317038179185e-03
8 1.0456973881919887e-01 -2.8794438516942737e-02 -1.0767233651978719e-02
9 -1.8842568296905022e-02 1.7092235687070117e-02 -1.4248069827099535e-02
10 1.8278642893333091e-04 4.1941592242561833e-04 9.8624961124789412e-05
11 2.5605560774882609e-02 -1.0270993584697993e-01 -1.1452001476172027e-01
12 -1.7190177814447427e+00 4.3689733469585290e+00 1.2971816478430653e+00
13 -1.4111675314136360e-01 4.2420896699357280e-02 -6.2611143561577903e-02
14 -4.1487669742220991e-02 -9.8022075695912841e-03 1.5709450788901494e-01
15 1.8296904307488300e+00 -4.3363458812930853e+00 -1.2757754576582720e+00
16 1.2006895887469950e-03 8.6089636114693700e-04 4.8031493704482794e-05
17 6.1663662317927603e-04 -2.9235023723413710e-04 3.5387343984160454e-03
18 -3.6474070486619430e-02 -2.2991550342271461e-02 1.6164948715920169e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.4513560799741325e-03 -1.0256280577842596e-03 -6.5496554886125568e-04
21 6.2575756583543186e-03 1.1941945764063822e-03 8.5597303165342044e-04
22 -8.0641944111498032e-04 1.1943078425573997e-04 -9.8433130461350886e-04
23 -2.2951590782728504e+00 1.1606385208734131e+00 -2.3600356755963983e+00
24 7.2701432066201874e-02 7.6224922982964613e-02 1.0727212071804378e-01
25 -1.1259499193288544e-01 1.6395143300644177e-02 -9.0439962448757879e-04
26 1.0304882476745540e-01 7.8050496446500431e-03 1.9255643852370385e-02
27 -4.3346315050731818e-02 -2.2993111981455712e-02 -1.1392238078737215e-02
run_vdwl: 0.0367182211916192
run_coul: 0
run_stress: ! |2-
3.2547164843322682e+00 1.2822272288053410e+01 2.9381904014809601e+00 -6.5755478139675869e+00 2.3052669923975402e+00 1.5140903126232870e+00
3.2545234036053254e+00 1.2817594721674045e+01 2.9383804763809485e+00 -6.5738604535059864e+00 2.3065705273606478e+00 1.5121755374667603e+00
run_forces: ! |2
1 5.8196197291626200e-02 7.7034242946349835e-02 1.0456829776646272e-01
2 -4.1552988555246327e-02 -1.1540194679885485e-02 1.5467509805429100e-01
3 -6.6558201800471870e-03 -9.8090417095619679e-04 -4.5870088099424396e-04
4 -1.6184658518487705e-02 -1.2048418808706472e-02 -1.2190156970179262e-02
5 -3.3316803862839403e-02 -2.1131495092834619e-02 1.4963198384123102e-02
6 4.7038346343352116e-04 -3.0307837184768524e-02 -4.3608938525025766e-02
7 1.0396652127649081e-03 7.5317021954246196e-04 4.1930480907776011e-05
8 6.4423457305213877e-03 1.1468835414812518e-03 6.8526191116136061e-04
9 1.7585178109228492e+00 -4.3823551563819683e+00 -1.2468505588145766e+00
10 1.3392700382201897e-01 -2.9471624976996433e-02 -9.0915578208389906e-03
11 7.9839096035526878e-04 6.5881060273477888e-05 1.3840230606223551e-04
12 1.8540711169000120e-04 3.9935505493718929e-04 1.2187452986589205e-04
13 -1.4075050026169886e-01 1.8204883476184704e-02 -2.0674384501958810e-03
14 -2.4452271394244688e-04 -5.7393635529214755e-05 1.6908863929994165e-05
15 -3.4577717467029730e-03 -9.3350658134032588e-04 -1.5696623820874889e-03
16 5.4762334392971725e+00 -3.3427628201841180e+00 5.8008787111000641e+00
17 -1.7123087875076635e-02 1.5861563712562771e-02 -1.3049347364948664e-02
18 1.1434938302418582e-01 8.6112942994291432e-03 2.0292892978926169e-02
19 -3.8092989065812294e-02 -1.9963855292591683e-02 -9.8568020007968223e-03
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.1446743356056415e-02 8.4986882092498734e-03 -7.1565549648589480e-02
22 -5.5269240273286284e+00 3.3665609275390551e+00 -5.7932229791606380e+00
23 1.5450935830291834e-02 -1.0557581147017484e-01 -1.1539704280712981e-01
24 -7.5222270706101302e-04 1.1318969997713211e-04 -9.2841932954143037e-04
25 1.1778375109233723e-03 -7.1254208855309762e-05 3.8035457799329405e-03
26 -1.6497221843622811e+00 4.4233264004281976e+00 1.2763820593322663e+00
27 -1.3345796635606474e-01 3.6623792481486360e-02 -5.6711027332452278e-02
1 4.1447236896055568e-02 8.4984633363665834e-03 -7.1565468520811443e-02
2 -2.4451051586476304e-04 -5.7391827303752204e-05 1.6907428582411947e-05
3 -1.6184661894786571e-02 -1.2047019655367870e-02 -1.2189740041384602e-02
4 4.7049123492354846e-04 -3.0305522972900800e-02 -4.3609470856651497e-02
5 7.9841068285770012e-04 6.5880880038333371e-05 1.3840659114887740e-04
6 5.4775302106929926e+00 -3.3435741912382899e+00 5.8021270645808745e+00
7 -3.4577504657755014e-03 -9.3349532033408413e-04 -1.5696706937008277e-03
8 1.3392885595460452e-01 -2.9467417351771795e-02 -9.0902656642780702e-03
9 -1.7123469245219562e-02 1.5862013571050645e-02 -1.3049795876411903e-02
10 1.8544095678649998e-04 3.9928059116526140e-04 1.2197284521083842e-04
11 1.5450042897697011e-02 -1.0557539868014379e-01 -1.1539698970598398e-01
12 -1.6489596487183065e+00 4.4215309402003911e+00 1.2758457405639687e+00
13 -1.3347630087741680e-01 3.6628535027741112e-02 -5.6719354360201757e-02
14 -4.1553429683017021e-02 -1.1540368471833904e-02 1.5467778628554116e-01
15 1.7577484997124488e+00 -4.3805631271698307e+00 -1.2463135748050815e+00
16 1.0396685127639720e-03 7.5317296594891241e-04 4.1929174719427522e-05
17 1.1778403499163440e-03 -7.1238534940179336e-05 3.8035065505820104e-03
18 -3.3315903593767987e-02 -2.1130575404712599e-02 1.4962162894607431e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.6557683033616174e-03 -9.8089615075005066e-04 -4.5869764220629459e-04
21 6.4422979624128504e-03 1.1468754348007283e-03 6.8525347184504311e-04
22 -7.5216145678175033e-04 1.1318396788931372e-04 -9.2834078115899713e-04
23 -5.5282158633416767e+00 3.3673691649465551e+00 -5.7944723588252627e+00
24 5.8210356160171874e-02 7.7028326238927733e-02 1.0457588922155152e-01
25 -1.4075156519665147e-01 1.8202248246942915e-02 -2.0674716237672891e-03
26 1.1435050190904306e-01 8.6100147519947955e-03 2.0290916519561587e-02
27 -3.8088820630046961e-02 -1.9961457381631506e-02 -9.8563367312923152e-03
...

121
unittest/force-styles/tests/atomic-pair-colloid_multi.yaml
View File

@ -1,12 +1,15 @@
---
lammps_version: 21 Jul 2020
date_generated: Mon Aug 3 16:06:24 202
date_generated: Fri Aug 7 11:57:40 202
epsilon: 5e-14
prerequisites: ! |
pair colloid
pre_commands: ! |
variable units index lj
variable atom_style index atomic
post_commands: ! |
mass 1 9.0
mass 2 1.0
neighbor 1.0 multi
comm_style brick
comm_modify mode multi
@ -18,68 +21,68 @@ pair_coeff: ! |
2 2 10.0 1.0 0.0 0.0 2.5
extract: ! ""
natoms: 27
init_vdwl: -0.250265293649414
init_vdwl: -0.250320882681727
init_coul: 0
init_stress: ! |2-
1.6693554149342913e+00 1.1752557173260270e+01 1.0084916384631617e+00 -5.5164068587184527e+00 3.4448938909043431e-01 2.8524984331506258e+00
1.6687682615057460e+00 1.1747891073356403e+01 1.0083450567161796e+00 -5.5144333681331261e+00 3.4548099128807747e-01 2.8508471154795298e+00
init_forces: ! |2
1 7.2698902368314056e-02 7.6228314794484775e-02 1.0726717051410484e-01
2 -4.1487528743785314e-02 -9.8021436494101673e-03 1.5709289685835706e-01
3 -6.4514052467860251e-03 -1.0256365723614392e-03 -6.5497018575823654e-04
4 -1.5879716299045554e-02 -1.3313626328419375e-02 -1.1960220488237115e-02
5 -3.6475105049685402e-02 -2.2992583864092257e-02 1.6166097759476136e-02
6 5.4769180740248721e-04 -2.4426732691233789e-02 -3.6444902428178980e-02
7 1.2006858444817004e-03 8.6089326153587114e-04 4.8033026817420977e-05
8 6.2576206692511894e-03 1.1942031910109577e-03 8.5598267089392085e-04
9 1.8304752931477772e+00 -4.3380838069003760e+00 -1.2763141904695217e+00
10 1.0456768773213965e-01 -2.8797785965171453e-02 -1.0768387599117879e-02
11 7.4457152913473112e-04 5.3043528596594623e-05 1.1848179655097060e-04
12 1.8275344545112628e-04 4.1948901318387846e-04 9.8529270880004928e-05
13 -1.1259391902077470e-01 1.6396787556459839e-02 -9.0455786475091124e-04
14 -2.4411189619188785e-04 -5.8284002050384879e-05 1.8046434825690733e-05
15 -2.9421019290851403e-03 -7.4831176382853144e-04 -1.2326238660052020e-03
16 2.2461127490925330e+00 -1.1352527056039741e+00 2.3679980564815386e+00
17 -1.8842133286566962e-02 1.7091738246445273e-02 -1.4247568288438442e-02
18 1.0304762560041125e-01 7.8062023784307467e-03 1.9257495383127633e-02
19 -4.3351162132476757e-02 -2.2995928402219718e-02 -1.1392806036292676e-02
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.2501672350519284e-02 6.8582908647593938e-03 -6.8677246013076107e-02
22 -2.2945818607502839e+00 1.1603382884260454e+00 -2.3594708301087310e+00
23 2.5606115583434265e-02 -1.0270997331622232e-01 -1.1451933998746673e-01
24 -8.0648597970214809e-04 1.1943689027481706e-04 -9.8441554909815322e-04
25 6.1663239781204458e-04 -2.9236784591333030e-04 3.5387752678593654e-03
26 -1.7197958284241512e+00 4.3707144027711022e+00 1.2977194913045054e+00
27 -1.4110864281012736e-01 4.2418795982941855e-02 -6.2606997884263993e-02
run_vdwl: 0.0367090517931737
1 4.2502214537152215e-02 6.8580119035896499e-03 -6.8677100496229715e-02
2 -2.4409971231138453e-04 -5.8282143696384295e-05 1.8044935805594862e-05
3 -1.5879676810899343e-02 -1.3312213962585376e-02 -1.1959790230130810e-02
4 5.4778343227770802e-04 -2.4424984321614867e-02 -3.6445490162239011e-02
5 7.4458982140296174e-04 5.3043116325023304e-05 1.1848540218007477e-04
6 2.2466945989092868e+00 -1.1355562227620619e+00 2.3685617789648643e+00
7 -2.9420817632825613e-03 -7.4830177547552715e-04 -1.2326317038179183e-03
8 1.0456973881919887e-01 -2.8794438516942737e-02 -1.0767233651978719e-02
9 -1.8842568296905022e-02 1.7092235687070117e-02 -1.4248069827099535e-02
10 1.8278642893333091e-04 4.1941592242561833e-04 9.8624961124789412e-05
11 2.5605560774882609e-02 -1.0270993584697993e-01 -1.1452001476172027e-01
12 -1.7190177814447429e+00 4.3689733469585290e+00 1.2971816478430653e+00
13 -1.4111675314136360e-01 4.2420896699357280e-02 -6.2611143561577903e-02
14 -4.1487669742220991e-02 -9.8022075695912841e-03 1.5709450788901494e-01
15 1.8296904307488300e+00 -4.3363458812930853e+00 -1.2757754576582720e+00
16 1.2006895887469950e-03 8.6089636114693700e-04 4.8031493704482794e-05
17 6.1663662317927744e-04 -2.9235023723413574e-04 3.5387343984160432e-03
18 -3.6474070486619430e-02 -2.2991550342271461e-02 1.6164948715920169e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.4513560799741325e-03 -1.0256280577842596e-03 -6.5496554886125568e-04
21 6.2575756583543186e-03 1.1941945764063822e-03 8.5597303165342044e-04
22 -8.0641944111498032e-04 1.1943078425573997e-04 -9.8433130461350886e-04
23 -2.2951590782728504e+00 1.1606385208734131e+00 -2.3600356755963983e+00
24 7.2701432066201888e-02 7.6224922982964599e-02 1.0727212071804376e-01
25 -1.1259499193288544e-01 1.6395143300644177e-02 -9.0439962448757879e-04
26 1.0304882476745540e-01 7.8050496446500431e-03 1.9255643852370385e-02
27 -4.3346315050731818e-02 -2.2993111981455712e-02 -1.1392238078737215e-02
run_vdwl: 0.0367182211916192
run_coul: 0
run_stress: ! |2-
3.2547164843322691e+00 1.2822272288053391e+01 2.9381904014809663e+00 -6.5755478139675931e+00 2.3052669923975411e+00 1.5140903126232879e+00
3.2545234036053201e+00 1.2817594721674030e+01 2.9383804763809369e+00 -6.5738604535059864e+00 2.3065705273606469e+00 1.5121755374667567e+00
run_forces: ! |2
1 5.8196197291625715e-02 7.7034242946349946e-02 1.0456829776646254e-01
2 -4.1552988555246327e-02 -1.1540194679885485e-02 1.5467509805429100e-01
3 -6.6558201800471870e-03 -9.8090417095619679e-04 -4.5870088099424396e-04
4 -1.6184658518487705e-02 -1.2048418808706472e-02 -1.2190156970179262e-02
5 -3.3316803862839438e-02 -2.1131495092834646e-02 1.4963198384123119e-02
6 4.7038346343352116e-04 -3.0307837184768524e-02 -4.3608938525025766e-02
7 1.0396652127649081e-03 7.5317021954246196e-04 4.1930480907776011e-05
8 6.4423457305213877e-03 1.1468835414812518e-03 6.8526191116136061e-04
9 1.7585178109228492e+00 -4.3823551563819700e+00 -1.2468505588145766e+00
10 1.3392700382201897e-01 -2.9471624976996433e-02 -9.0915578208389906e-03
11 7.9839096035526857e-04 6.5881060273478105e-05 1.3840230606223543e-04
12 1.8540711169000120e-04 3.9935505493718929e-04 1.2187452986589205e-04
13 -1.4075050026169886e-01 1.8204883476184704e-02 -2.0674384501958810e-03
14 -2.4452271394244688e-04 -5.7393635529214755e-05 1.6908863929994165e-05
15 -3.4577717467029726e-03 -9.3350658134032599e-04 -1.5696623820874889e-03
16 5.4762334392971725e+00 -3.3427628201841184e+00 5.8008787111000641e+00
17 -1.7123087875076635e-02 1.5861563712562771e-02 -1.3049347364948664e-02
18 1.1434938302418582e-01 8.6112942994291432e-03 2.0292892978926169e-02
19 -3.8092989065812294e-02 -1.9963855292591683e-02 -9.8568020007968223e-03
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.1446743356056415e-02 8.4986882092498717e-03 -7.1565549648589466e-02
22 -5.5269240273286284e+00 3.3665609275390551e+00 -5.7932229791606371e+00
23 1.5450935830291834e-02 -1.0557581147017484e-01 -1.1539704280712981e-01
24 -7.5222270706101302e-04 1.1318969997713211e-04 -9.2841932954143037e-04
25 1.1778375109233736e-03 -7.1254208855308434e-05 3.8035457799329383e-03
26 -1.6497221843622811e+00 4.4233264004281976e+00 1.2763820593322663e+00
27 -1.3345796635606427e-01 3.6623792481486235e-02 -5.6711027332452084e-02
1 4.1447236896055561e-02 8.4984633363665869e-03 -7.1565468520811457e-02
2 -2.4451051586476304e-04 -5.7391827303752204e-05 1.6907428582411947e-05
3 -1.6184661894786571e-02 -1.2047019655367870e-02 -1.2189740041384602e-02
4 4.7049123492354846e-04 -3.0305522972900800e-02 -4.3609470856651497e-02
5 7.9841068285770055e-04 6.5880880038332991e-05 1.3840659114887751e-04
6 5.4775302106929935e+00 -3.3435741912382895e+00 5.8021270645808736e+00
7 -3.4577504657755018e-03 -9.3349532033408424e-04 -1.5696706937008275e-03
8 1.3392885595460452e-01 -2.9467417351771795e-02 -9.0902656642780702e-03
9 -1.7123469245219562e-02 1.5862013571050645e-02 -1.3049795876411903e-02
10 1.8544095678649998e-04 3.9928059116526140e-04 1.2197284521083842e-04
11 1.5450042897697011e-02 -1.0557539868014379e-01 -1.1539698970598398e-01
12 -1.6489596487183067e+00 4.4215309402003902e+00 1.2758457405639687e+00
13 -1.3347630087741644e-01 3.6628535027741022e-02 -5.6719354360201611e-02
14 -4.1553429683017021e-02 -1.1540368471833904e-02 1.5467778628554116e-01
15 1.7577484997124486e+00 -4.3805631271698307e+00 -1.2463135748050815e+00
16 1.0396685127639720e-03 7.5317296594891241e-04 4.1929174719427522e-05
17 1.1778403499163440e-03 -7.1238534940179363e-05 3.8035065505820104e-03
18 -3.3315903593767938e-02 -2.1130575404712564e-02 1.4962162894607408e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.6557683033616174e-03 -9.8089615075005066e-04 -4.5869764220629459e-04
21 6.4422979624128504e-03 1.1468754348007283e-03 6.8525347184504311e-04
22 -7.5216145678175033e-04 1.1318396788931372e-04 -9.2834078115899713e-04
23 -5.5282158633416767e+00 3.3673691649465547e+00 -5.7944723588252627e+00
24 5.8210356160171492e-02 7.7028326238927844e-02 1.0457588922155134e-01
25 -1.4075156519665147e-01 1.8202248246942915e-02 -2.0674716237672891e-03
26 1.1435050190904306e-01 8.6100147519947955e-03 2.0290916519561587e-02
27 -3.8088820630046961e-02 -1.9961457381631506e-02 -9.8563367312923152e-03
...

121
unittest/force-styles/tests/atomic-pair-colloid_multi_tri.yaml
View File

@ -1,12 +1,15 @@
---
lammps_version: 21 Jul 2020
date_generated: Mon Aug 3 22:40:52 202
date_generated: Fri Aug 7 11:57:40 202
epsilon: 5e-13
prerequisites: ! |
pair colloid
pre_commands: ! |
variable units index lj
variable atom_style index atomic
post_commands: ! |
mass 1 9.0
mass 2 1.0
neighbor 1.0 multi
comm_style brick
comm_modify mode multi
@ -19,68 +22,68 @@ pair_coeff: ! |
2 2 10.0 1.0 0.0 0.0 2.5
extract: ! ""
natoms: 27
init_vdwl: -0.250265293649414
init_vdwl: -0.250320882681715
init_coul: 0
init_stress: ! |2-
1.6693554149342911e+00 1.1752557173260259e+01 1.0084916384631786e+00 -5.5164068587184580e+00 3.4448938909043375e-01 2.8524984331506325e+00
1.6687682615058241e+00 1.1747891073356824e+01 1.0083450567162384e+00 -5.5144333681333002e+00 3.4548099128804161e-01 2.8508471154796613e+00
init_forces: ! |2
1 7.2698902368313792e-02 7.6228314794484928e-02 1.0726717051410495e-01
2 -4.1487528743785314e-02 -9.8021436494101673e-03 1.5709289685835706e-01
3 -6.4514052467860251e-03 -1.0256365723614392e-03 -6.5497018575823654e-04
4 -1.5879716299045554e-02 -1.3313626328419375e-02 -1.1960220488237115e-02
5 -3.6475105049685402e-02 -2.2992583864092257e-02 1.6166097759476136e-02
6 5.4769180740244731e-04 -2.4426732691233955e-02 -3.6444902428179146e-02
7 1.2006858444817004e-03 8.6089326153587114e-04 4.8033026817420977e-05
8 6.2576206692511894e-03 1.1942031910109575e-03 8.5598267089392129e-04
9 1.8304752931477772e+00 -4.3380838069003751e+00 -1.2763141904695217e+00
10 1.0456768773213966e-01 -2.8797785965171449e-02 -1.0768387599117872e-02
11 7.4457152913473068e-04 5.3043528596594895e-05 1.1848179655097051e-04
12 1.8275344545112628e-04 4.1948901318387846e-04 9.8529270880004928e-05
13 -1.1259391902077470e-01 1.6396787556459836e-02 -9.0455786475091167e-04
14 -2.4411189619188755e-04 -5.8284002050384791e-05 1.8046434825690706e-05
15 -2.9421019290851403e-03 -7.4831176382853144e-04 -1.2326238660052020e-03
16 2.2461127490925330e+00 -1.1352527056039741e+00 2.3679980564815390e+00
17 -1.8842133286566962e-02 1.7091738246445273e-02 -1.4247568288438442e-02
18 1.0304762560041145e-01 7.8062023784307155e-03 1.9257495383127543e-02
19 -4.3351162132476757e-02 -2.2995928402219718e-02 -1.1392806036292676e-02
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.2501672350519284e-02 6.8582908647593982e-03 -6.8677246013076107e-02
22 -2.2945818607502835e+00 1.1603382884260451e+00 -2.3594708301087310e+00
23 2.5606115583434334e-02 -1.0270997331622238e-01 -1.1451933998746676e-01
24 -8.0648597970214874e-04 1.1943689027481719e-04 -9.8441554909815452e-04
25 6.1663239781204458e-04 -2.9236784591333030e-04 3.5387752678593654e-03
26 -1.7197958284241512e+00 4.3707144027711022e+00 1.2977194913045054e+00
27 -1.4110864281012736e-01 4.2418795982941855e-02 -6.2606997884263993e-02
run_vdwl: 0.0367090517931635
1 4.2502214537152139e-02 6.8580119035896794e-03 -6.8677100496229632e-02
2 -2.4409971231138493e-04 -5.8282143696384329e-05 1.8044935805594828e-05
3 -1.5879676810899381e-02 -1.3312213962585501e-02 -1.1959790230130862e-02
4 5.4778343227773751e-04 -2.4424984321614905e-02 -3.6445490162238962e-02
5 7.4458982140296185e-04 5.3043116325022762e-05 1.1848540218007482e-04
6 2.2466945989093152e+00 -1.1355562227620701e+00 2.3685617789648945e+00
7 -2.9420817632825631e-03 -7.4830177547552694e-04 -1.2326317038179192e-03
8 1.0456973881919868e-01 -2.8794438516942699e-02 -1.0767233651978712e-02
9 -1.8842568296904984e-02 1.7092235687070075e-02 -1.4248069827099497e-02
10 1.8278642893333129e-04 4.1941592242561638e-04 9.8624961124791906e-05
11 2.5605560774882380e-02 -1.0270993584698010e-01 -1.1452001476172038e-01
12 -1.7190177814448016e+00 4.3689733469586800e+00 1.2971816478431131e+00
13 -1.4111675314136379e-01 4.2420896699357342e-02 -6.2611143561578070e-02
14 -4.1487669742221053e-02 -9.8022075695912807e-03 1.5709450788901488e-01
15 1.8296904307488884e+00 -4.3363458812932363e+00 -1.2757754576583198e+00
16 1.2006895887469950e-03 8.6089636114693700e-04 4.8031493704482679e-05
17 6.1663662317927549e-04 -2.9235023723413927e-04 3.5387343984160471e-03
18 -3.6474070486619471e-02 -2.2991550342271464e-02 1.6164948715920183e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.4513560799741325e-03 -1.0256280577842611e-03 -6.5496554886125568e-04
21 6.2575756583543186e-03 1.1941945764063838e-03 8.5597303165342044e-04
22 -8.0641944111497544e-04 1.1943078425573959e-04 -9.8433130461350322e-04
23 -2.2951590782728784e+00 1.1606385208734213e+00 -2.3600356755964280e+00
24 7.2701432066201735e-02 7.6224922982964224e-02 1.0727212071804387e-01
25 -1.1259499193288525e-01 1.6395143300644142e-02 -9.0439962448757977e-04
26 1.0304882476745569e-01 7.8050496446502114e-03 1.9255643852370365e-02
27 -4.3346315050731436e-02 -2.2993111981455479e-02 -1.1392238078737104e-02
run_vdwl: 0.0367182211916219
run_coul: 0
run_stress: ! |2-
3.2547164843321998e+00 1.2822272288053016e+01 2.9381904014809583e+00 -6.5755478139674493e+00 2.3052669923975677e+00 1.5140903126232053e+00
3.2545234036053388e+00 1.2817594721674228e+01 2.9383804763809755e+00 -6.5738604535060468e+00 2.3065705273606301e+00 1.5121755374668360e+00
run_forces: ! |2
1 5.8196197291626575e-02 7.7034242946349710e-02 1.0456829776646283e-01
2 -4.1552988555246369e-02 -1.1540194679885474e-02 1.5467509805429086e-01
3 -6.6558201800471888e-03 -9.8090417095619570e-04 -4.5870088099424298e-04
4 -1.6184658518487716e-02 -1.2048418808706531e-02 -1.2190156970179285e-02
5 -3.3316803862839500e-02 -2.1131495092834723e-02 1.4963198384123161e-02
6 4.7038346343351856e-04 -3.0307837184768479e-02 -4.3608938525025828e-02
7 1.0396652127649044e-03 7.5317021954245870e-04 4.1930480907775929e-05
8 6.4423457305213851e-03 1.1468835414812505e-03 6.8526191116135909e-04
9 1.7585178109228026e+00 -4.3823551563818404e+00 -1.2468505588145364e+00
10 1.3392700382201883e-01 -2.9471624976996413e-02 -9.0915578208389941e-03
11 7.9839096035526651e-04 6.5881060273478051e-05 1.3840230606223500e-04
12 1.8540711168999979e-04 3.9935505493719320e-04 1.2187452986588761e-04
13 -1.4075050026169872e-01 1.8204883476184680e-02 -2.0674384501958849e-03
14 -2.4452271394244732e-04 -5.7393635529214925e-05 1.6908863929994199e-05
15 -3.4577717467029808e-03 -9.3350658134032805e-04 -1.5696623820874909e-03
16 5.4762334392971521e+00 -3.3427628201841029e+00 5.8008787111000384e+00
17 -1.7123087875076739e-02 1.5861563712562882e-02 -1.3049347364948762e-02
18 1.1434938302418579e-01 8.6112942994289246e-03 2.0292892978925857e-02
19 -3.8092989065812349e-02 -1.9963855292591714e-02 -9.8568020007968171e-03
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.1446743356056477e-02 8.4986882092497468e-03 -7.1565549648589286e-02
22 -5.5269240273286080e+00 3.3665609275390400e+00 -5.7932229791606122e+00
23 1.5450935830291862e-02 -1.0557581147017453e-01 -1.1539704280712934e-01
24 -7.5222270706101790e-04 1.1318969997713260e-04 -9.2841932954143699e-04
25 1.1778375109233756e-03 -7.1254208855309979e-05 3.8035457799329474e-03
26 -1.6497221843622341e+00 4.4233264004280688e+00 1.2763820593322264e+00
27 -1.3345796635606499e-01 3.6623792481486409e-02 -5.6711027332452466e-02
1 4.1447236896055638e-02 8.4984633363665851e-03 -7.1565468520811498e-02
2 -2.4451051586476380e-04 -5.7391827303752326e-05 1.6907428582411937e-05
3 -1.6184661894786685e-02 -1.2047019655368075e-02 -1.2189740041384702e-02
4 4.7049123492350683e-04 -3.0305522972901157e-02 -4.3609470856651830e-02
5 7.9841068285770012e-04 6.5880880038332612e-05 1.3840659114887713e-04
6 5.4775302106929926e+00 -3.3435741912382899e+00 5.8021270645808745e+00
7 -3.4577504657755053e-03 -9.3349532033408424e-04 -1.5696706937008277e-03
8 1.3392885595460410e-01 -2.9467417351771743e-02 -9.0902656642780511e-03
9 -1.7123469245219562e-02 1.5862013571050645e-02 -1.3049795876411906e-02
10 1.8544095678649792e-04 3.9928059116526535e-04 1.2197284521083560e-04
11 1.5450042897696470e-02 -1.0557539868014387e-01 -1.1539698970598380e-01
12 -1.6489596487183285e+00 4.4215309402004559e+00 1.2758457405639874e+00
13 -1.3347630087741638e-01 3.6628535027741001e-02 -5.6719354360201764e-02
14 -4.1553429683017180e-02 -1.1540368471833928e-02 1.5467778628554141e-01
15 1.7577484997124713e+00 -4.3805631271698946e+00 -1.2463135748050993e+00
16 1.0396685127639722e-03 7.5317296594891231e-04 4.1929174719427319e-05
17 1.1778403499163354e-03 -7.1238534940187549e-05 3.8035065505820221e-03
18 -3.3315903593768015e-02 -2.1130575404712623e-02 1.4962162894607469e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.6557683033616182e-03 -9.8089615075005370e-04 -4.5869764220629415e-04
21 6.4422979624128504e-03 1.1468754348007311e-03 6.8525347184504181e-04
22 -7.5216145678175326e-04 1.1318396788931422e-04 -9.2834078115900081e-04
23 -5.5282158633416776e+00 3.3673691649465551e+00 -5.7944723588252627e+00
24 5.8210356160171721e-02 7.7028326238927830e-02 1.0457588922155144e-01
25 -1.4075156519665108e-01 1.8202248246942852e-02 -2.0674716237672969e-03
26 1.1435050190904378e-01 8.6100147519948857e-03 2.0290916519561344e-02
27 -3.8088820630047288e-02 -1.9961457381631711e-02 -9.8563367312924106e-03
...

121
unittest/force-styles/tests/atomic-pair-colloid_tiled.yaml
View File

@ -1,12 +1,15 @@
---
lammps_version: 21 Jul 2020
date_generated: Mon Aug 3 16:06:24 202
date_generated: Fri Aug 7 11:57:40 202
epsilon: 5e-14
prerequisites: ! |
pair colloid
pre_commands: ! |
variable units index lj
variable atom_style index atomic
post_commands: ! |
mass 1 9.0
mass 2 1.0
neighbor 1.0 multi
comm_style tiled
comm_modify mode multi
@ -18,68 +21,68 @@ pair_coeff: ! |
2 2 10.0 1.0 0.0 0.0 2.5
extract: ! ""
natoms: 27
init_vdwl: -0.250265293649414
init_vdwl: -0.250320882681727
init_coul: 0
init_stress: ! |2-
1.6693554149342913e+00 1.1752557173260270e+01 1.0084916384631617e+00 -5.5164068587184527e+00 3.4448938909043431e-01 2.8524984331506258e+00
1.6687682615057460e+00 1.1747891073356403e+01 1.0083450567161796e+00 -5.5144333681331261e+00 3.4548099128807747e-01 2.8508471154795298e+00
init_forces: ! |2
1 7.2698902368314056e-02 7.6228314794484775e-02 1.0726717051410484e-01
2 -4.1487528743785314e-02 -9.8021436494101673e-03 1.5709289685835706e-01
3 -6.4514052467860251e-03 -1.0256365723614392e-03 -6.5497018575823654e-04
4 -1.5879716299045554e-02 -1.3313626328419375e-02 -1.1960220488237115e-02
5 -3.6475105049685402e-02 -2.2992583864092257e-02 1.6166097759476136e-02
6 5.4769180740248721e-04 -2.4426732691233789e-02 -3.6444902428178980e-02
7 1.2006858444817004e-03 8.6089326153587114e-04 4.8033026817420977e-05
8 6.2576206692511894e-03 1.1942031910109577e-03 8.5598267089392085e-04
9 1.8304752931477772e+00 -4.3380838069003760e+00 -1.2763141904695217e+00
10 1.0456768773213965e-01 -2.8797785965171453e-02 -1.0768387599117879e-02
11 7.4457152913473112e-04 5.3043528596594623e-05 1.1848179655097060e-04
12 1.8275344545112628e-04 4.1948901318387846e-04 9.8529270880004928e-05
13 -1.1259391902077470e-01 1.6396787556459839e-02 -9.0455786475091124e-04
14 -2.4411189619188785e-04 -5.8284002050384879e-05 1.8046434825690733e-05
15 -2.9421019290851403e-03 -7.4831176382853144e-04 -1.2326238660052020e-03
16 2.2461127490925330e+00 -1.1352527056039741e+00 2.3679980564815386e+00
17 -1.8842133286566962e-02 1.7091738246445273e-02 -1.4247568288438442e-02
18 1.0304762560041125e-01 7.8062023784307467e-03 1.9257495383127633e-02
19 -4.3351162132476757e-02 -2.2995928402219718e-02 -1.1392806036292676e-02
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.2501672350519284e-02 6.8582908647593938e-03 -6.8677246013076107e-02
22 -2.2945818607502839e+00 1.1603382884260454e+00 -2.3594708301087310e+00
23 2.5606115583434265e-02 -1.0270997331622232e-01 -1.1451933998746673e-01
24 -8.0648597970214809e-04 1.1943689027481706e-04 -9.8441554909815322e-04
25 6.1663239781204458e-04 -2.9236784591333030e-04 3.5387752678593654e-03
26 -1.7197958284241512e+00 4.3707144027711022e+00 1.2977194913045054e+00
27 -1.4110864281012736e-01 4.2418795982941855e-02 -6.2606997884263993e-02
run_vdwl: 0.0367090517931737
1 4.2502214537152215e-02 6.8580119035896499e-03 -6.8677100496229715e-02
2 -2.4409971231138453e-04 -5.8282143696384295e-05 1.8044935805594862e-05
3 -1.5879676810899343e-02 -1.3312213962585376e-02 -1.1959790230130810e-02
4 5.4778343227770802e-04 -2.4424984321614867e-02 -3.6445490162239011e-02
5 7.4458982140296174e-04 5.3043116325023304e-05 1.1848540218007477e-04
6 2.2466945989092868e+00 -1.1355562227620619e+00 2.3685617789648643e+00
7 -2.9420817632825613e-03 -7.4830177547552715e-04 -1.2326317038179183e-03
8 1.0456973881919887e-01 -2.8794438516942737e-02 -1.0767233651978719e-02
9 -1.8842568296905022e-02 1.7092235687070117e-02 -1.4248069827099535e-02
10 1.8278642893333091e-04 4.1941592242561833e-04 9.8624961124789412e-05
11 2.5605560774882609e-02 -1.0270993584697993e-01 -1.1452001476172027e-01
12 -1.7190177814447429e+00 4.3689733469585290e+00 1.2971816478430653e+00
13 -1.4111675314136360e-01 4.2420896699357280e-02 -6.2611143561577903e-02
14 -4.1487669742220991e-02 -9.8022075695912841e-03 1.5709450788901494e-01
15 1.8296904307488300e+00 -4.3363458812930853e+00 -1.2757754576582720e+00
16 1.2006895887469950e-03 8.6089636114693700e-04 4.8031493704482794e-05
17 6.1663662317927744e-04 -2.9235023723413574e-04 3.5387343984160432e-03
18 -3.6474070486619430e-02 -2.2991550342271461e-02 1.6164948715920169e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.4513560799741325e-03 -1.0256280577842596e-03 -6.5496554886125568e-04
21 6.2575756583543186e-03 1.1941945764063822e-03 8.5597303165342044e-04
22 -8.0641944111498032e-04 1.1943078425573997e-04 -9.8433130461350886e-04
23 -2.2951590782728504e+00 1.1606385208734131e+00 -2.3600356755963983e+00
24 7.2701432066201888e-02 7.6224922982964599e-02 1.0727212071804376e-01
25 -1.1259499193288544e-01 1.6395143300644177e-02 -9.0439962448757879e-04
26 1.0304882476745540e-01 7.8050496446500431e-03 1.9255643852370385e-02
27 -4.3346315050731818e-02 -2.2993111981455712e-02 -1.1392238078737215e-02
run_vdwl: 0.0367182211916192
run_coul: 0
run_stress: ! |2-
3.2547164843322691e+00 1.2822272288053391e+01 2.9381904014809663e+00 -6.5755478139675931e+00 2.3052669923975411e+00 1.5140903126232879e+00
3.2545234036053201e+00 1.2817594721674030e+01 2.9383804763809369e+00 -6.5738604535059864e+00 2.3065705273606469e+00 1.5121755374667567e+00
run_forces: ! |2
1 5.8196197291625715e-02 7.7034242946349946e-02 1.0456829776646254e-01
2 -4.1552988555246327e-02 -1.1540194679885485e-02 1.5467509805429100e-01
3 -6.6558201800471870e-03 -9.8090417095619679e-04 -4.5870088099424396e-04
4 -1.6184658518487705e-02 -1.2048418808706472e-02 -1.2190156970179262e-02
5 -3.3316803862839438e-02 -2.1131495092834646e-02 1.4963198384123119e-02
6 4.7038346343352116e-04 -3.0307837184768524e-02 -4.3608938525025766e-02
7 1.0396652127649081e-03 7.5317021954246196e-04 4.1930480907776011e-05
8 6.4423457305213877e-03 1.1468835414812518e-03 6.8526191116136061e-04
9 1.7585178109228492e+00 -4.3823551563819700e+00 -1.2468505588145766e+00
10 1.3392700382201897e-01 -2.9471624976996433e-02 -9.0915578208389906e-03
11 7.9839096035526857e-04 6.5881060273478105e-05 1.3840230606223543e-04
12 1.8540711169000120e-04 3.9935505493718929e-04 1.2187452986589205e-04
13 -1.4075050026169886e-01 1.8204883476184704e-02 -2.0674384501958810e-03
14 -2.4452271394244688e-04 -5.7393635529214755e-05 1.6908863929994165e-05
15 -3.4577717467029726e-03 -9.3350658134032599e-04 -1.5696623820874889e-03
16 5.4762334392971725e+00 -3.3427628201841184e+00 5.8008787111000641e+00
17 -1.7123087875076635e-02 1.5861563712562771e-02 -1.3049347364948664e-02
18 1.1434938302418582e-01 8.6112942994291432e-03 2.0292892978926169e-02
19 -3.8092989065812294e-02 -1.9963855292591683e-02 -9.8568020007968223e-03
20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
21 4.1446743356056415e-02 8.4986882092498717e-03 -7.1565549648589466e-02
22 -5.5269240273286284e+00 3.3665609275390551e+00 -5.7932229791606371e+00
23 1.5450935830291834e-02 -1.0557581147017484e-01 -1.1539704280712981e-01
24 -7.5222270706101302e-04 1.1318969997713211e-04 -9.2841932954143037e-04
25 1.1778375109233736e-03 -7.1254208855308434e-05 3.8035457799329383e-03
26 -1.6497221843622811e+00 4.4233264004281976e+00 1.2763820593322663e+00
27 -1.3345796635606427e-01 3.6623792481486235e-02 -5.6711027332452084e-02
1 4.1447236896055561e-02 8.4984633363665869e-03 -7.1565468520811457e-02
2 -2.4451051586476304e-04 -5.7391827303752204e-05 1.6907428582411947e-05
3 -1.6184661894786571e-02 -1.2047019655367870e-02 -1.2189740041384602e-02
4 4.7049123492354846e-04 -3.0305522972900800e-02 -4.3609470856651497e-02
5 7.9841068285770055e-04 6.5880880038332991e-05 1.3840659114887751e-04
6 5.4775302106929935e+00 -3.3435741912382895e+00 5.8021270645808736e+00
7 -3.4577504657755018e-03 -9.3349532033408424e-04 -1.5696706937008275e-03
8 1.3392885595460452e-01 -2.9467417351771795e-02 -9.0902656642780702e-03
9 -1.7123469245219562e-02 1.5862013571050645e-02 -1.3049795876411903e-02
10 1.8544095678649998e-04 3.9928059116526140e-04 1.2197284521083842e-04
11 1.5450042897697011e-02 -1.0557539868014379e-01 -1.1539698970598398e-01
12 -1.6489596487183067e+00 4.4215309402003902e+00 1.2758457405639687e+00
13 -1.3347630087741644e-01 3.6628535027741022e-02 -5.6719354360201611e-02
14 -4.1553429683017021e-02 -1.1540368471833904e-02 1.5467778628554116e-01
15 1.7577484997124486e+00 -4.3805631271698307e+00 -1.2463135748050815e+00
16 1.0396685127639720e-03 7.5317296594891241e-04 4.1929174719427522e-05
17 1.1778403499163440e-03 -7.1238534940179363e-05 3.8035065505820104e-03
18 -3.3315903593767938e-02 -2.1130575404712564e-02 1.4962162894607408e-02
19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
20 -6.6557683033616174e-03 -9.8089615075005066e-04 -4.5869764220629459e-04
21 6.4422979624128504e-03 1.1468754348007283e-03 6.8525347184504311e-04
22 -7.5216145678175033e-04 1.1318396788931372e-04 -9.2834078115899713e-04
23 -5.5282158633416767e+00 3.3673691649465547e+00 -5.7944723588252627e+00
24 5.8210356160171492e-02 7.7028326238927844e-02 1.0457588922155134e-01
25 -1.4075156519665147e-01 1.8202248246942915e-02 -2.0674716237672891e-03
26 1.1435050190904306e-01 8.6100147519947955e-03 2.0290916519561587e-02
27 -3.8088820630046961e-02 -1.9961457381631506e-02 -9.8563367312923152e-03
...

Some files were not shown because too many files have changed in this diff Show More