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

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

Browse Source 
Collected small changes
This commit is contained in:
Axel Kohlmeyer
2022-03-10 23:42:04 -05:00
committed by GitHub
parent c6b4c418e7 3a1921b3ed
commit 9b989b1860
31 changed files with 859 additions and 900 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
cmake/Modules/Packages/GPU.cmake
View File

@ -30,7 +30,15 @@ file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/[^.]*.cpp)
file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu)
if(GPU_API STREQUAL "CUDA")
find_package(CUDA REQUIRED)
find_package(CUDA QUIET)
# augment search path for CUDA toolkit libraries to include the stub versions. Needed to find libcuda.so on machines without a CUDA driver installation
if(CUDA_FOUND)
set(CMAKE_LIBRARY_PATH "${CUDA_TOOLKIT_ROOT_DIR}/lib64/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib64;${CUDA_TOOLKIT_ROOT_DIR}/lib;${CMAKE_LIBRARY_PATH}")
find_package(CUDA REQUIRED)
else()
message(FATAL_ERROR "CUDA Toolkit not found")
endif()
find_program(BIN2C bin2c)
if(NOT BIN2C)
message(FATAL_ERROR "Could not find bin2c, use -DBIN2C=/path/to/bin2c to help cmake finding it.")

18
doc/src/Build_settings.rst
View File

@ -4,15 +4,15 @@ Optional build settings
LAMMPS can be built with several optional settings. Each sub-section
explain how to do this for building both with CMake and make.
* :ref:`C++11 standard compliance <cxx11>` when building all of LAMMPS
* :ref:`FFT library <fft>` for use with the :doc:`kspace_style pppm <kspace_style>` command
* :ref:`Size of LAMMPS integer types <size>`
* :ref:`Read or write compressed files <gzip>`
* :ref:`Output of JPG and PNG files <graphics>` via the :doc:`dump image <dump_image>` command
* :ref:`Output of movie files <graphics>` via the :doc:`dump_movie <dump_image>` command
* :ref:`Memory allocation alignment <align>`
* :ref:`Workaround for long long integers <longlong>`
* :ref:`Error handling exceptions <exceptions>` when using LAMMPS as a library
* `C++11 standard compliance`_ when building all of LAMMPS
* `FFT library`_ for use with the :doc:`kspace_style pppm <kspace_style>` command
* `Size of LAMMPS integer types and size limits`_
* `Read or write compressed files`_
* `Output of JPG, PNG, and move files` via the :doc:`dump image <dump_image>` or :doc:`dump movie <dump_image>` commands
* `Memory allocation alignment`_
* `Workaround for long long integers`_
* `Exception handling when using LAMMPS as a library`_ to capture errors
* `Trigger selected floating-point exceptions`_
----------

6
doc/src/Commands_compute.rst
View File

@ -142,11 +142,11 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`sph/t/atom <compute_sph_t_atom>`
* :doc:`spin <compute_spin>`
* :doc:`stress/atom <compute_stress_atom>`
* :doc:`stress/cartesian <compute_stress_cartesian>`
* :doc:`stress/cylinder <compute_stress_cartesian>`
* :doc:`stress/cartesian <compute_stress_profile>`
* :doc:`stress/cylinder <compute_stress_profile>`
* :doc:`stress/mop <compute_stress_mop>`
* :doc:`stress/mop/profile <compute_stress_mop>`
* :doc:`stress/spherical <compute_stress_cartesian>`
* :doc:`stress/spherical <compute_stress_profile>`
* :doc:`stress/tally <compute_tally>`
* :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>`
* :doc:`temp (k) <compute_temp>`

23
doc/src/Developer_notes.rst
View File

@ -25,18 +25,17 @@ following benefits:
- transparent support for translating unsupported UTF-8 characters to their ASCII equivalents
(the text-to-value conversion functions **only** accept ASCII characters)
In most cases (e.g. potential files) the same data is needed on all
MPI ranks. Then it is best to do the reading and parsing only on MPI
rank 0, and communicate the data later with one or more
``MPI_Bcast()`` calls. For reading generic text and potential
parameter files the custom classes :cpp:class:`TextFileReader
<LAMMPS_NS::TextFileReader>` and :cpp:class:`PotentialFileReader
<LAMMPS_NS::PotentialFileReader>` are available. They allow reading
the file as individual lines for which they can return a tokenizer
class (see below) for parsing the line. Or they can return blocks of
numbers as a vector directly. The documentation on `File reader
classes <file-reader-classes>`_ contains an example for a typical
case.
In most cases (e.g. potential files) the same data is needed on all MPI
ranks. Then it is best to do the reading and parsing only on MPI rank
0, and communicate the data later with one or more ``MPI_Bcast()``
calls. For reading generic text and potential parameter files the
custom classes :cpp:class:`TextFileReader <LAMMPS_NS::TextFileReader>`
and :cpp:class:`PotentialFileReader <LAMMPS_NS::PotentialFileReader>`
are available. They allow reading the file as individual lines for which
they can return a tokenizer class (see below) for parsing the line. Or
they can return blocks of numbers as a vector directly. The
documentation on :ref:`File reader classes <file-reader-classes>`
contains an example for a typical case.
When reading per-atom data, the data on each line of the file usually
needs to include an atom ID so it can be associated with a particular

6
doc/src/Developer_org.rst
View File

@ -252,12 +252,6 @@ follows:
- The Timer class logs timing information, output at the end
of a run.
.. TODO section on "Spatial decomposition and parallel operations"
.. diagram of 3d processor grid, brick vs. tiled. local vs. ghost
.. atoms, 6-way communication with pack/unpack functions,
.. PBC as part of the communication, forward and reverse communication
.. rendezvous communication, ring communication.
.. TODO section on "Fixes, Computes, and Variables"
.. how and when data is computed and provided and how it is
.. referenced. flags in Fix/Compute/Variable classes tell

2
doc/src/Developer_utils.rst
View File

@ -396,7 +396,7 @@ A typical code segment would look like this:
----------
.. file-reader-classes:
.. _file-reader-classes:
File reader classes
-------------------

2
doc/src/Packages_details.rst
View File

@ -9,7 +9,7 @@ gives links to documentation, example scripts, and pictures/movies (if
available) that illustrate use of the package.
The majority of packages can be included in a LAMMPS build with a
single setting (``-D PGK_<NAME>=on`` for CMake) or command
single setting (``-D PKG_<NAME>=on`` for CMake) or command
(``make yes-<name>`` for make). See the :doc:`Build package <Build_package>`
page for more info. A few packages may require additional steps;
this is indicated in the descriptions below. The :doc:`Build extras <Build_extras>`

6
doc/src/compute.rst
View File

@ -288,11 +288,11 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`sph/t/atom <compute_sph_t_atom>` - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms
* :doc:`spin <compute_spin>` - magnetic quantities for a system of atoms having spins
* :doc:`stress/atom <compute_stress_atom>` - stress tensor for each atom
* :doc:`stress/cartesian <compute_stress_cartesian>` - stress tensor in cartesian coordinates
* :doc:`stress/cylinder <compute_stress_cartesian>` - stress tensor in cylindrical coordinates
* :doc:`stress/cartesian <compute_stress_profile>` - stress tensor in cartesian coordinates
* :doc:`stress/cylinder <compute_stress_profile>` - stress tensor in cylindrical coordinates
* :doc:`stress/mop <compute_stress_mop>` - normal components of the local stress tensor using the method of planes
* :doc:`stress/mop/profile <compute_stress_mop>` - profile of the normal components of the local stress tensor using the method of planes
* :doc:`stress/spherical <compute_stress_cartesian>` - stress tensor in spherical coordinates
* :doc:`stress/spherical <compute_stress_profile>` - stress tensor in spherical coordinates
* :doc:`stress/tally <compute_tally>` - stress between two groups of atoms via the tally callback mechanism
* :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>` - per-atom chemical concentration of a specified species for each tDPD particle
* :doc:`temp <compute_temp>` - temperature of group of atoms

9
doc/src/compute_momentum.rst
View File

@ -23,11 +23,10 @@ Examples
Description
"""""""""""
Define a computation that calculates the translational momentum
of a group of particles.
The momentum of each particles is computed as m v, where m and v are
the mass and velocity of the particle.
Define a computation that calculates the translational momentum *p*
of a group of particles. It is computed as the sum :math:`\vec{p} = \sum_i m_i \cdot \vec{v}_i`
over all particles in the compute group, where *m* and *v* are
the mass and velocity vector of the particle, respectively.
Output info
"""""""""""

4
doc/src/compute_stress_mop.rst
View File

@ -73,7 +73,7 @@ command, since those are contributions to the global system pressure.
NOTE 3: The local stress profile generated by compute *stress/mop/profile*
is similar to that obtained by compute
:doc:`stress/cartesian <compute_stress_cartesian>`.
:doc:`stress/cartesian <compute_stress_profile>`.
A key difference
is that compute *stress/mop/profile* considers particles
crossing a set of planes,
@ -122,7 +122,7 @@ intra-molecular interactions, and long range (kspace) interactions.
Related commands
""""""""""""""""
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`, :doc:`compute stress/cartesian <compute_stress_cartesian>`, :doc:`compute stress/cylinder <compute_stress_cartesian>`, :doc:`compute stress/spherical <compute_stress_cartesian>`
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`, :doc:`compute stress/cartesian <compute_stress_profile>`, :doc:`compute stress/cylinder <compute_stress_profile>`, :doc:`compute stress/spherical <compute_stress_profile>`
Default
"""""""

0
doc/src/compute_stress_cartesian.rst → doc/src/compute_stress_profile.rst
View File

60
doc/src/compute_tally.rst
View File

@ -31,7 +31,7 @@ Syntax
compute ID group-ID style group2-ID
* ID, group-ID are documented in :doc:`compute <compute>` command
* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or * or *pe/tally* or *pe/mol/tally* or *stress/tally*
* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or *pe/tally* or *pe/mol/tally* or *stress/tally*
* group2-ID = group ID of second (or same) group
Examples
@ -61,7 +61,7 @@ mechanism. Compute *pe/mol/tally* is one such style, that can
- through using this mechanism - separately tally intermolecular
and intramolecular energies. Something that would otherwise be
impossible without integrating this as a core functionality into
the based classes of LAMMPS.
the base classes of LAMMPS.
----------
@ -148,30 +148,38 @@ pairwise property computations.
Output info
"""""""""""
Compute *pe/tally* calculates a global scalar (the energy) and a per
atom scalar (the contributions of the single atom to the global
scalar). Compute *pe/mol/tally* calculates a global 4-element vector
containing (in this order): *evdwl* and *ecoul* for intramolecular pairs
and *evdwl* and *ecoul* for intermolecular pairs. Since molecules are
identified by their molecule IDs, the partitioning does not have to be
related to molecules, but the energies are tallied into the respective
slots depending on whether the molecule IDs of a pair are the same or
different. Compute *force/tally* calculates a global scalar (the force
magnitude) and a per atom 3-element vector (force contribution from
each atom). Compute *stress/tally* calculates a global scalar
(average of the diagonal elements of the stress tensor) and a per atom
vector (the 6 elements of stress tensor contributions from the
individual atom). As in :doc:`compute heat/flux <compute_heat_flux>`,
compute *heat/flux/tally* calculates a global vector of length 6,
where the first 3 components are the :math:`x`, :math:`y`, :math:`z`
components of the full heat flow vector,
and the next 3 components are the corresponding components
of just the convective portion of the flow, i.e. the
first term in the equation for :math:`\mathbf{Q}`.
Compute *heat/flux/virial/tally* calculates a global scalar (heat flow)
and a per atom 3-element vector
(contribution to the force acting over atoms in the first group
from individual atoms in both groups).
- Compute *pe/tally* calculates a global scalar (the energy) and a per
atom scalar (the contributions of the single atom to the global
scalar).
- Compute *pe/mol/tally* calculates a global 4-element vector containing
(in this order): *evdwl* and *ecoul* for intramolecular pairs and
*evdwl* and *ecoul* for intermolecular pairs. Since molecules are
identified by their molecule IDs, the partitioning does not have to be
related to molecules, but the energies are tallied into the respective
slots depending on whether the molecule IDs of a pair are the same or
different.
- Compute *force/tally* calculates a global scalar (the force magnitude)
and a per atom 3-element vector (force contribution from each atom).
- Compute *stress/tally* calculates a global scalar
(average of the diagonal elements of the stress tensor) and a per atom
vector (the 6 elements of stress tensor contributions from the
individual atom).
- As in :doc:`compute heat/flux <compute_heat_flux>`,
compute *heat/flux/tally* calculates a global vector of length 6,
where the first 3 components are the :math:`x`, :math:`y`, :math:`z`
components of the full heat flow vector,
and the next 3 components are the corresponding components
of just the convective portion of the flow, i.e. the
first term in the equation for :math:`\mathbf{Q}`.
- Compute *heat/flux/virial/tally* calculates a global scalar (heat flow)
and a per atom 3-element vector
(contribution to the force acting over atoms in the first group
from individual atoms in both groups).
Both the scalar and vector values calculated by this compute are
"extensive".

13
src/.gitignore vendored
View File

@ -1,6 +1,7 @@
/Makefile.package
/Makefile.package.settings
/MAKE/MINE
/ADIOS/Makefile.lammps
/Make.py.last
/lmp_*
@ -267,6 +268,8 @@
/fix_drag.h
/fix_numdiff.cpp
/fix_numdiff.h
/fix_numdiff_virial.cpp
/fix_numdiff_virial.h
/fix_spring_rg.cpp
/fix_spring_rg.h
/fix_temp_csld.cpp
@ -934,6 +937,8 @@
/msm_cg.h
/neb.cpp
/neb.h
/netcdf_units.cpp
/netcdf_units.h
/pair_adp.cpp
/pair_adp.h
/pair_agni.cpp
@ -1069,6 +1074,8 @@
/pair_hdnnp.h
/pair_ilp_graphene_hbn.cpp
/pair_ilp_graphene_hbn.h
/pair_ilp_tmd.cpp
/pair_ilp_tmd.h
/pair_kolmogorov_crespi_full.cpp
/pair_kolmogorov_crespi_full.h
/pair_kolmogorov_crespi_z.cpp
@ -1175,8 +1182,8 @@
/pair_nm_cut_coul_cut.h
/pair_nm_cut_coul_long.cpp
/pair_nm_cut_coul_long.h
/pait_nm_cut_split.cpp
/pait_nm_cut_split.h
/pair_nm_cut_split.cpp
/pair_nm_cut_split.h
/pair_oxdna_*.cpp
/pair_oxdna_*.h
/pair_oxdna2_*.cpp
@ -1202,6 +1209,8 @@
/pair_rebo.h
/pair_resquared.cpp
/pair_resquared.h
/pair_saip_metal.cpp
/pair_saip_metal.h
/pair_sdpd_taitwater_isothermal.cpp
/pair_sdpd_taitwater_isothermal.h
/pair_sph_heatconduction.cpp

391
src/CLASS2/pair_lj_class2.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -65,13 +64,13 @@ PairLJClass2::~PairLJClass2()
void PairLJClass2::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -104,34 +103,32 @@ void PairLJClass2::compute(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fpair = factor_lj * forcelj * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype];
evdwl *= factor_lj;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
@ -144,10 +141,10 @@ void PairLJClass2::compute(int eflag, int vflag)
void PairLJClass2::compute_inner()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x;
double **f = atom->f;
@ -165,8 +162,8 @@ void PairLJClass2::compute_inner()
double cut_out_off = cut_respa[1];
double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms
@ -187,28 +184,28 @@ void PairLJClass2::compute_inner()
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
jtype = type[j];
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fpair = factor_lj * forcelj * r2inv;
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw);
}
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
}
@ -219,10 +216,10 @@ void PairLJClass2::compute_inner()
void PairLJClass2::compute_middle()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x;
double **f = atom->f;
@ -243,10 +240,10 @@ void PairLJClass2::compute_middle()
double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on * cut_in_on;
double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms
@ -267,32 +264,32 @@ void PairLJClass2::compute_middle()
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
jtype = type[j];
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fpair = factor_lj * forcelj * r2inv;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
}
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0);
}
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
}
@ -303,13 +300,13 @@ void PairLJClass2::compute_middle()
void PairLJClass2::compute_outer(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -327,8 +324,8 @@ void PairLJClass2::compute_outer(int eflag, int vflag)
double cut_in_on = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on * cut_in_on;
// loop over neighbors of my atoms
@ -349,55 +346,54 @@ void PairLJClass2::compute_outer(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
if (rsq > cut_in_off_sq) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fpair = factor_lj * forcelj * r2inv;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
}
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
if (eflag) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype];
evdwl *= factor_lj;
}
if (vflag) {
if (rsq <= cut_in_off_sq) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fpair = factor_lj * forcelj * r2inv;
} else if (rsq < cut_in_on_sq)
fpair = factor_lj*forcelj*r2inv;
fpair = factor_lj * forcelj * r2inv;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
@ -409,23 +405,21 @@ void PairLJClass2::compute_outer(int eflag, int vflag)
void PairLJClass2::allocate()
{
allocated = 1;
int n = atom->ntypes;
const int np1 = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(cut, np1, np1, "pair:cut");
memory->create(epsilon, np1, np1, "pair:epsilon");
memory->create(sigma, np1, np1, "pair:sigma");
memory->create(lj1, np1, np1, "pair:lj1");
memory->create(lj2, np1, np1, "pair:lj2");
memory->create(lj3, np1, np1, "pair:lj3");
memory->create(lj4, np1, np1, "pair:lj4");
memory->create(offset, np1, np1, "pair:offset");
}
/* ----------------------------------------------------------------------
@ -434,14 +428,14 @@ void PairLJClass2::allocate()
void PairLJClass2::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
if (narg != 1) error->all(FLERR, "Illegal pair_style command");
cut_global = utils::numeric(FLERR,arg[0],false,lmp);
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
@ -454,22 +448,22 @@ void PairLJClass2::settings(int narg, char **arg)
void PairLJClass2::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
double cut_one = cut_global;
if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp);
if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
@ -478,7 +472,7 @@ void PairLJClass2::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -492,12 +486,12 @@ void PairLJClass2::init_style()
int irequest;
int respa = 0;
if (update->whichflag == 1 && utils::strmatch(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;
}
irequest = neighbor->request(this,instance_me);
irequest = neighbor->request(this, instance_me);
if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1;
@ -507,10 +501,11 @@ void PairLJClass2::init_style()
// set rRESPA cutoffs
if (utils::strmatch(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 = nullptr;
else
cut_respa = nullptr;
}
/* ----------------------------------------------------------------------
@ -523,23 +518,22 @@ double PairLJClass2::init_one(int i, int j)
// mix distance via user-defined rule
if (setflag[i][j] == 0) {
epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) *
pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) /
(pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0));
sigma[i][j] =
pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) *
pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0));
sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0);
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
}
lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0);
lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0);
lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
if (offset_flag && (cut[i][j] > 0.0)) {
double ratio = sigma[i][j] / cut[i][j];
offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0));
} else offset[i][j] = 0.0;
offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0));
} else
offset[i][j] = 0.0;
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
@ -550,7 +544,7 @@ double PairLJClass2::init_one(int i, int j)
// check interior rRESPA cutoff
if (cut_respa && cut[i][j] < cut_respa[3])
error->all(FLERR,"Pair cutoff < Respa interior cutoff");
error->all(FLERR, "Pair cutoff < Respa interior cutoff");
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
@ -559,22 +553,21 @@ double PairLJClass2::init_one(int i, int j)
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
double count[2], all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world);
double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3;
etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 2.0*rc3) / rc6;
double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j];
double sig6 = sig3 * sig3;
double rc3 = cut[i][j] * cut[i][j] * cut[i][j];
double rc6 = rc3 * rc3;
etail_ij =
2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6);
ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6;
}
return cut[i][j];
@ -588,14 +581,14 @@ void PairLJClass2::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&sigma[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
}
}
}
@ -609,21 +602,21 @@ void PairLJClass2::read_restart(FILE *fp)
read_restart_settings(fp);
allocate();
int i,j;
int i, j;
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,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
@ -634,10 +627,10 @@ void PairLJClass2::read_restart(FILE *fp)
void PairLJClass2::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&tail_flag,sizeof(int),1,fp);
fwrite(&cut_global, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -648,26 +641,24 @@ void PairLJClass2::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairLJClass2::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]);
}
/* ----------------------------------------------------------------------
@ -678,27 +669,25 @@ void PairLJClass2::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]);
fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj,
double &fforce)
double /*factor_coul*/, double factor_lj, double &fforce)
{
double r2inv,rinv,r3inv,r6inv,forcelj,philj;
double r2inv, rinv, r3inv, r6inv, forcelj, philj;
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
fforce = factor_lj*forcelj*r2inv;
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
fforce = factor_lj * forcelj * r2inv;
philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
return factor_lj*philj;
philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype];
return factor_lj * philj;
}
/* ---------------------------------------------------------------------- */
@ -706,7 +695,7 @@ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double r
void *PairLJClass2::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma;
if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
if (strcmp(str, "sigma") == 0) return (void *) sigma;
return nullptr;
}

297
src/CLASS2/pair_lj_class2_coul_cut.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -14,17 +13,17 @@
#include "pair_lj_class2_coul_cut.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -65,14 +64,14 @@ PairLJClass2CoulCut::~PairLJClass2CoulCut()
void PairLJClass2CoulCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
double factor_coul,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj;
double factor_coul, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = ecoul = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -110,47 +109,49 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
else forcecoul = 0.0;
forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv);
else
forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
} else
forcelj = 0.0;
fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
fpair = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
if (rsq < cut_coulsq[itype][jtype])
ecoul = factor_coul * qqrd2e * qtmp*q[j]*sqrt(r2inv);
else ecoul = 0.0;
ecoul = factor_coul * qqrd2e * qtmp * q[j] * sqrt(r2inv);
else
ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
} else
evdwl = 0.0;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz);
}
}
}
@ -165,26 +166,25 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag)
void PairLJClass2CoulCut::allocate()
{
allocated = 1;
int n = atom->ntypes;
const int np1 = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
memory->create(cut_coul,n+1,n+1,"pair:cut_coul");
memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(cut_lj, np1, np1, "pair:cut_lj");
memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq");
memory->create(cut_coul, np1, np1, "pair:cut_coul");
memory->create(cut_coulsq, np1, np1, "pair:cut_coulsq");
memory->create(epsilon, np1, np1, "pair:epsilon");
memory->create(sigma, np1, np1, "pair:sigma");
memory->create(lj1, np1, np1, "pair:lj1");
memory->create(lj2, np1, np1, "pair:lj2");
memory->create(lj3, np1, np1, "pair:lj3");
memory->create(lj4, np1, np1, "pair:lj4");
memory->create(offset, np1, np1, "pair:offset");
}
/* ----------------------------------------------------------------------
@ -193,16 +193,18 @@ void PairLJClass2CoulCut::allocate()
void PairLJClass2CoulCut::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = utils::numeric(FLERR,arg[1],false,lmp);
cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 1)
cut_coul_global = cut_lj_global;
else
cut_coul_global = utils::numeric(FLERR, arg[1], false, lmp);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) {
@ -218,26 +220,25 @@ void PairLJClass2CoulCut::settings(int narg, char **arg)
void PairLJClass2CoulCut::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6)
error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 4 || narg > 6) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
double cut_lj_one = cut_lj_global;
double cut_coul_one = cut_coul_global;
if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp);
if (narg == 6) cut_coul_one = utils::numeric(FLERR,arg[5],false,lmp);
if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR, arg[4], false, lmp);
if (narg == 6) cut_coul_one = utils::numeric(FLERR, arg[5], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut_lj[i][j] = cut_lj_one;
@ -247,7 +248,7 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -256,10 +257,9 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg)
void PairLJClass2CoulCut::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/class2/coul/cut requires atom attribute q");
if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/cut requires atom attribute q");
neighbor->request(this,instance_me);
neighbor->request(this, instance_me);
}
/* ----------------------------------------------------------------------
@ -272,28 +272,27 @@ double PairLJClass2CoulCut::init_one(int i, int j)
// mix distance via user-defined rule
if (setflag[i][j] == 0) {
epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) *
pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) /
(pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0));
sigma[i][j] =
pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0);
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]);
epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) *
pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0));
sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0);
cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]);
cut_coul[i][j] = mix_distance(cut_coul[i][i], cut_coul[j][j]);
}
double cut = MAX(cut_lj[i][j],cut_coul[i][j]);
double cut = MAX(cut_lj[i][j], cut_coul[i][j]);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j];
lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0);
lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0);
lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
if (offset_flag && (cut_lj[i][j] > 0.0)) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0));
} else offset[i][j] = 0.0;
offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0));
} else
offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[i][j];
cut_coulsq[j][i] = cut_coulsq[i][j];
@ -310,22 +309,21 @@ double PairLJClass2CoulCut::init_one(int i, int j)
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
double count[2], all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world);
double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3;
etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 2.0*rc3) / rc6;
double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j];
double sig6 = sig3 * sig3;
double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j];
double rc6 = rc3 * rc3;
etail_ij =
2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6);
ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6;
}
return cut;
@ -339,15 +337,15 @@ void PairLJClass2CoulCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&sigma[i][j], sizeof(double), 1, fp);
fwrite(&cut_lj[i][j], sizeof(double), 1, fp);
fwrite(&cut_coul[i][j], sizeof(double), 1, fp);
}
}
}
@ -361,23 +359,23 @@ void PairLJClass2CoulCut::read_restart(FILE *fp)
read_restart_settings(fp);
allocate();
int i,j;
int i, j;
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,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_coul[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut_coul[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut_coul[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
@ -388,11 +386,11 @@ void PairLJClass2CoulCut::read_restart(FILE *fp)
void PairLJClass2CoulCut::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&tail_flag,sizeof(int),1,fp);
fwrite(&cut_lj_global, sizeof(double), 1, fp);
fwrite(&cut_coul_global, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -402,17 +400,17 @@ void PairLJClass2CoulCut::write_restart_settings(FILE *fp)
void PairLJClass2CoulCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_coul_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut_coul_global, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut_coul_global, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
}
/* ----------------------------------------------------------------------
@ -421,8 +419,7 @@ void PairLJClass2CoulCut::read_restart_settings(FILE *fp)
void PairLJClass2CoulCut::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]);
}
/* ----------------------------------------------------------------------
@ -433,39 +430,38 @@ void PairLJClass2CoulCut::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]);
fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
{
double r2inv,rinv,r3inv,r6inv,forcecoul,forcelj,phicoul,philj;
double r2inv, rinv, r3inv, r6inv, forcecoul, forcelj, phicoul, philj;
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
else forcecoul = 0.0;
forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv);
else
forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
r3inv = r2inv * rinv;
r6inv = r3inv * r3inv;
forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]);
} else
forcelj = 0.0;
fforce = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv;
double eng = 0.0;
if (rsq < cut_coulsq[itype][jtype]) {
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
eng += factor_coul*phicoul;
phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv);
eng += factor_coul * phicoul;
}
if (rsq < cut_ljsq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
eng += factor_lj*philj;
philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype];
eng += factor_lj * philj;
}
return eng;
@ -476,9 +472,8 @@ double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype,
void *PairLJClass2CoulCut::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma;
if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
if (strcmp(str, "sigma") == 0) return (void *) sigma;
return nullptr;
}

647
src/CLASS2/pair_lj_class2_coul_long.cpp
View File

File diff suppressed because it is too large Load Diff

7
src/DPD-REACT/pair_multi_lucy.cpp
View File

@ -37,6 +37,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -104,7 +105,6 @@ void PairMultiLucy::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double pi = MathConst::MY_PI;
double A_i;
double A_j;
double fraction_i,fraction_j;
@ -198,7 +198,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
evdwl,0.0,0.0,0.0,0.0,0.0);
@ -733,7 +733,6 @@ void PairMultiLucy::computeLocalDensity()
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double pi = MathConst::MY_PI;
double *rho = atom->rho;
// zero out density
@ -766,7 +765,7 @@ void PairMultiLucy::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
factor= (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;

10
src/DPD-REACT/pair_multi_lucy_rx.cpp
View File

@ -39,6 +39,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -127,7 +128,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
double pi = MathConst::MY_PI;
double A_i, A_j;
double fraction_i,fraction_j;
int jtable;
@ -276,7 +276,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwlOld = mixWtSite1old_i*evdwl;
evdwl = mixWtSite1_i*evdwl;
@ -866,15 +866,13 @@ void PairMultiLucyRX::computeLocalDensity()
const int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
const double pi = MathConst::MY_PI;
const bool newton_pair = force->newton_pair;
const bool one_type = (atom->ntypes == 1);
// Special cut-off values for when there's only one type.
const double cutsq_type11 = cutsq[1][1];
const double rcut_type11 = sqrt(cutsq_type11);
const double factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11);
const double factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11);
double *rho = atom->rho;
@ -925,7 +923,7 @@ void PairMultiLucyRX::computeLocalDensity()
const double rcut = sqrt(cutsq[itype][jtype]);
const double tmpFactor = 1.0-sqrt(rsq)/rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho_i += factor;
if (newton_pair || j < nlocal)
rho[j] += factor;

35
src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp
View File

@ -23,20 +23,24 @@
------------------------------------------------------------------------------------------- */
#include "pair_multi_lucy_rx_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kokkos.h"
#include "math_const.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include <cmath>
#include <cstring>
#include "math_const.h"
#include "atom_kokkos.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
#include "neigh_request.h"
#include "kokkos.h"
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -282,7 +286,6 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEI
double mixWtSite2old_i,mixWtSite2old_j;
double mixWtSite1_i;
double pi = MathConst::MY_PI;
double A_i, A_j;
double fraction_i,fraction_j;
int jtable;
@ -417,7 +420,7 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEI
evdwl = d_table_const.e(tidx,itable) + fraction_i*d_table_const.de(tidx,itable);
} else k_error_flag.template view<DeviceType>()() = 3;
evdwl *=(pi*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0;
evdwl *=(MY_PI*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0;
evdwlOld = mixWtSite1old_i*evdwl;
evdwl = mixWtSite1_i*evdwl;
@ -458,15 +461,13 @@ void PairMultiLucyRXKokkos<DeviceType>::computeLocalDensity()
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
const double pi = MathConst::MY_PI;
const bool newton_pair = force->newton_pair;
const bool one_type = (atom->ntypes == 1);
// Special cut-off values for when there's only one type.
cutsq_type11 = cutsq[1][1];
rcut_type11 = sqrt(cutsq_type11);
factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11);
factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11);
// zero out density
int m = nlocal;
@ -548,8 +549,6 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXComputeLoca
const int itype = type[i];
const int jnum = d_numneigh[i];
const double pi = MathConst::MY_PI;
for (int jj = 0; jj < jnum; jj++) {
const int j = (d_neighbors(i,jj) & NEIGHMASK);
const int jtype = type[j];
@ -574,7 +573,7 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXComputeLoca
const double rcut = sqrt(d_cutsq(itype,jtype));
const double tmpFactor = 1.0-sqrt(rsq)/rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho_i_contrib += factor;
if ((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR || j < nlocal))
a_rho[j] += factor;

61
src/LATTE/fix_latte.cpp
View File

@ -17,20 +17,18 @@
------------------------------------------------------------------------- */
#include "fix_latte.h"
#include <cstdio>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "update.h"
#include "neighbor.h"
#include "domain.h"
#include "force.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "modify.h"
#include "compute.h"
#include "memory.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -82,13 +80,12 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[3],"NULL") != 0) {
coulomb = 1;
error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation "
"of a Coulomb potential");
error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation of a Coulomb potential");
id_pe = utils::strdup(arg[3]);
int ipe = modify->find_compute(id_pe);
if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID");
if (modify->compute[ipe]->peatomflag == 0)
c_pe = modify->get_compute_by_id(id_pe);
if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe);
if (c_pe->peatomflag == 0)
error->all(FLERR,"Fix latte compute ID does not compute pe/atom");
}
@ -105,7 +102,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) :
FixLatte::~FixLatte()
{
delete [] id_pe;
delete[] id_pe;
memory->destroy(qpotential);
memory->destroy(flatte);
}
@ -134,9 +131,8 @@ void FixLatte::init()
if (atom->q_flag == 0 || force->pair == nullptr || force->kspace == nullptr)
error->all(FLERR,"Fix latte cannot compute Coulomb potential");
int ipe = modify->find_compute(id_pe);
if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID");
c_pe = modify->compute[ipe];
c_pe = modify->get_compute_by_id(id_pe);
if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe);
}
// must be fully periodic or fully non-periodic
@ -153,33 +149,6 @@ void FixLatte::init()
memory->create(qpotential,atom->nlocal,"latte:qpotential");
memory->create(flatte,atom->nlocal,3,"latte:flatte");
}
/*
// warn if any integrate fix comes after this one
// is it actually necessary for q(n) update to come after x,v update ??
int after = 0;
int flag = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(id,modify->fix[i]->id) == 0) after = 1;
else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1;
}
if (flag && comm->me == 0)
error->warning(FLERR,"Fix latte should come after all other "
"integration fixes");
*/
/*
// need a full neighbor list
// could we use a half list?
// perpetual list, built whenever re-neighboring occurs
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
*/
}
/* ---------------------------------------------------------------------- */

2
src/MANYBODY/pair_sw.cpp
View File

@ -292,7 +292,7 @@ void PairSW::read_file(char *file)
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "sw", unit_convert_flag);
char * line;
char *line;
// transparently convert units for supported conversions

7
src/MANYBODY/pair_tersoff_table.cpp
View File

@ -27,6 +27,7 @@
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
@ -37,6 +38,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
#define MAXLINE 1024
#define DELTA 4
@ -543,7 +545,6 @@ void PairTersoffTable::allocateGrids()
double deltaArgumentCutoffFunction, deltaArgumentExponential, deltaArgumentBetaZetaPower;
double deltaArgumentGtetaFunction;
double r, minMu, maxLambda, maxCutoff;
double const PI=acos(-1.0);
deallocateGrids();
@ -652,8 +653,8 @@ void PairTersoffTable::allocateGrids()
}
for (l = numGridPointsOneCutoffFunction; l < numGridPointsCutoffFunction; l++) {
cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (PI * (r - cutoffR)/(cutoffS-cutoffR)) ;
cutoffFunctionDerived[i][j][l] = -0.5 * PI * sin (PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR) ;
cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) ;
cutoffFunctionDerived[i][j][l] = -0.5 * MY_PI * sin (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR);
r += deltaArgumentCutoffFunction;
}
}

1
src/MEAM/meam_setup_param.cpp
View File

@ -233,7 +233,6 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3
// 21 = theta
// see alloyparams(void) in meam_setup_done.cpp
case 21:
// const double PI = 3.141592653589793238463;
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;

10
src/OPENMP/fix_rigid_nh_omp.cpp
View File

@ -246,12 +246,11 @@ void FixRigidNHOMP::compute_forces_and_torques()
if (rstyle == SINGLE) {
// we have just one rigid body. use OpenMP reduction to get sum[]
double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0;
int i;
#if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE private(i) reduction(+:s0,s1,s2,s3,s4,s5)
#pragma omp parallel for LMP_DEFAULT_NONE reduction(+:s0,s1,s2,s3,s4,s5)
#endif
for (i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
const int ibody = body[i];
if (ibody < 0) continue;
@ -285,12 +284,11 @@ void FixRigidNHOMP::compute_forces_and_torques()
for (int ib = 0; ib < nbody; ++ib) {
double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0;
int i;
#if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE private(i) LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5)
#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5)
#endif
for (i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
const int ibody = body[i];
if (ibody != ib) continue;

9
src/OPENMP/pair_comb_omp.cpp
View File

@ -384,7 +384,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr)
double PairCombOMP::yasu_char(double *qf_fix, int &igroup)
{
int ii;
double potal,fac11,fac11e;
const double * const * const x = atom->x;
@ -401,7 +400,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup)
const int groupbit = group->bitmask[igroup];
qf = qf_fix;
for (ii = 0; ii < inum; ii++) {
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit)
qf[i] = 0.0;
@ -417,9 +416,9 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup)
// loop over full neighbor list of my atoms
#if defined(_OPENMP)
#pragma omp parallel for private(ii) LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e)
#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e)
#endif
for (ii = 0; ii < inum; ii ++) {
for (int ii = 0; ii < inum; ii ++) {
double fqi,fqij,fqji,fqjj,delr1[3];
double sr1,sr2,sr3;
int mr1,mr2,mr3;
@ -533,7 +532,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup)
// sum charge force on each node and return it
double eneg = 0.0;
for (ii = 0; ii < inum; ii++) {
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit)
eneg += qf[i];

19
src/OPENMP/pair_reaxff_omp.cpp
View File

@ -402,8 +402,7 @@ int PairReaxFFOMP::estimate_reax_lists()
int PairReaxFFOMP::write_reax_lists()
{
int itr_i, itr_j, i, j, num_mynbrs;
int *jlist;
int num_mynbrs;
double d_sqr, dist, cutoff_sqr;
rvec dvec;
@ -425,19 +424,19 @@ int PairReaxFFOMP::write_reax_lists()
num_nbrs = 0;
for (itr_i = 0; itr_i < numall; ++itr_i) {
i = ilist[itr_i];
for (int itr_i = 0; itr_i < numall; ++itr_i) {
int i = ilist[itr_i];
num_nbrs_offset[i] = num_nbrs;
num_nbrs += numneigh[i];
}
#if defined(_OPENMP)
#pragma omp parallel for schedule(dynamic,50) default(shared) \
private(itr_i, itr_j, i, j, jlist, cutoff_sqr, num_mynbrs, d_sqr, dvec, dist)
private(cutoff_sqr, num_mynbrs, d_sqr, dvec, dist)
#endif
for (itr_i = 0; itr_i < numall; ++itr_i) {
i = ilist[itr_i];
jlist = firstneigh[i];
for (int itr_i = 0; itr_i < numall; ++itr_i) {
int i = ilist[itr_i];
auto jlist = firstneigh[i];
Set_Start_Index(i, num_nbrs_offset[i], far_nbrs);
if (i < inum)
@ -447,8 +446,8 @@ int PairReaxFFOMP::write_reax_lists()
num_mynbrs = 0;
for (itr_j = 0; itr_j < numneigh[i]; ++itr_j) {
j = jlist[itr_j];
for (int itr_j = 0; itr_j < numneigh[i]; ++itr_j) {
int j = jlist[itr_j];
j &= NEIGHMASK;
get_distance(x[j], x[i], &d_sqr, &dvec);

8
src/RIGID/fix_rigid_nh_small.cpp
View File

@ -1192,13 +1192,7 @@ void FixRigidNHSmall::compute_dof()
for (int k = 0; k < dimension; k++)
if (fabs(b->inertia[k]) < EPSILON) nf_r--;
}
} else if (dimension == 2) {
nf_r = nlocal_body;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
if (fabs(b->inertia[2]) < EPSILON) nf_r--;
}
}
} else if (dimension == 2) nf_r = nlocal_body;
double nf[2], nfall[2];
nf[0] = nf_t;

76
src/RIGID/fix_shake.cpp
View File

@ -51,12 +51,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
type_flag(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr),
loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr),
vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr),
shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr),
list(nullptr), b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr),
b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr),
a_count(nullptr), a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr),
a_ave_all(nullptr), a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr),
onemols(nullptr)
shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr),
b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), b_max(nullptr),
b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
@ -172,8 +171,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix shake does not exist");
if (atom->molecules[imol]->nset > 1 && comm->me == 0)
error->warning(FLERR,"Molecule template for "
"fix shake has multiple molecules");
error->warning(FLERR,"Molecule template for fix shake has multiple molecules");
onemols = &atom->molecules[imol];
nmol = onemols[0]->nset;
iarg += 2;
@ -199,6 +197,8 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
int nb = atom->nbondtypes + 1;
b_count = new int[nb];
b_count_all = new int[nb];
b_atom = new int[nb];
b_atom_all = new int[nb];
b_ave = new double[nb];
b_ave_all = new double[nb];
b_max = new double[nb];
@ -228,8 +228,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
find_clusters();
if (comm->me == 0)
utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n",
platform::walltime()-time1);
utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n",platform::walltime()-time1);
// initialize list of SHAKE clusters to constrain
@ -281,32 +280,34 @@ FixShake::~FixShake()
memory->destroy(vtmp);
delete [] bond_flag;
delete [] angle_flag;
delete [] type_flag;
delete [] mass_list;
delete[] bond_flag;
delete[] angle_flag;
delete[] type_flag;
delete[] mass_list;
delete [] bond_distance;
delete [] angle_distance;
delete[] bond_distance;
delete[] angle_distance;
if (output_every) {
delete [] b_count;
delete [] b_count_all;
delete [] b_ave;
delete [] b_ave_all;
delete [] b_max;
delete [] b_max_all;
delete [] b_min;
delete [] b_min_all;
delete[] b_count;
delete[] b_count_all;
delete[] b_atom;
delete[] b_atom_all;
delete[] b_ave;
delete[] b_ave_all;
delete[] b_max;
delete[] b_max_all;
delete[] b_min;
delete[] b_min_all;
delete [] a_count;
delete [] a_count_all;
delete [] a_ave;
delete [] a_ave_all;
delete [] a_max;
delete [] a_max_all;
delete [] a_min;
delete [] a_min_all;
delete[] a_count;
delete[] a_count_all;
delete[] a_ave;
delete[] a_ave_all;
delete[] a_max;
delete[] a_max_all;
delete[] a_min;
delete[] a_min_all;
}
memory->destroy(list);
@ -2473,6 +2474,7 @@ void FixShake::stats()
b_count[i] = 0;
b_ave[i] = b_max[i] = 0.0;
b_min[i] = BIG;
b_atom[i] = -1;
}
for (i = 0; i < na; i++) {
a_count[i] = 0;
@ -2504,6 +2506,7 @@ void FixShake::stats()
m = shake_type[i][j-1];
b_count[m]++;
b_atom[m] = n;
b_ave[m] += r;
b_max[m] = MAX(b_max[m],r);
b_min[m] = MIN(b_min[m],r);
@ -2547,6 +2550,7 @@ void FixShake::stats()
// sum across all procs
MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world);
MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world);
MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world);
@ -2563,16 +2567,16 @@ void FixShake::stats()
auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n",
update->ntimestep);
for (i = 1; i < nb; i++) {
const auto bcnt = b_count_all[i]/2;
const auto bcnt = b_count_all[i];
if (bcnt)
mesg += fmt::format("Bond: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width,
b_ave_all[i]/bcnt/2.0,b_max_all[i]-b_min_all[i],bcnt);
b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]);
}
for (i = 1; i < na; i++) {
const auto acnt = a_count_all[i]/3;
const auto acnt = a_count_all[i];
if (acnt)
mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width,
a_ave_all[i]/acnt/3.0,a_max_all[i]-a_min_all[i],acnt);
a_ave_all[i]/acnt,a_max_all[i]-a_min_all[i],acnt/3);
}
utils::logmesg(lmp,mesg);
}

8
src/RIGID/fix_shake.h
View File

@ -108,10 +108,10 @@ class FixShake : public Fix {
int nlist, maxlist; // size and max-size of list
// stat quantities
int *b_count, *b_count_all; // counts for each bond type
double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type
double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays
int *a_count, *a_count_all; // ditto for angle types
int *b_count, *b_count_all, *b_atom, *b_atom_all; // counts for each bond type, atoms in bond cluster
double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type
double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays
int *a_count, *a_count_all; // ditto for angle types
double *a_ave, *a_max, *a_min;
double *a_ave_all, *a_max_all, *a_min_all;

7
src/thermo.cpp
View File

@ -905,7 +905,8 @@ void Thermo::parse_fields(char *str)
if (argi.get_type() == ArgInfo::COMPUTE) {
auto icompute = modify->get_compute_by_id(argi.get_name());
if (!icompute) error->all(FLERR,"Could not find thermo custom compute ID");
if (!icompute)
error->all(FLERR,"Could not find thermo custom compute ID: {}", argi.get_name());
if (argindex1[nfield] == 0 && icompute->scalar_flag == 0)
error->all(FLERR,"Thermo compute does not compute scalar");
if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
@ -935,7 +936,7 @@ void Thermo::parse_fields(char *str)
} else if (argi.get_type() == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(argi.get_name());
if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID");
if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID: {}", argi.get_name());
if (argindex1[nfield] == 0 && ifix->scalar_flag == 0)
error->all(FLERR,"Thermo fix does not compute scalar");
if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
@ -961,7 +962,7 @@ void Thermo::parse_fields(char *str)
} else if (argi.get_type() == ArgInfo::VARIABLE) {
int n = input->variable->find(argi.get_name());
if (n < 0)
error->all(FLERR,"Could not find thermo custom variable name");
error->all(FLERR,"Could not find thermo custom variable name: {}", argi.get_name());
if (argindex1[nfield] == 0 && input->variable->equalstyle(n) == 0)
error->all(FLERR,"Thermo custom variable is not equal-style variable");
if (argindex1[nfield] && input->variable->vectorstyle(n) == 0)

5
unittest/utils/test_platform.cpp
View File

@ -21,7 +21,7 @@ TEST(Platform, clock)
// spend some time computing pi
constexpr double known_pi = 3.141592653589793238462643;
constexpr int n = 10000000;
constexpr int n = 100000000;
constexpr double h = 1.0 / (double)n;
double my_pi = 0.0, x;
for (int i = 0; i < n; ++i) {
@ -34,7 +34,8 @@ TEST(Platform, clock)
ASSERT_NEAR(my_pi, known_pi, 1e-12);
ASSERT_GT(wt_used, 1e-4);
ASSERT_GT(ct_used, 1e-4);
// windows sometimes incorrectly reports used CPU time as 0.0
if (ct_used != 0.0) ASSERT_GT(ct_used, 1e-4);
}
TEST(Platform, putenv_unsetenv)