diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index c48d329998..e76d040c09 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -77,7 +77,7 @@ check_for_autogen_files(${LAMMPS_SOURCE_DIR}) include(CheckIncludeFileCXX) # set required compiler flags and compiler/CPU arch specific optimizations -if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") +if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -restrict") if(CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.3 OR CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.4) set(CMAKE_TUNE_DEFAULT "-xCOMMON-AVX512") diff --git a/cmake/Modules/OpenCLLoader.cmake b/cmake/Modules/OpenCLLoader.cmake index ea83f8e447..239fba5ba9 100644 --- a/cmake/Modules/OpenCLLoader.cmake +++ b/cmake/Modules/OpenCLLoader.cmake @@ -1,6 +1,6 @@ message(STATUS "Downloading and building OpenCL loader library") -set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2021.05.02.tar.gz" CACHE STRING "URL for OpenCL loader tarball") -set(OPENCL_LOADER_MD5 "29180b05056578afda92f0d394c3a373" CACHE STRING "MD5 checksum of OpenCL loader tarball") +set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2021.06.30.tar.gz" CACHE STRING "URL for OpenCL loader tarball") +set(OPENCL_LOADER_MD5 "f9e55dd550cfbf77f46507adf7cb8fd2" CACHE STRING "MD5 checksum of OpenCL loader tarball") mark_as_advanced(OPENCL_LOADER_URL) mark_as_advanced(OPENCL_LOADER_MD5) diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index aa34516d52..a110d93a8d 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -188,6 +188,8 @@ formulas for the meaning of these parameters: +------------------------------------------------------------------------------+--------------------------------------------------+-------------+ | :doc:`reax/c ` | chi, eta, gamma | type global | +------------------------------------------------------------------------------+--------------------------------------------------+-------------+ +| :doc:`snap ` | scale | type pairs | ++------------------------------------------------------------------------------+--------------------------------------------------+-------------+ | :doc:`spin/dmi ` | coulombic_cutoff | type global | +------------------------------------------------------------------------------+--------------------------------------------------+-------------+ | :doc:`spin/exchange ` | coulombic_cutoff | type global | diff --git a/doc/src/fix_adapt_fep.rst b/doc/src/fix_adapt_fep.rst index 10a1faf202..0f669eade3 100644 --- a/doc/src/fix_adapt_fep.rst +++ b/doc/src/fix_adapt_fep.rst @@ -173,10 +173,12 @@ styles and their energy formulas for the meaning of these parameters: +------------------------------------------------------------------------------+-------------------------+------------+ | :doc:`nm/cut/coul/cut, nm/cut/coul/long ` | e0,r0,nn,mm | type pairs | +------------------------------------------------------------------------------+-------------------------+------------+ -| :doc:`ufm ` | epsilon,sigma,scale | type pairs | +| :doc:`snap ` | scale | type pairs | +------------------------------------------------------------------------------+-------------------------+------------+ | :doc:`soft ` | a | type pairs | +------------------------------------------------------------------------------+-------------------------+------------+ +| :doc:`ufm ` | epsilon,sigma,scale | type pairs | ++------------------------------------------------------------------------------+-------------------------+------------+ .. note:: diff --git a/doc/src/pair_rann.rst b/doc/src/pair_rann.rst index b02a81e4b8..a9ad5d9236 100644 --- a/doc/src/pair_rann.rst +++ b/doc/src/pair_rann.rst @@ -382,6 +382,8 @@ Pair style *rann* is part of the ML-RANN package. It is only enabled if LAMMPS package. Additionally, if any spin fingerprint styles are used LAMMPS must be built with the SPIN package as well. +Pair style *rann* does not support computing per-atom stress or using :doc:`pair_modify nofdotr `. + Defaults """""""""""" diff --git a/doc/src/pair_tersoff.rst b/doc/src/pair_tersoff.rst index 0f0750e9df..4fa77fc032 100644 --- a/doc/src/pair_tersoff.rst +++ b/doc/src/pair_tersoff.rst @@ -16,9 +16,6 @@ pair_style tersoff/table command Accelerator Variants: *tersoff/table/omp* -pair_style tersoff/shift command -================================ - Syntax """""" @@ -309,6 +306,9 @@ The *shift* keyword is not supported by the *tersoff/gpu*, *tersoff/intel*, *tersoff/kk*, *tersoff/table* or *tersoff/table/omp* variants. +The *tersoff/intel* pair style is only available when compiling LAMMPS +with the Intel compilers. + The Tersoff potential files provided with LAMMPS (see the potentials directory) are parameterized for :doc:`"metal" units `. In addition the pair style supports converting potential parameters on-the-fly between diff --git a/doc/utils/requirements.txt b/doc/utils/requirements.txt index 9b8e106875..7e4563a1ec 100644 --- a/doc/utils/requirements.txt +++ b/doc/utils/requirements.txt @@ -1,4 +1,4 @@ -Sphinx +Sphinx==4.0.3 sphinxcontrib-spelling git+git://github.com/akohlmey/sphinx-fortran@parallel-read sphinx_tabs diff --git a/examples/plugins/CMakeLists.txt b/examples/plugins/CMakeLists.txt index 1169ab490a..59c2802b45 100644 --- a/examples/plugins/CMakeLists.txt +++ b/examples/plugins/CMakeLists.txt @@ -31,6 +31,11 @@ endif() set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) +# Need -restrict with Intel compilers +if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -restrict") +endif() + # bail out on windows if(CMAKE_SYSTEM_NAME STREQUAL Windows) message(FATAL_ERROR "LAMMPS plugins are currently not supported on Windows") diff --git a/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapcoeff b/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapcoeff new file mode 120000 index 0000000000..9401393007 --- /dev/null +++ b/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapcoeff @@ -0,0 +1 @@ +../../potentials/Ni_Zuo_JPCA2020.quadratic.snapcoeff \ No newline at end of file diff --git a/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapparam b/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapparam new file mode 120000 index 0000000000..a47ada019b --- /dev/null +++ b/examples/snap/Ni_Zuo_JPCA2020.quadratic.snapparam @@ -0,0 +1 @@ +../../potentials/Ni_Zuo_JPCA2020.quadratic.snapparam \ No newline at end of file diff --git a/examples/snap/Ni_Zuo_JPCA2020.snapcoeff b/examples/snap/Ni_Zuo_JPCA2020.snapcoeff new file mode 120000 index 0000000000..7211856b43 --- /dev/null +++ b/examples/snap/Ni_Zuo_JPCA2020.snapcoeff @@ -0,0 +1 @@ +../../potentials/Ni_Zuo_JPCA2020.snapcoeff \ No newline at end of file diff --git a/examples/snap/Ni_Zuo_JPCA2020.snapparam b/examples/snap/Ni_Zuo_JPCA2020.snapparam new file mode 120000 index 0000000000..b7fc095fab --- /dev/null +++ b/examples/snap/Ni_Zuo_JPCA2020.snapparam @@ -0,0 +1 @@ +../../potentials/Ni_Zuo_JPCA2020.snapparam \ No newline at end of file diff --git a/examples/snap/in.snap.scale.Ni_Zuo_JCPA2020 b/examples/snap/in.snap.scale.Ni_Zuo_JCPA2020 new file mode 100644 index 0000000000..ee195c11f3 --- /dev/null +++ b/examples/snap/in.snap.scale.Ni_Zuo_JCPA2020 @@ -0,0 +1,53 @@ +# Toy demonstration of SNAP "scale" parameter, using fix/adapt and hybrid/overlay +# Mixing linear and quadratic SNAP Ni potentials by Zuo et al. JCPA 2020 + +# mixing parameter + +variable lambda equal 0.2 + + +# Initialize simulation + +variable nsteps index 100 +variable nrep equal 3 +variable a equal 3.52 +units metal + +# generate the box and atom positions using a FCC lattice +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice fcc $a +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 34. + +# choose bundled SNAP Ni potential from Zuo et al. JCPA 2020 +pair_style hybrid/overlay snap snap +pair_coeff * * snap 1 Ni_Zuo_JPCA2020.snapcoeff Ni_Zuo_JPCA2020.snapparam Ni +pair_coeff * * snap 2 Ni_Zuo_JPCA2020.quadratic.snapcoeff Ni_Zuo_JPCA2020.quadratic.snapparam Ni + +# scale according to mixing parameter +variable l1 equal ${lambda} +variable l2 equal 1.0-${lambda} +fix scale1 all adapt 1 pair snap:1 scale * * v_l1 +fix scale2 all adapt 1 pair snap:2 scale * * v_l2 + +# Setup output +thermo 1 +thermo_modify norm yes + +# Set up NVE run +timestep 0.5e-3 +neighbor 1.0 bin +neigh_modify every 1 delay 0 check yes + +# Run MD +velocity all create 300.0 4928459 loop geom +fix 1 all nve +run ${nsteps} diff --git a/src/INTEL/pair_airebo_intel.cpp b/src/INTEL/pair_airebo_intel.cpp index 8b2eadbe72..94874153aa 100644 --- a/src/INTEL/pair_airebo_intel.cpp +++ b/src/INTEL/pair_airebo_intel.cpp @@ -180,12 +180,10 @@ PairAIREBOIntel::~PairAIREBOIntel() void PairAIREBOIntel::init_style() { PairAIREBO::init_style(); - neighbor->requests[neighbor->nrequest-1]->intel = 1; + neighbor->find_request(this)->intel = 1; - const int nrequest = neighbor->nrequest; - for (int i = 0; i < nrequest; ++i) - if (neighbor->requests[i]->skip) - error->all(FLERR, "Cannot yet use airebo/intel with hybrid."); + if (utils::strmatch(force->pair_style,"^hybrid")) + error->all(FLERR, "Cannot yet use airebo/intel with hybrid."); int ifix = modify->find_fix("package_intel"); if (ifix < 0) @@ -294,6 +292,8 @@ void PairAIREBOIntel::compute( ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); pvector[0] = pvector[1] = pvector[2] = 0.0; diff --git a/src/INTEL/pair_buck_coul_cut_intel.cpp b/src/INTEL/pair_buck_coul_cut_intel.cpp index 1859edb732..99905bfaa0 100644 --- a/src/INTEL/pair_buck_coul_cut_intel.cpp +++ b/src/INTEL/pair_buck_coul_cut_intel.cpp @@ -16,25 +16,25 @@ Contributing author: Rodrigo Canales (RWTH Aachen University) ------------------------------------------------------------------------- */ -#include -#include - -#include #include "pair_buck_coul_cut_intel.h" + #include "atom.h" #include "comm.h" +#include "error.h" +#include "force.h" #include "force.h" #include "group.h" #include "kspace.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "math_const.h" #include "memory.h" -#include "error.h" -#include "suffix.h" -#include "force.h" #include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "suffix.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -77,6 +77,8 @@ void PairBuckCoulCutIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -403,11 +405,12 @@ void PairBuckCoulCutIntel::eval(const int offload, const int vflag, void PairBuckCoulCutIntel::init_style() { PairBuckCoulCut::init_style(); + auto request = neighbor->find_request(this); if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_buck_coul_long_intel.cpp b/src/INTEL/pair_buck_coul_long_intel.cpp index dc5eed7521..1566ec23b6 100644 --- a/src/INTEL/pair_buck_coul_long_intel.cpp +++ b/src/INTEL/pair_buck_coul_long_intel.cpp @@ -16,25 +16,25 @@ Contributing author: Rodrigo Canales (RWTH Aachen University) ------------------------------------------------------------------------- */ -#include -#include - -#include #include "pair_buck_coul_long_intel.h" + #include "atom.h" #include "comm.h" +#include "error.h" +#include "force.h" #include "force.h" #include "group.h" #include "kspace.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "math_const.h" #include "memory.h" -#include "error.h" -#include "suffix.h" -#include "force.h" #include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "suffix.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -77,6 +77,8 @@ void PairBuckCoulLongIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -479,11 +481,13 @@ void PairBuckCoulLongIntel::eval(const int offload, const int vflag, void PairBuckCoulLongIntel::init_style() { PairBuckCoulLong::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_buck_intel.cpp b/src/INTEL/pair_buck_intel.cpp index 23f7852486..26ef13be9a 100644 --- a/src/INTEL/pair_buck_intel.cpp +++ b/src/INTEL/pair_buck_intel.cpp @@ -70,6 +70,8 @@ void PairBuckIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -363,11 +365,14 @@ void PairBuckIntel::eval(const int offload, const int vflag, void PairBuckIntel::init_style() { PairBuck::init_style(); + + // augment neighbor list request + auto request = neighbor->find_request(this); if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_dpd_intel.cpp b/src/INTEL/pair_dpd_intel.cpp index 0e872efdf6..089396afa3 100644 --- a/src/INTEL/pair_dpd_intel.cpp +++ b/src/INTEL/pair_dpd_intel.cpp @@ -14,17 +14,20 @@ Shun Xu (Computer Network Information Center, CAS) ------------------------------------------------------------------------- */ -#include #include "pair_dpd_intel.h" + #include "atom.h" #include "comm.h" #include "force.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include "suffix.h" + +#include + using namespace LAMMPS_NS; #define LMP_MKL_RNG VSL_BRNG_MT19937 @@ -86,6 +89,8 @@ void PairDPDIntel::compute(int eflag, int vflag, ev_init(eflag, vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -469,11 +474,13 @@ void PairDPDIntel::settings(int narg, char **arg) { void PairDPDIntel::init_style() { PairDPD::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_eam_intel.cpp b/src/INTEL/pair_eam_intel.cpp index 04724f599c..a86040b6b3 100644 --- a/src/INTEL/pair_eam_intel.cpp +++ b/src/INTEL/pair_eam_intel.cpp @@ -82,6 +82,8 @@ void PairEAMIntel::compute(int eflag, int vflag, ev_init(eflag, vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -638,11 +640,13 @@ void PairEAMIntel::eval(const int offload, const int vflag, void PairEAMIntel::init_style() { PairEAM::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_gayberne_intel.cpp b/src/INTEL/pair_gayberne_intel.cpp index 8b81b1ea81..d7becc7585 100644 --- a/src/INTEL/pair_gayberne_intel.cpp +++ b/src/INTEL/pair_gayberne_intel.cpp @@ -25,16 +25,16 @@ #endif #include "atom.h" -#include "comm.h" #include "atom_vec_ellipsoid.h" +#include "comm.h" #include "force.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" - +#include "neighbor.h" #include "suffix.h" + using namespace LAMMPS_NS; #define FC_PACKED1_T typename ForceConst::fc_packed1 @@ -76,6 +76,8 @@ void PairGayBerneIntel::compute(int eflag, int vflag, ev_init(eflag, vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nall = atom->nlocal + atom->nghost; @@ -874,11 +876,13 @@ void PairGayBerneIntel::eval(const int offload, const int vflag, void PairGayBerneIntel::init_style() { PairGayBerne::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_lj_charmm_coul_charmm_intel.cpp b/src/INTEL/pair_lj_charmm_coul_charmm_intel.cpp index f67f5156b8..ad8ef4d84f 100644 --- a/src/INTEL/pair_lj_charmm_coul_charmm_intel.cpp +++ b/src/INTEL/pair_lj_charmm_coul_charmm_intel.cpp @@ -73,6 +73,8 @@ void PairLJCharmmCoulCharmmIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -448,11 +450,13 @@ void PairLJCharmmCoulCharmmIntel::eval(const int offload, const int vflag, void PairLJCharmmCoulCharmmIntel::init_style() { PairLJCharmmCoulCharmm::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_lj_charmm_coul_long_intel.cpp b/src/INTEL/pair_lj_charmm_coul_long_intel.cpp index 40a18be7d5..a910c74acb 100644 --- a/src/INTEL/pair_lj_charmm_coul_long_intel.cpp +++ b/src/INTEL/pair_lj_charmm_coul_long_intel.cpp @@ -13,20 +13,23 @@ Contributing author: W. Michael Brown (Intel) ------------------------------------------------------------------------- */ -#include #include "pair_lj_charmm_coul_long_intel.h" + #include "atom.h" #include "comm.h" #include "force.h" #include "group.h" #include "kspace.h" #include "memory.h" +#include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" -#include "memory.h" +#include "neighbor.h" #include "suffix.h" + +#include + using namespace LAMMPS_NS; #define LJ_T typename IntelBuffers::vec2_t @@ -74,6 +77,8 @@ void PairLJCharmmCoulLongIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -510,11 +515,13 @@ void PairLJCharmmCoulLongIntel::eval(const int offload, const int vflag, void PairLJCharmmCoulLongIntel::init_style() { PairLJCharmmCoulLong::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_lj_cut_coul_long_intel.cpp b/src/INTEL/pair_lj_cut_coul_long_intel.cpp index 47d14186f6..51e208314b 100644 --- a/src/INTEL/pair_lj_cut_coul_long_intel.cpp +++ b/src/INTEL/pair_lj_cut_coul_long_intel.cpp @@ -76,6 +76,8 @@ void PairLJCutCoulLongIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -476,11 +478,13 @@ void PairLJCutCoulLongIntel::eval(const int offload, const int vflag, void PairLJCutCoulLongIntel::init_style() { PairLJCutCoulLong::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_lj_cut_intel.cpp b/src/INTEL/pair_lj_cut_intel.cpp index 71e12d7b24..eb60e0442f 100644 --- a/src/INTEL/pair_lj_cut_intel.cpp +++ b/src/INTEL/pair_lj_cut_intel.cpp @@ -13,18 +13,20 @@ Contributing author: W. Michael Brown (Intel) ------------------------------------------------------------------------- */ -#include #include "pair_lj_cut_intel.h" + #include "atom.h" #include "comm.h" #include "force.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" - +#include "neighbor.h" #include "suffix.h" + +#include + using namespace LAMMPS_NS; #define FC_PACKED1_T typename ForceConst::fc_packed1 @@ -66,6 +68,8 @@ void PairLJCutIntel::compute(int eflag, int vflag, ev_init(eflag, vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -380,11 +384,13 @@ void PairLJCutIntel::eval(const int offload, const int vflag, void PairLJCutIntel::init_style() { PairLJCut::init_style(); + auto request = neighbor->find_request(this); + if (force->newton_pair == 0) { - neighbor->requests[neighbor->nrequest-1]->half = 0; - neighbor->requests[neighbor->nrequest-1]->full = 1; + request->half = 0; + request->full = 1; } - neighbor->requests[neighbor->nrequest-1]->intel = 1; + request->intel = 1; int ifix = modify->find_fix("package_intel"); if (ifix < 0) diff --git a/src/INTEL/pair_sw_intel.cpp b/src/INTEL/pair_sw_intel.cpp index 421e8b2fc5..17dffa2843 100644 --- a/src/INTEL/pair_sw_intel.cpp +++ b/src/INTEL/pair_sw_intel.cpp @@ -97,6 +97,8 @@ void PairSWIntel::compute(int eflag, int vflag, ev_init(eflag, vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -1102,7 +1104,8 @@ void PairSWIntel::allocate() void PairSWIntel::init_style() { PairSW::init_style(); - neighbor->requests[neighbor->nrequest-1]->intel = 1; + neighbor->find_request(this)->intel = 1; + map[0] = map[1]; int ifix = modify->find_fix("package_intel"); diff --git a/src/INTEL/pair_tersoff_intel.cpp b/src/INTEL/pair_tersoff_intel.cpp index 707e8404ff..248c3420e6 100644 --- a/src/INTEL/pair_tersoff_intel.cpp +++ b/src/INTEL/pair_tersoff_intel.cpp @@ -16,11 +16,8 @@ Contributing author: Markus Höhnerbach (RWTH) ------------------------------------------------------------------------- */ -#include -#include - -#include #include "pair_tersoff_intel.h" + #include "atom.h" #include "neighbor.h" #include "neigh_list.h" @@ -30,30 +27,13 @@ #include "memory.h" #include "error.h" -// Currently Intel compiler is required for this pair style. -// For convenience, base class routines are called if not using Intel compiler. -#ifndef __INTEL_COMPILER +#include +#include + using namespace LAMMPS_NS; -PairTersoffIntel::PairTersoffIntel(LAMMPS *lmp) : PairTersoff(lmp) -{ -} - -void PairTersoffIntel::compute(int eflag, int vflag) -{ - PairTersoff::compute(eflag, vflag); -} - -void PairTersoffIntel::init_style() -{ - if (comm->me == 0) { - error->warning(FLERR, "Tersoff/intel currently requires intel compiler. " - "Using MANYBODY version."); - } - PairTersoff::init_style(); -} - -#else +// Currently the Intel compiler is required for this pair style. +#ifdef __INTEL_COMPILER #ifdef _LMP_INTEL_OFFLOAD #pragma offload_attribute(push,target(mic)) @@ -111,6 +91,8 @@ void PairTersoffIntel::compute(int eflag, int vflag, ev_init(eflag,vflag); if (vflag_atom) error->all(FLERR,"INTEL package does not support per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"INTEL package does not support pair_modify nofdotr"); const int inum = list->inum; const int nthreads = comm->nthreads; @@ -1439,4 +1421,4 @@ void IntelKernelTersoff::attractive_vector( #pragma offload_attribute(pop) #endif -#endif +#endif // __INTEL_COMPILER diff --git a/src/INTEL/pair_tersoff_intel.h b/src/INTEL/pair_tersoff_intel.h index b613a00194..0a40844622 100644 --- a/src/INTEL/pair_tersoff_intel.h +++ b/src/INTEL/pair_tersoff_intel.h @@ -14,13 +14,18 @@ #ifdef PAIR_CLASS // clang-format off +// Currently the Intel compilers are required for this pair style. +#ifdef __INTEL_COMPILER PairStyle(tersoff/intel,PairTersoffIntel); +#endif // clang-format on #else #ifndef LMP_PAIR_TERSOFF_INTEL_H #define LMP_PAIR_TERSOFF_INTEL_H +#ifdef __INTEL_COMPILER + #include "pair.h" #include "fix_intel.h" #include "pair_tersoff.h" @@ -33,7 +38,6 @@ class PairTersoffIntel : public PairTersoff { virtual void compute(int, int); void init_style(); -#ifdef __INTEL_COMPILER protected: typedef struct { float x,y,z; int w; } sng4_t; @@ -88,10 +92,10 @@ class PairTersoffIntel : public PairTersoff { template void pack_force_const(ForceConst &fc, IntelBuffers *buffers); -#endif // __INTEL_COMPILER }; } +#endif // __INTEL_COMPILER #endif #endif diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index b18dbfef90..67ff81d4a3 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -85,10 +85,10 @@ PairAIREBO::~PairAIREBO() { memory->destroy(REBO_numneigh); memory->sfree(REBO_firstneigh); - delete [] ipage; + delete[] ipage; memory->destroy(nC); memory->destroy(nH); - delete [] pvector; + delete[] pvector; if (allocated) { memory->destroy(setflag); @@ -100,7 +100,7 @@ PairAIREBO::~PairAIREBO() memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); - delete [] map; + delete[] map; } } @@ -112,9 +112,9 @@ void PairAIREBO::compute(int eflag, int vflag) pvector[0] = pvector[1] = pvector[2] = 0.0; REBO_neigh(); - FREBO(eflag,vflag); - if (ljflag) FLJ(eflag,vflag); - if (torflag) TORSION(eflag,vflag); + FREBO(eflag); + if (ljflag) FLJ(eflag); + if (torflag) TORSION(eflag); if (vflag_fdotr) virial_fdotr_compute(); } @@ -252,7 +252,7 @@ void PairAIREBO::init_style() if (oneatom != neighbor->oneatom) create = 1; if (create) { - delete [] ipage; + delete[] ipage; pgsize = neighbor->pgsize; oneatom = neighbor->oneatom; @@ -421,7 +421,7 @@ void PairAIREBO::REBO_neigh() REBO forces and energy ------------------------------------------------------------------------- */ -void PairAIREBO::FREBO(int eflag, int /*vflag*/) +void PairAIREBO::FREBO(int eflag) { int i,j,k,m,ii,inum,itype,jtype; tagint itag,jtag; @@ -497,7 +497,7 @@ void PairAIREBO::FREBO(int eflag, int /*vflag*/) del[0] = delx; del[1] = dely; del[2] = delz; - bij = bondorder(i,j,del,rij,VA,f,vflag_atom); + bij = bondorder(i,j,del,rij,VA,f); dVAdi = bij*dVA; fpair = -(dVRdi+dVAdi) / rij; @@ -520,7 +520,7 @@ void PairAIREBO::FREBO(int eflag, int /*vflag*/) find 3- and 4-step paths between atoms I,J via REBO neighbor lists ------------------------------------------------------------------------- */ -void PairAIREBO::FLJ(int eflag, int /*vflag*/) +void PairAIREBO::FLJ(int eflag) { int i,j,k,m,ii,jj,kk,mm,inum,jnum,itype,jtype,ktype,mtype; int atomi,atomj,atomk,atomm; @@ -786,8 +786,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/) delscale[0] = scale * delij[0]; delscale[1] = scale * delij[1]; delscale[2] = scale * delij[2]; - Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA, - delij,rij,f,vflag_atom); + Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,f); } else Stb = 0.0; fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) + @@ -815,7 +814,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/) f[atomj][1] -= delij[1]*fpair; f[atomj][2] -= delij[2]*fpair; - if (vflag_atom) v_tally2(atomi,atomj,fpair,delij); + if (vflag_either) v_tally2(atomi,atomj,fpair,delij); } else if (npath == 3) { fpair1 = dC*dwikS*wkjS / rikS; @@ -837,7 +836,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/) f[atomk][1] -= fi[1] + fj[1]; f[atomk][2] -= fi[2] + fj[2]; - if (vflag_atom) + if (vflag_either) v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS); } else if (npath == 4) { @@ -873,7 +872,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/) f[atomm][1] += fm[1]; f[atomm][2] += fm[2]; - if (vflag_atom) { + if (vflag_either) { delimS[0] = delikS[0] + delkmS[0]; delimS[1] = delikS[1] + delkmS[1]; delimS[2] = delikS[2] + delkmS[2]; @@ -889,7 +888,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/) torsional forces and energy ------------------------------------------------------------------------- */ -void PairAIREBO::TORSION(int eflag, int /*vflag*/) +void PairAIREBO::TORSION(int eflag) { int i,j,k,l,ii,inum; tagint itag,jtag; @@ -1249,9 +1248,7 @@ void PairAIREBO::TORSION(int eflag, int /*vflag*/) Bij function ------------------------------------------------------------------------- */ -double PairAIREBO::bondorder(int i, int j, double rij[3], - double rijmag, double VA, - double **f, int vflag_atom) +double PairAIREBO::bondorder(int i, int j, double rij[3], double rijmag, double VA, double **f) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; @@ -1435,7 +1432,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; - if (vflag_atom) { + if (vflag_either) { rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); @@ -1577,7 +1574,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; - if (vflag_atom) { + if (vflag_either) { rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); } @@ -1613,7 +1610,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -1627,7 +1624,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -1649,7 +1646,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); + if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -1680,7 +1677,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -1694,7 +1691,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -1716,7 +1713,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); + if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -1928,7 +1925,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; f[atom4][2] += f4[2]; - if (vflag_atom) { + if (vflag_either) { r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); @@ -1964,7 +1961,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -1978,7 +1975,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -2000,7 +1997,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); + if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2031,7 +2028,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -2045,7 +2042,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -2067,7 +2064,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); + if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -2113,8 +2110,7 @@ but of the vector r_ij. */ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rijmag_mod, - double VA, double rij[3], double rijmag, - double **f, int vflag_atom) + double VA, double rij[3], double rijmag, double **f) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; @@ -2438,7 +2434,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; - if (vflag_atom) { + if (vflag_either) { rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); @@ -2542,7 +2538,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; - if (vflag_atom) { + if (vflag_either) { rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); } @@ -2577,7 +2573,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -2591,7 +2587,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -2613,7 +2609,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); + if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2644,7 +2640,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -2658,7 +2654,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -2680,7 +2676,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); + if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -2885,7 +2881,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; f[atom4][2] += f4[2]; - if (vflag_atom) { + if (vflag_either) { r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); @@ -2919,7 +2915,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -2933,7 +2929,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -2955,7 +2951,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); + if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2986,7 +2982,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -3000,7 +2996,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -3022,7 +3018,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); + if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln); } } } diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index fa1572b1e4..fc511b7bb0 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -96,12 +96,12 @@ class PairAIREBO : public Pair { double Tf[5][5][10], Tdfdx[5][5][10], Tdfdy[5][5][10], Tdfdz[5][5][10]; void REBO_neigh(); - void FREBO(int, int); - void FLJ(int, int); - void TORSION(int, int); + void FREBO(int); + void FLJ(int); + void TORSION(int); - double bondorder(int, int, double *, double, double, double **, int); - double bondorderLJ(int, int, double *, double, double, double *, double, double **, int); + double bondorder(int, int, double *, double, double, double **); + double bondorderLJ(int, int, double *, double, double, double *, double, double **); double gSpline(double, double, int, double *, double *); double PijSpline(double, double, int, int, double *); diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp index 0dbcb1de15..e2a677c2a0 100644 --- a/src/MANYBODY/pair_bop.cpp +++ b/src/MANYBODY/pair_bop.cpp @@ -251,10 +251,8 @@ void PairBOP::compute(int eflag, int vflag) f[j][1] -= ftmp2; f[j][2] -= ftmp3; dE = pl_ij.rep - 2.0*pl_ij.betaS*sigB_0 - 2.0*pl_ij.betaP*piB_0; - if (evflag) { - ev_tally(i,j,nlocal,newton_pair, dE, 0.0, dpr1, - pl_ij.dis[0],pl_ij.dis[1],pl_ij.dis[2]); - } + if (evflag) ev_tally(i,j,nlocal,newton_pair, dE, 0.0, -dpr1, + pl_ij.dis[0],pl_ij.dis[1],pl_ij.dis[2]); } nlisti = BOP_total2[i]; for (jj = 0; jj < nlisti; jj++) { @@ -274,10 +272,8 @@ void PairBOP::compute(int eflag, int vflag) f[j][1] -= ftmp2; f[j][2] -= ftmp3; dE = -p2_ij.rep; - if (evflag) { - ev_tally(i,j,nlocal,newton_pair, dE, 0.0, dpr2, - p2_ij.dis[0],p2_ij.dis[1],p2_ij.dis[2]); - } + if (evflag) ev_tally(i,j,nlocal,newton_pair, dE, 0.0, -dpr2, + p2_ij.dis[0],p2_ij.dis[1],p2_ij.dis[2]); } } if (vflag_fdotr) virial_fdotr_compute(); @@ -1142,10 +1138,8 @@ double PairBOP::SigmaBo(int itmp, int jtmp) f[bt_i][n] -= ftmp[n]; f[bt_j][n] += ftmp[n]; } - if (evflag) { - ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0, - ftmp[0],ftmp[1],ftmp[2],xtmp[0],xtmp[1],xtmp[2]); - } + if (evflag) ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0, + -ftmp[0],-ftmp[1],-ftmp[2],xtmp[0],xtmp[1],xtmp[2]); } else { for (int n = 0; n < 3; n++) { bt_sg[loop].dSigB[n] = dsigB*part2*bt_sg[loop].dSigB1[n] - @@ -1157,10 +1151,8 @@ double PairBOP::SigmaBo(int itmp, int jtmp) f[bt_i][n] -= ftmp[n]; f[bt_j][n] += ftmp[n]; } - if (evflag) { - ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0, - ftmp[0],ftmp[1],ftmp[2],xtmp[0],xtmp[1],xtmp[2]); - } + if (evflag) ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0, + -ftmp[0],-ftmp[1],-ftmp[2],xtmp[0],xtmp[1],xtmp[2]); } } return(sigB); @@ -1841,10 +1833,8 @@ double PairBOP::PiBo(int itmp, int jtmp) f[bt_i][n] -= ftmp[n]; f[bt_j][n] += ftmp[n]; } - if (evflag) { - ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0,ftmp[0],ftmp[1], - ftmp[2],xtmp[0],xtmp[1],xtmp[2]); - } + if (evflag) ev_tally_xyz(bt_i,bt_j,nlocal,newton_pair,0.0,0.0, + -ftmp[0],-ftmp[1],-ftmp[2],xtmp[0],xtmp[1],xtmp[2]); } return(piB); } diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index fd278d1488..689ba992d8 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -394,7 +394,7 @@ void PairComb::compute(int eflag, int vflag) if (evflag) ev_tally(i,j,nlocal,newton_pair,elp_ij,0.0,0.0,0.0,0.0,0.0); - if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2); } } diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp index c3d9c3043f..eee97849b9 100644 --- a/src/MANYBODY/pair_comb3.cpp +++ b/src/MANYBODY/pair_comb3.cpp @@ -1294,7 +1294,7 @@ void PairComb3::compute(int eflag, int vflag) if (evflag) ev_tally(i,j,nlocal,newton_pair,ep6p_ij,0.0,0.0,0.0,0.0,0.0); - if (vflag_atom) + if (vflag_either) v_tally3(i,j,k,fj,fk,delrj,delrk); } // k-loop @@ -1387,7 +1387,7 @@ void PairComb3::compute(int eflag, int vflag) } } - if (vflag_atom) + if (vflag_either) v_tally3(j,i,l,fi,fl,delrl,delrk); } } diff --git a/src/MANYBODY/pair_gw.cpp b/src/MANYBODY/pair_gw.cpp index 71a0038bb6..98826f98f6 100644 --- a/src/MANYBODY/pair_gw.cpp +++ b/src/MANYBODY/pair_gw.cpp @@ -227,7 +227,7 @@ void PairGW::compute(int eflag, int vflag) f[k][1] += fk[1]; f[k][2] += fk[2]; - if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2); } // kk } // jj } // ii diff --git a/src/MANYBODY/pair_lcbop.cpp b/src/MANYBODY/pair_lcbop.cpp index 1af884d814..f1223e1eba 100644 --- a/src/MANYBODY/pair_lcbop.cpp +++ b/src/MANYBODY/pair_lcbop.cpp @@ -67,7 +67,7 @@ PairLCBOP::~PairLCBOP() { memory->destroy(SR_numneigh); memory->sfree(SR_firstneigh); - delete [] ipage; + delete[] ipage; memory->destroy(N); memory->destroy(M); @@ -167,7 +167,7 @@ void PairLCBOP::init_style() if (oneatom != neighbor->oneatom) create = 1; if (create) { - delete [] ipage; + delete[] ipage; pgsize = neighbor->pgsize; oneatom = neighbor->oneatom; @@ -386,7 +386,7 @@ void PairLCBOP::FSR(int eflag, int /*vflag*/) del[0] = delx; del[1] = dely; del[2] = delz; - Bij = bondorder(i,j,del,rijmag,VA,f,vflag_atom); + Bij = bondorder(i,j,del,rijmag,VA,f); dVAdi = Bij*dVA; // F = (dVRdi+dVAdi)*(-grad rijmag) @@ -402,8 +402,7 @@ void PairLCBOP::FSR(int eflag, int /*vflag*/) double evdwl=0.0; if (eflag) evdwl = VR - Bij*VA; - 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); } } } @@ -493,8 +492,7 @@ void PairLCBOP::FLR(int eflag, int /*vflag*/) double evdwl=0.0; if (eflag) evdwl = V; - 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); } } } @@ -503,7 +501,7 @@ void PairLCBOP::FLR(int eflag, int /*vflag*/) forces for Nij and Mij ------------------------------------------------------------------------- */ -void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) { +void PairLCBOP::FNij( int i, int j, double factor, double **f) { int atomi = i; int atomj = j; int *SR_neighs = SR_firstneigh[i]; @@ -532,7 +530,7 @@ void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) { f[atomk][1] -= rik[1]*fpair; f[atomk][2] -= rik[2]*fpair; - if (vflag_atom) v_tally2(atomi,atomk,fpair,rik); + if (vflag_either) v_tally2(atomi,atomk,fpair,rik); } } } @@ -540,7 +538,7 @@ void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) { /* ---------------------------------------------------------------------- */ -void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) { +void PairLCBOP::FMij( int i, int j, double factor, double **f) { int atomi = i; int atomj = j; int *SR_neighs = SR_firstneigh[i]; @@ -573,12 +571,12 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) { f[atomk][0] -= rik[0]*fpair; f[atomk][1] -= rik[1]*fpair; f[atomk][2] -= rik[2]*fpair; - if (vflag_atom) v_tally2(atomi,atomk,fpair,rik); + if (vflag_either) v_tally2(atomi,atomk,fpair,rik); } if (dF > TOL) { double factor2 = factor*f_c_ik*dF; - FNij( atomk, atomi, factor2, f, vflag_atom ); + FNij(atomk, atomi, factor2, f); } } } @@ -588,17 +586,15 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) { Bij function ------------------------------------------------------------------------- */ -double PairLCBOP::bondorder(int i, int j, double rij[3], - double rijmag, double VA, - double **f, int vflag_atom) +double PairLCBOP::bondorder(int i, int j, double rij[3],double rijmag, double VA,double **f) { double bij, bji; /* bij & bji */{ double rji[3]; rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; - bij = b(i,j,rij,rijmag,VA,f,vflag_atom); - bji = b(j,i,rji,rijmag,VA,f,vflag_atom); + bij = b(i,j,rij,rijmag,VA,f); + bji = b(j,i,rji,rijmag,VA,f); } double Fij_conj; @@ -663,31 +659,30 @@ double PairLCBOP::bondorder(int i, int j, double rij[3], } double dF_dNij, dF_dNji, dF_dNconj; - Fij_conj = F_conj( Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj ); + Fij_conj = F_conj(Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj); /*forces for Nij*/ if (3-Nij > TOL) { - double factor = -VA*0.5*( dF_dNij + dF_dNconj*( dNconj_dNij + dNconj_dNel*dNij_el_dNij ) ); - FNij( i, j, factor, f, vflag_atom ); + double factor = -VA*0.5*(dF_dNij + dF_dNconj*(dNconj_dNij + dNconj_dNel*dNij_el_dNij)); + FNij(i, j, factor, f); } /*forces for Nji*/ if (3-Nji > TOL) { - double factor = -VA*0.5*( dF_dNji + dF_dNconj*( dNconj_dNji + dNconj_dNel*dNji_el_dNji ) ); - FNij( j, i, factor, f, vflag_atom ); + double factor = -VA*0.5*(dF_dNji + dF_dNconj*(dNconj_dNji + dNconj_dNel*dNji_el_dNji)); + FNij(j, i, factor, f); } /*forces for Mij*/ if (3-Mij > TOL) { - double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNij_el_dMij ); - FMij( i, j, factor, f, vflag_atom ); + double factor = -VA*0.5*(dF_dNconj*dNconj_dNel*dNij_el_dMij); + FMij(i, j, factor, f); } if (3-Mji > TOL) { - double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNji_el_dMji ); - FMij( j, i, factor, f, vflag_atom ); + double factor = -VA*0.5*(dF_dNconj*dNconj_dNel*dNji_el_dMji); + FMij(j,i,factor,f); } } - - double Bij = 0.5*( bij + bji + Fij_conj ); + double Bij = 0.5*(bij + bji + Fij_conj); return Bij; } @@ -695,9 +690,7 @@ double PairLCBOP::bondorder(int i, int j, double rij[3], bij function ------------------------------------------------------------------------- */ -double PairLCBOP::b(int i, int j, double rij[3], - double rijmag, double VA, - double **f, int vflag_atom) { +double PairLCBOP::b(int i, int j, double rij[3], double rijmag, double VA, double **f) { int *SR_neighs = SR_firstneigh[i]; double **x = atom->x; int atomi = i; @@ -818,7 +811,7 @@ double PairLCBOP::b(int i, int j, double rij[3], f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; - if (vflag_atom) { + if (vflag_either) { double rji[3], rki[3]; rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; diff --git a/src/MANYBODY/pair_lcbop.h b/src/MANYBODY/pair_lcbop.h index 9350283f18..836e66e160 100644 --- a/src/MANYBODY/pair_lcbop.h +++ b/src/MANYBODY/pair_lcbop.h @@ -76,10 +76,10 @@ class PairLCBOP : public Pair { void FSR(int, int); void FLR(int, int); - void FNij(int, int, double, double **, int); - void FMij(int, int, double, double **, int); - double bondorder(int, int, double *, double, double, double **, int); - double b(int, int, double *, double, double, double **, int); + void FNij(int, int, double, double **); + void FMij(int, int, double, double **); + double bondorder(int, int, double *, double, double, double **); + double b(int, int, double *, double, double, double **); double gSpline(double, double *); double hSpline(double, double *); diff --git a/src/MANYBODY/pair_polymorphic.cpp b/src/MANYBODY/pair_polymorphic.cpp index 319f0dde8f..aa269ae6ea 100644 --- a/src/MANYBODY/pair_polymorphic.cpp +++ b/src/MANYBODY/pair_polymorphic.cpp @@ -293,8 +293,7 @@ void PairPolymorphic::compute(int eflag, int vflag) f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; - 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); } if (eta == 1) { @@ -352,7 +351,7 @@ void PairPolymorphic::compute(int eflag, int vflag) f[k][1] -= delr2[1]*fpair; f[k][2] -= delr2[2]*fpair; - if (vflag_atom) v_tally2(i, k, -fpair, delr2); + if (vflag_either) v_tally2(i, k, -fpair, delr2); } } @@ -419,8 +418,8 @@ void PairPolymorphic::compute(int eflag, int vflag) f[j][1] -= delr1[1]*fpair; f[j][2] -= delr1[2]*fpair; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]); + if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0, + -fpair,-delr1[0],-delr1[1],-delr1[2]); // attractive term via loop over k @@ -451,7 +450,7 @@ void PairPolymorphic::compute(int eflag, int vflag) f[k][1] += fk[1]; f[k][2] += fk[2]; - if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2); } } } diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 33682ce3fa..f21c3bffca 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -88,10 +88,10 @@ void PairTersoff::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(); + if (vflag_either) eval<1,1,1,1>(); else eval<1,1,1,0>(); } else { - if (vflag_atom) eval<1,1,0,1>(); + if (vflag_either) eval<1,1,0,1>(); else eval<1,1,0,0>(); } } else eval<1,0,0,0>(); @@ -100,17 +100,17 @@ void PairTersoff::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(); + if (vflag_either) eval<0,1,1,1>(); else eval<0,1,1,0>(); } else { - if (vflag_atom) eval<0,1,0,1>(); + if (vflag_either) eval<0,1,0,1>(); else eval<0,1,0,0>(); } } else eval<0,0,0,0>(); } } -template +template void PairTersoff::eval() { int i,j,k,ii,jj,kk,inum,jnum; @@ -315,7 +315,7 @@ void PairTersoff::eval() f[k][1] += fk[1]; f[k][2] += fk[2]; - if (VFLAG_ATOM) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (VFLAG_EITHER) v_tally3(i,j,k,fj,fk,delr1,delr2); } f[j][0] += fjxtmp; f[j][1] += fjytmp; diff --git a/src/MEAM/meam.h b/src/MEAM/meam.h index d41e5f6aed..5205e01ec8 100644 --- a/src/MEAM/meam.h +++ b/src/MEAM/meam.h @@ -294,10 +294,10 @@ class MEAM { void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double *eng_vdwl, double *eatom, int ntype, int *type, int *fmap, double **scale, int &errorflag); - void meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, + void meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom, double *eng_vdwl, double *eatom, int ntype, int *type, int *fmap, double **scale, double **x, int numneigh, int *firstneigh, int numneigh_full, - int *firstneigh_full, int fnoffset, double **f, double **vatom); + int *firstneigh_full, int fnoffset, double **f, double **vatom, double *virial); }; // Functions we need for compat diff --git a/src/MEAM/meam_force.cpp b/src/MEAM/meam_force.cpp index ab72d3c6c8..88b6140f80 100644 --- a/src/MEAM/meam_force.cpp +++ b/src/MEAM/meam_force.cpp @@ -21,12 +21,15 @@ using namespace LAMMPS_NS; void -MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl, - double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh, - int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom) +MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom, + double* eng_vdwl, double* eatom, int /*ntype*/, int* type, int* fmap, + double** scale, double** x, int numneigh, int* firstneigh, int numneigh_full, + int* firstneigh_full, int fnoffset, double** f, double** vatom, double *virial) { int j, jn, k, kn, kk, m, n, p, q; int nv2, nv3, elti, eltj, eltk, ind; + int eflag_either = eflag_atom || eflag_global; + int vflag_either = vflag_atom || vflag_global; double xitmp, yitmp, zitmp, delij[3], rij2, rij, rij3; double v[6], fi[3], fj[3]; double third, sixth; @@ -414,7 +417,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int // Tabulate per-atom virial as symmetrized stress tensor - if (vflag_atom != 0) { + if (vflag_either) { fi[0] = delij[0] * force + dUdrijm[0]; fi[1] = delij[1] * force + dUdrijm[1]; fi[2] = delij[2] * force + dUdrijm[2]; @@ -425,9 +428,16 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int v[4] = -0.25 * (delij[0] * fi[2] + delij[2] * fi[0]); v[5] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]); - for (m = 0; m < 6; m++) { - vatom[i][m] = vatom[i][m] + v[m]; - vatom[j][m] = vatom[j][m] + v[m]; + if (vflag_global) { + for (m = 0; m < 6; m++) { + virial[m] += 2.0*v[m]; + } + } + if (vflag_atom) { + for (m = 0; m < 6; m++) { + vatom[i][m] += v[m]; + vatom[j][m] += v[m]; + } } } @@ -499,7 +509,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int // Tabulate per-atom virial as symmetrized stress tensor - if (vflag_atom != 0) { + if (vflag_either) { fi[0] = force1 * dxik; fi[1] = force1 * dyik; fi[2] = force1 * dzik; @@ -513,10 +523,18 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int v[4] = -sixth * (dxik * fi[2] + dxjk * fj[2] + dzik * fi[0] + dzjk * fj[0]); v[5] = -sixth * (dyik * fi[2] + dyjk * fj[2] + dzik * fi[1] + dzjk * fj[1]); - for (m = 0; m < 6; m++) { - vatom[i][m] = vatom[i][m] + v[m]; - vatom[j][m] = vatom[j][m] + v[m]; - vatom[k][m] = vatom[k][m] + v[m]; + if (vflag_global) { + for (m = 0; m < 6; m++) { + virial[m] += 3.0*v[m]; + } + } + + if (vflag_atom) { + for (m = 0; m < 6; m++) { + vatom[i][m] += v[m]; + vatom[j][m] += v[m]; + vatom[k][m] += v[m]; + } } } } diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index 0d607b8fb2..2acf58f738 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -154,17 +154,16 @@ void PairMEAM::compute(int eflag, int vflag) // vptr is first value in vatom if it will be used by meam_force() // else vatom may not exist, so pass dummy ptr - double **vptr; + double **vptr = nullptr; if (vflag_atom) vptr = vatom; - else vptr = nullptr; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; - meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom, - vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x, - numneigh_half[i],firstneigh_half[i], - numneigh_full[i],firstneigh_full[i], - offset,f,vptr); + meam_inst->meam_force(i,eflag_global,eflag_atom,vflag_global, + vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x, + numneigh_half[i],firstneigh_half[i], + numneigh_full[i],firstneigh_full[i], + offset,f,vptr,virial); offset += numneigh_half[i]; } diff --git a/src/ML-RANN/pair_rann.cpp b/src/ML-RANN/pair_rann.cpp index 2fe8f4d0c1..f0d1dec876 100644 --- a/src/ML-RANN/pair_rann.cpp +++ b/src/ML-RANN/pair_rann.cpp @@ -701,15 +701,22 @@ bool PairRANN::check_potential() { void PairRANN::compute(int eflag, int vflag) { - //perform force/energy computation_ + //perform force/energy computation_ if (dospin) { if (strcmp(update->unit_style,"metal") != 0) error->one(FLERR,"Spin pair styles require metal units"); if (!atom->sp_flag) error->one(FLERR,"Spin pair styles requires atom/spin style"); } - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = vflag_atom = 0; + ev_init(eflag,vflag); + + // only global virial via fdotr is supported by this pair style + + if (vflag_atom) + error->all(FLERR,"Pair style rann does not support computing per-atom stress"); + if (vflag && !vflag_fdotr) + error->all(FLERR,"Pair style rann does not support 'pair_modify nofdotr'"); + int ii,i,j; int nn = 0; sims = new Simulation[1]; @@ -724,13 +731,11 @@ void PairRANN::compute(int eflag, int vflag) sims->s = atom->sp; } int itype,f,jnum,len; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; if (eflag_global) {eng_vdwl=0;eng_coul=0;} double energy=0; double **force = atom->f; double **fm = atom->fm; - double **virial = vatom; + //loop over atoms for (ii=0;iiinum;ii++) { i = sims->ilist[ii]; @@ -800,12 +805,8 @@ void PairRANN::compute(int eflag, int vflag) sx[j]=sy[j]=sz[j]=0; } } - if (doscreen) { - screen(ii,0,jnum-1); - } - if (allscreen) { - screen_neighbor_list(&jnum); - } + if (doscreen) screening(ii,0,jnum-1); + if (allscreen) screen_neighbor_list(&jnum); //do fingerprints for atom type len = fingerprintperelement[itype]; for (j=0;jilist; @@ -1026,8 +1026,8 @@ void PairRANN::propagateforward(double * energy,double **force,double ** /*viria sum[j] = activation[itype][i]->activation_function(sum[j]); if (i==L-1) { energy[j] = sum[j]; - if (eflag_atom)eatom[i1]=sum[j]; - if (eflag_global) {eng_vdwl +=sum[j];} + if (eflag_atom) eatom[i1]=sum[j]; + if (eflag_global) eng_vdwl +=sum[j]; } //force propagation for (jj=0;jjilist; @@ -1103,8 +1103,8 @@ void PairRANN::propagateforwardspin(double * energy,double **force,double **fm,d sum[j] = activation[itype][i]->activation_function(sum[j]); if (i==L-1) { energy[j] = sum[j]; - if (eflag_atom)eatom[i1]=sum[j]; - if (eflag_global) {eng_vdwl +=sum[j];} + if (eflag_atom) eatom[i1]=sum[j]; + if (eflag_global) eng_vdwl +=sum[j]; } //force propagation for (jj=0;jjrequest(this,instance_me); diff --git a/src/ML-RANN/pair_rann.h b/src/ML-RANN/pair_rann.h index a89e0265b6..ceaa5f445c 100644 --- a/src/ML-RANN/pair_rann.h +++ b/src/ML-RANN/pair_rann.h @@ -149,9 +149,9 @@ namespace LAMMPS_NS { void read_screening(std::vector,std::vector,char*,int); void read_mass(const std::vector &, const std::vector &,const char*,int); bool check_potential();//after finishing reading potential file - void propagateforward(double *,double **,double **,int,int);//called by compute to get force and energy - void propagateforwardspin(double *,double **,double **,double**,int,int);//called by compute to get force and energy - void screen(int,int,int); + void propagateforward(double *,double **,int,int);//called by compute to get force and energy + void propagateforwardspin(double *,double **,double**,int,int);//called by compute to get force and energy + void screening(int,int,int); void cull_neighbor_list(int *,int,int); void screen_neighbor_list(int *); }; diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index b8d8333902..dd023e096c 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -70,6 +70,7 @@ PairSNAP::~PairSNAP() if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(scale); } } @@ -164,6 +165,7 @@ void PairSNAP::compute(int eflag, int vflag) // for neighbors of I within cutoff: // compute Fij = dEi/dRj = -dEi/dRi // add to Fi, subtract from Fj + // scaling is that for type I snaptr->compute_yi(beta[ii]); @@ -178,12 +180,12 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->compute_deidrj(fij); - f[i][0] += fij[0]; - f[i][1] += fij[1]; - f[i][2] += fij[2]; - f[j][0] -= fij[0]; - f[j][1] -= fij[1]; - f[j][2] -= fij[2]; + f[i][0] += fij[0]*scale[itype][itype]; + f[i][1] += fij[1]*scale[itype][itype]; + f[i][2] += fij[2]*scale[itype][itype]; + f[j][0] -= fij[0]*scale[itype][itype]; + f[j][1] -= fij[1]*scale[itype][itype]; + f[j][2] -= fij[2]*scale[itype][itype]; // tally per-atom virial contribution @@ -224,6 +226,7 @@ void PairSNAP::compute(int eflag, int vflag) } } } + evdwl *= scale[itype][itype]; ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0); } @@ -351,9 +354,9 @@ void PairSNAP::allocate() { allocated = 1; int n = atom->ntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(scale,n+1,n+1,"pair:scale"); map = new int[n+1]; } @@ -408,11 +411,16 @@ void PairSNAP::coeff(int narg, char **arg) } // Calculate maximum cutoff for all elements - rcutmax = 0.0; for (int ielem = 0; ielem < nelements; ielem++) rcutmax = MAX(2.0*radelem[ielem]*rcutfac,rcutmax); + // set default scaling + int n = atom->ntypes; + for (int ii = 0; ii < n+1; ii++) + for (int jj = 0; jj < n+1; jj++) + scale[ii][jj] = 1.0; + } /* ---------------------------------------------------------------------- @@ -441,6 +449,7 @@ void PairSNAP::init_style() double PairSNAP::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + scale[j][i] = scale[i][j]; return (radelem[map[i]] + radelem[map[j]])*rcutfac; } @@ -711,6 +720,7 @@ double PairSNAP::memory_usage() int n = atom->ntypes+1; bytes += (double)n*n*sizeof(int); // setflag bytes += (double)n*n*sizeof(double); // cutsq + bytes += (double)n*n*sizeof(double); // scale bytes += (double)n*sizeof(int); // map bytes += (double)beta_max*ncoeff*sizeof(double); // bispectrum bytes += (double)beta_max*ncoeff*sizeof(double); // beta @@ -720,3 +730,9 @@ double PairSNAP::memory_usage() return bytes; } +void *PairSNAP::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + return nullptr; +} diff --git a/src/ML-SNAP/pair_snap.h b/src/ML-SNAP/pair_snap.h index 0cc38ef9cb..afce4bf895 100644 --- a/src/ML-SNAP/pair_snap.h +++ b/src/ML-SNAP/pair_snap.h @@ -34,6 +34,7 @@ class PairSNAP : public Pair { virtual void init_style(); virtual double init_one(int, int); virtual double memory_usage(); + virtual void *extract(const char *, int &); double rcutfac, quadraticflag; // declared public to workaround gcc 4.9 int ncoeff; // compiler bug, manifest in KOKKOS package @@ -55,6 +56,7 @@ class PairSNAP : public Pair { double **coeffelem; // element bispectrum coefficients double **beta; // betas for all atoms in list double **bispectrum; // bispectrum components for all atoms in list + double **scale; // for thermodynamic integration int twojmax, switchflag, bzeroflag, bnormflag; int chemflag, wselfallflag; int chunksize; diff --git a/src/OPENMP/pair_airebo_omp.cpp b/src/OPENMP/pair_airebo_omp.cpp index 2e3b35b799..6e7491787a 100644 --- a/src/OPENMP/pair_airebo_omp.cpp +++ b/src/OPENMP/pair_airebo_omp.cpp @@ -70,9 +70,9 @@ void PairAIREBOOMP::compute(int eflag, int vflag) thr->timer(Timer::START); ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); - FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv0,thr); - if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv1,thr); - if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,&pv2,thr); + FREBO_thr(ifrom,ito,eflag,&pv0,thr); + if (ljflag) FLJ_thr(ifrom,ito,eflag,&pv1,thr); + if (torflag) TORSION_thr(ifrom,ito,eflag,&pv2,thr); thr->timer(Timer::PAIR); reduce_thr(this, eflag, vflag, thr); @@ -184,8 +184,7 @@ void PairAIREBOOMP::REBO_neigh_thr() REBO forces and energy ------------------------------------------------------------------------- */ -void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag, - int vflag_atom, double *pv0, ThrData * const thr) +void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData * const thr) { int i,j,k,m,ii,itype,jtype; tagint itag,jtag; @@ -259,7 +258,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag, del[0] = delx; del[1] = dely; del[2] = delz; - bij = bondorder_thr(i,j,del,rij,VA,vflag_atom,thr); + bij = bondorder_thr(i,j,del,rij,VA,thr); dVAdi = bij*dVA; fpair = -(dVRdi+dVAdi) / rij; @@ -282,8 +281,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag, find 3- and 4-step paths between atoms I,J via REBO neighbor lists ------------------------------------------------------------------------- */ -void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, - int vflag_atom, double *pv1, ThrData * const thr) +void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData * const thr) { int i,j,k,m,ii,jj,kk,mm,jnum,itype,jtype,ktype,mtype; tagint itag,jtag; @@ -547,8 +545,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, delscale[0] = scale * delij[0]; delscale[1] = scale * delij[1]; delscale[2] = scale * delij[2]; - Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA, - delij,rij,vflag_atom,thr); + Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,thr); } else Stb = 0.0; fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) + @@ -576,7 +573,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, f[atomj][1] -= delij[1]*fpair; f[atomj][2] -= delij[2]*fpair; - if (vflag_atom) v_tally2_thr(atomi,atomj,fpair,delij,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomj,fpair,delij,thr); } else if (npath == 3) { fpair1 = dC*dwikS*wkjS / rikS; @@ -598,8 +595,8 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, f[atomk][1] -= fi[1] + fj[1]; f[atomk][2] -= fi[2] + fj[2]; - if (vflag_atom) - v_tally3_thr(atomi,atomj,atomk,fi,fj,delikS,deljkS,thr); + if (vflag_either) + v_tally3_thr(this,atomi,atomj,atomk,fi,fj,delikS,deljkS,thr); } else if (npath == 4) { fpair1 = dC*dwikS*wkmS*wmjS / rikS; @@ -634,11 +631,11 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, f[atomm][1] += fm[1]; f[atomm][2] += fm[2]; - if (vflag_atom) { + if (vflag_either) { delimS[0] = delikS[0] + delkmS[0]; delimS[1] = delikS[1] + delkmS[1]; delimS[2] = delikS[2] + delkmS[2]; - v_tally4_thr(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr); + v_tally4_thr(this,atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr); } } } @@ -650,8 +647,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, torsional forces and energy ------------------------------------------------------------------------- */ -void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag, - double *pv2, ThrData * const thr) +void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData * const thr) { int i,j,k,l,ii; tagint itag,jtag; @@ -1011,7 +1007,7 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag, ------------------------------------------------------------------------- */ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, - double VA, int vflag_atom, ThrData * const thr) + double VA, ThrData * const thr) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; @@ -1197,10 +1193,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; - if (vflag_atom) { + if (vflag_either) { rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; - v_tally3_thr(atomi,atomj,atomk,fj,fk,rji,rki,thr); + v_tally3_thr(this,atomi,atomj,atomk,fj,fk,rji,rki,thr); } } } @@ -1342,9 +1338,9 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; - if (vflag_atom) { + if (vflag_either) { rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; - v_tally3_thr(atomi,atomj,atoml,fi,fl,rij,rlj,thr); + v_tally3_thr(this,atomi,atomj,atoml,fi,fl,rij,rlj,thr); } } } @@ -1378,7 +1374,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -1392,7 +1388,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -1414,7 +1410,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr); + if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr); } } } @@ -1445,7 +1441,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -1459,7 +1455,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -1481,7 +1477,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr); + if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr); } } } @@ -1693,10 +1689,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; f[atom4][2] += f4[2]; - if (vflag_atom) { + if (vflag_either) { r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; - v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr); + v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr); } } } @@ -1729,7 +1725,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -1743,7 +1739,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -1765,7 +1761,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr); + if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr); } } } @@ -1796,7 +1792,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -1810,7 +1806,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -1832,7 +1828,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr); + if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr); } } } @@ -1870,8 +1866,7 @@ there probably also need to be performed here. */ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], double rijmag_mod, - double VA, double rij[3], double rijmag, - int vflag_atom, ThrData * const thr) + double VA, double rij[3], double rijmag, ThrData * const thr) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; @@ -2196,10 +2191,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; - if (vflag_atom) { + if (vflag_either) { rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; - v_tally3_thr(atomi,atomj,atomk,fj,fk,rji,rki,thr); + v_tally3_thr(this,atomi,atomj,atomk,fj,fk,rji,rki,thr); } } } @@ -2300,9 +2295,9 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; - if (vflag_atom) { + if (vflag_either) { rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; - v_tally3_thr(atomi,atomj,atoml,fi,fl,rij,rlj,thr); + v_tally3_thr(this,atomi,atomj,atoml,fi,fl,rij,rlj,thr); } } } @@ -2335,7 +2330,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -2349,7 +2344,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -2371,7 +2366,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr); + if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr); } } } @@ -2402,7 +2397,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -2416,7 +2411,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -2438,7 +2433,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr); + if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr); } } } @@ -2643,10 +2638,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; f[atom4][2] += f4[2]; - if (vflag_atom) { + if (vflag_either) { r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; - v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr); + v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr); } } } @@ -2677,7 +2672,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); // due to kronecker(ktype, 0) term in contribution // to NconjtmpI and later Nijconj @@ -2691,7 +2686,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; - if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr); if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; @@ -2713,7 +2708,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; - if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr); + if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr); } } } @@ -2744,7 +2739,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); // due to kronecker(ltype, 0) term in contribution // to NconjtmpJ and later Nijconj @@ -2758,7 +2753,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; - if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr); if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; @@ -2780,7 +2775,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; - if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr); + if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr); } } } diff --git a/src/OPENMP/pair_airebo_omp.h b/src/OPENMP/pair_airebo_omp.h index 9c5b29545f..355614552e 100644 --- a/src/OPENMP/pair_airebo_omp.h +++ b/src/OPENMP/pair_airebo_omp.h @@ -33,16 +33,13 @@ class PairAIREBOOMP : public PairAIREBO, public ThrOMP { virtual double memory_usage(); protected: - double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, int vflag_atom, - ThrData *const thr); + double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, ThrData *const thr); double bondorderLJ_thr(int i, int j, double rij[3], double rijmag, double VA, double rij0[3], - double rijmag0, int vflag_atom, ThrData *const thr); + double rijmag0, ThrData *const thr); - void FREBO_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv0, - ThrData *const thr); - void FLJ_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv1, - ThrData *const thr); - void TORSION_thr(int ifrom, int ito, int evflag, int eflag, double *pv2, ThrData *const thr); + void FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData *const thr); + void FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData *const thr); + void TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData *const thr); void REBO_neigh_thr(); }; diff --git a/src/OPENMP/pair_comb_omp.cpp b/src/OPENMP/pair_comb_omp.cpp index 7c0c14f7e5..169704b289 100644 --- a/src/OPENMP/pair_comb_omp.cpp +++ b/src/OPENMP/pair_comb_omp.cpp @@ -70,10 +70,10 @@ void PairCombOMP::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1>(ifrom, ito, thr); else eval<1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,0,1>(ifrom, ito, thr); else eval<1,0,0>(ifrom, ito, thr); } } else eval<0,0,0>(ifrom, ito, thr); @@ -83,7 +83,7 @@ void PairCombOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum,iparam_i; @@ -365,7 +365,7 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1, elp_ij,0.0,0.0,0.0,0.0,0.0, thr); - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j][0] += fjxtmp; f[j][1] += fjytmp; diff --git a/src/OPENMP/pair_edip_omp.cpp b/src/OPENMP/pair_edip_omp.cpp index a23764f959..da025d75d5 100644 --- a/src/OPENMP/pair_edip_omp.cpp +++ b/src/OPENMP/pair_edip_omp.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -19,22 +18,22 @@ #include "comm.h" #include "neigh_list.h" #include "suffix.h" -using namespace LAMMPS_NS; #include #include "omp_compat.h" -#define GRIDDENSITY 8000 -#define GRIDSTART 0.1 +using namespace LAMMPS_NS; // max number of interaction per atom for f(Z) environment potential static constexpr int leadDimInteractionList = 64; +#define GRIDDENSITY 8000 +#define GRIDSTART 0.1 + /* ---------------------------------------------------------------------- */ -PairEDIPOMP::PairEDIPOMP(LAMMPS *lmp) : - PairEDIP(lmp), ThrOMP(lmp, THR_PAIR) +PairEDIPOMP::PairEDIPOMP(LAMMPS *lmp) : PairEDIP(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; respa_enable = 0; @@ -44,14 +43,14 @@ PairEDIPOMP::PairEDIPOMP(LAMMPS *lmp) : void PairEDIPOMP::compute(int eflag, int vflag) { - ev_init(eflag,vflag); + ev_init(eflag, vflag); const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; const int inum = list->inum; #if defined(_OPENMP) -#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) #endif { int ifrom, ito, tid; @@ -63,26 +62,31 @@ void PairEDIPOMP::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1>(ifrom, ito, thr); - else eval<1,1,0>(ifrom, ito, thr); + if (vflag_atom) + eval<1, 1, 1>(ifrom, ito, thr); + else + eval<1, 1, 0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,0,1>(ifrom, ito, thr); - else eval<1,0,0>(ifrom, ito, thr); + if (vflag_atom) + eval<1, 0, 1>(ifrom, ito, thr); + else + eval<1, 0, 0>(ifrom, ito, thr); } - } else eval<0,0,0>(ifrom, ito, thr); + } else + eval<0, 0, 0>(ifrom, ito, thr); thr->timer(Timer::PAIR); reduce_thr(this, eflag, vflag, thr); - } // end of omp parallel region + } // end of omp parallel region } template -void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) +void PairEDIPOMP::eval(int iifrom, int iito, ThrData *const thr) { - int i,j,k,ii,jnum; - int itype,jtype,ktype,ijparam,ikparam; - double xtmp,ytmp,ztmp,evdwl; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, k, ii, jnum; + int itype, jtype, ktype, ijparam, ikparam; + double xtmp, ytmp, ztmp, evdwl; + int *ilist, *jlist, *numneigh, **firstneigh; int preForceCoord_counter; double invR_ij; @@ -147,9 +151,9 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) double *pre_thrPow2B_ij = prePow2B_ij + tid * leadDimInteractionList; double *pre_thrForceCoord = preForceCoord + tid * leadDimInteractionList; - const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; - dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; - const int * _noalias const type = atom->type; + const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0]; + const int *_noalias const type = atom->type; const int nlocal = atom->nlocal; ilist = list->ilist; @@ -174,86 +178,82 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) // pre-loop to compute environment coordination f(Z) for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { - j = jlist[neighbor_j]; - j &= NEIGHMASK; + j = jlist[neighbor_j]; + j &= NEIGHMASK; - double dr_ij[3], r_ij; + double dr_ij[3], r_ij; - dr_ij[0] = xtmp - x[j].x; - dr_ij[1] = ytmp - x[j].y; - dr_ij[2] = ztmp - x[j].z; - r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + dr_ij[0] = xtmp - x[j].x; + dr_ij[1] = ytmp - x[j].y; + dr_ij[2] = ztmp - x[j].z; + r_ij = dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2]; - jtype = map[type[j]]; - ijparam = elem3param[itype][jtype][jtype]; - if (r_ij > params[ijparam].cutsq) continue; + jtype = map[type[j]]; + ijparam = elem3param[itype][jtype][jtype]; + if (r_ij > params[ijparam].cutsq) continue; - r_ij = sqrt(r_ij); + r_ij = sqrt(r_ij); - invR_ij = 1.0 / r_ij; - pre_thrInvR_ij[neighbor_j] = invR_ij; + invR_ij = 1.0 / r_ij; + pre_thrInvR_ij[neighbor_j] = invR_ij; - invRMinusCutoffA = 1.0 / (r_ij - cutoffA); - sigmaInvRMinusCutoffA = sigma * invRMinusCutoffA; - gammInvRMinusCutoffA = gamm * invRMinusCutoffA; + invRMinusCutoffA = 1.0 / (r_ij - cutoffA); + sigmaInvRMinusCutoffA = sigma * invRMinusCutoffA; + gammInvRMinusCutoffA = gamm * invRMinusCutoffA; - interpolDeltaX = r_ij - GRIDSTART; - interpolTMP = (interpolDeltaX * GRIDDENSITY); - interpolIDX = (int) interpolTMP; + interpolDeltaX = r_ij - GRIDSTART; + interpolTMP = (interpolDeltaX * GRIDDENSITY); + interpolIDX = (int) interpolTMP; - interpolY1 = exp3B[interpolIDX]; - interpolY2 = exp3B[interpolIDX+1]; - exp3B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = exp3B[interpolIDX]; + interpolY2 = exp3B[interpolIDX + 1]; + exp3B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - exp3BDerived_ij = - exp3B_ij * gammInvRMinusCutoffA * invRMinusCutoffA; + exp3BDerived_ij = -exp3B_ij * gammInvRMinusCutoffA * invRMinusCutoffA; - pre_thrExp3B_ij[neighbor_j] = exp3B_ij; - pre_thrExp3BDerived_ij[neighbor_j] = exp3BDerived_ij; + pre_thrExp3B_ij[neighbor_j] = exp3B_ij; + pre_thrExp3BDerived_ij[neighbor_j] = exp3BDerived_ij; - interpolY1 = exp2B[interpolIDX]; - interpolY2 = exp2B[interpolIDX+1]; - exp2B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = exp2B[interpolIDX]; + interpolY2 = exp2B[interpolIDX + 1]; + exp2B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - exp2BDerived_ij = - exp2B_ij * sigmaInvRMinusCutoffA * invRMinusCutoffA; + exp2BDerived_ij = -exp2B_ij * sigmaInvRMinusCutoffA * invRMinusCutoffA; - pre_thrExp2B_ij[neighbor_j] = exp2B_ij; - pre_thrExp2BDerived_ij[neighbor_j] = exp2BDerived_ij; + pre_thrExp2B_ij[neighbor_j] = exp2B_ij; + pre_thrExp2BDerived_ij[neighbor_j] = exp2BDerived_ij; - interpolY1 = pow2B[interpolIDX]; - interpolY2 = pow2B[interpolIDX+1]; - pow2B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = pow2B[interpolIDX]; + interpolY2 = pow2B[interpolIDX + 1]; + pow2B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - pre_thrPow2B_ij[neighbor_j] = pow2B_ij; + pre_thrPow2B_ij[neighbor_j] = pow2B_ij; - // zeta and its derivative + // zeta and its derivative - if (r_ij < cutoffC) zeta_i += 1.0; - else { - interpolY1 = cutoffFunction[interpolIDX]; - interpolY2 = cutoffFunction[interpolIDX+1]; - cutoffFunction_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + if (r_ij < cutoffC) + zeta_i += 1.0; + else { + interpolY1 = cutoffFunction[interpolIDX]; + interpolY2 = cutoffFunction[interpolIDX + 1]; + cutoffFunction_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - zeta_i += cutoffFunction_ij; + zeta_i += cutoffFunction_ij; - interpolY1 = cutoffFunctionDerived[interpolIDX]; - interpolY2 = cutoffFunctionDerived[interpolIDX+1]; - zeta_iDerived = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = cutoffFunctionDerived[interpolIDX]; + interpolY2 = cutoffFunctionDerived[interpolIDX + 1]; + zeta_iDerived = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - zeta_iDerivedInvR_ij = zeta_iDerived * invR_ij; + zeta_iDerivedInvR_ij = zeta_iDerived * invR_ij; - preForceCoord_counter=numForceCoordPairs*5; - pre_thrForceCoord[preForceCoord_counter+0]=zeta_iDerivedInvR_ij; - pre_thrForceCoord[preForceCoord_counter+1]=dr_ij[0]; - pre_thrForceCoord[preForceCoord_counter+2]=dr_ij[1]; - pre_thrForceCoord[preForceCoord_counter+3]=dr_ij[2]; - pre_thrForceCoord[preForceCoord_counter+4]=j; - numForceCoordPairs++; - } + preForceCoord_counter = numForceCoordPairs * 5; + pre_thrForceCoord[preForceCoord_counter + 0] = zeta_iDerivedInvR_ij; + pre_thrForceCoord[preForceCoord_counter + 1] = dr_ij[0]; + pre_thrForceCoord[preForceCoord_counter + 2] = dr_ij[1]; + pre_thrForceCoord[preForceCoord_counter + 3] = dr_ij[2]; + pre_thrForceCoord[preForceCoord_counter + 4] = j; + numForceCoordPairs++; + } } // quantities depending on zeta_i @@ -263,24 +263,20 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) interpolIDX = (int) interpolTMP; interpolY1 = expMinusBetaZeta_iZeta_iGrid[interpolIDX]; - interpolY2 = expMinusBetaZeta_iZeta_iGrid[interpolIDX+1]; - expMinusBetaZeta_iZeta_i = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = expMinusBetaZeta_iZeta_iGrid[interpolIDX + 1]; + expMinusBetaZeta_iZeta_i = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = qFunctionGrid[interpolIDX]; - interpolY2 = qFunctionGrid[interpolIDX+1]; - qFunction = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = qFunctionGrid[interpolIDX + 1]; + qFunction = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = tauFunctionGrid[interpolIDX]; - interpolY2 = tauFunctionGrid[interpolIDX+1]; - tauFunction = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = tauFunctionGrid[interpolIDX + 1]; + tauFunction = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = tauFunctionDerivedGrid[interpolIDX]; - interpolY2 = tauFunctionDerivedGrid[interpolIDX+1]; - tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = tauFunctionDerivedGrid[interpolIDX + 1]; + tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); forceModCoord_factor = 2.0 * beta * zeta_i * expMinusBetaZeta_iZeta_i; @@ -297,7 +293,7 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) dr_ij[0] = x[j].x - xtmp; dr_ij[1] = x[j].y - ytmp; dr_ij[2] = x[j].z - ztmp; - r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + r_ij = dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2]; jtype = map[type[j]]; ijparam = elem3param[itype][jtype][jtype]; @@ -312,13 +308,12 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) exp2B_ij = pre_thrExp2B_ij[neighbor_j]; - pow2BDerived_ij = - rho * invR_ij * pow2B_ij; + pow2BDerived_ij = -rho * invR_ij * pow2B_ij; - forceModCoord += (forceModCoord_factor*exp2B_ij); + forceModCoord += (forceModCoord_factor * exp2B_ij); exp2BDerived_ij = pre_thrExp2BDerived_ij[neighbor_j]; - forceMod2B = exp2BDerived_ij * potential2B_factor + - exp2B_ij * pow2BDerived_ij; + forceMod2B = exp2BDerived_ij * potential2B_factor + exp2B_ij * pow2BDerived_ij; directorCos_ij_x = invR_ij * dr_ij[0]; directorCos_ij_y = invR_ij * dr_ij[1]; @@ -343,136 +338,123 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr) evdwl = (exp2B_ij * potential2B_factor); - if (EVFLAG) ev_tally_thr(this,i, j, nlocal, /* newton_pair */ 1, evdwl, 0.0, - -forceMod2B*invR_ij, dr_ij[0], dr_ij[1], dr_ij[2],thr); + if (EVFLAG) + ev_tally_thr(this, i, j, nlocal, /* newton_pair */ 1, evdwl, 0.0, -forceMod2B * invR_ij, + dr_ij[0], dr_ij[1], dr_ij[2], thr); // three-body Forces for (int neighbor_k = neighbor_j + 1; neighbor_k < jnum; neighbor_k++) { - double dr_ik[3], r_ik, f_ik[3]; + double dr_ik[3], r_ik, f_ik[3]; - k = jlist[neighbor_k]; - k &= NEIGHMASK; - ktype = map[type[k]]; - ikparam = elem3param[itype][ktype][ktype]; + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem3param[itype][ktype][ktype]; - dr_ik[0] = x[k].x - xtmp; - dr_ik[1] = x[k].y - ytmp; - dr_ik[2] = x[k].z - ztmp; - r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + dr_ik[0] = x[k].x - xtmp; + dr_ik[1] = x[k].y - ytmp; + dr_ik[2] = x[k].z - ztmp; + r_ik = dr_ik[0] * dr_ik[0] + dr_ik[1] * dr_ik[1] + dr_ik[2] * dr_ik[2]; - if (r_ik > params[ikparam].cutsq) continue; + if (r_ik > params[ikparam].cutsq) continue; - r_ik = sqrt(r_ik); + r_ik = sqrt(r_ik); - invR_ik = pre_thrInvR_ij[neighbor_k]; + invR_ik = pre_thrInvR_ij[neighbor_k]; - directorCos_ik_x = invR_ik * dr_ik[0]; - directorCos_ik_y = invR_ik * dr_ik[1]; - directorCos_ik_z = invR_ik * dr_ik[2]; + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; - cosTeta = directorCos_ij_x * directorCos_ik_x + - directorCos_ij_y * directorCos_ik_y + + cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + directorCos_ij_z * directorCos_ik_z; - cosTetaDiff = cosTeta + tauFunction; - cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; - qFunctionCosTetaDiffCosTetaDiff = cosTetaDiffCosTetaDiff * qFunction; - expMinusQFunctionCosTetaDiffCosTetaDiff = - exp(-qFunctionCosTetaDiffCosTetaDiff); + cosTetaDiff = cosTeta + tauFunction; + cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; + qFunctionCosTetaDiffCosTetaDiff = cosTetaDiffCosTetaDiff * qFunction; + expMinusQFunctionCosTetaDiffCosTetaDiff = exp(-qFunctionCosTetaDiffCosTetaDiff); - potentia3B_factor = lambda * + potentia3B_factor = lambda * ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + eta * qFunctionCosTetaDiffCosTetaDiff); - exp3B_ik = pre_thrExp3B_ij[neighbor_k]; - exp3BDerived_ik = pre_thrExp3BDerived_ij[neighbor_k]; + exp3B_ik = pre_thrExp3B_ij[neighbor_k]; + exp3BDerived_ik = pre_thrExp3BDerived_ij[neighbor_k]; - forceMod3B_factor1_ij = - exp3BDerived_ij * exp3B_ik * - potentia3B_factor; - forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * - qFunction * cosTetaDiff * + forceMod3B_factor1_ij = -exp3BDerived_ij * exp3B_ik * potentia3B_factor; + forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff * (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); - forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; + forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; - f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_x - directorCos_ik_x); - f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_y - directorCos_ik_y); - f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_z - directorCos_ik_z); + f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); + f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); + f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); - forceMod3B_factor1_ik = - exp3BDerived_ik * exp3B_ij * - potentia3B_factor; - forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; + forceMod3B_factor1_ik = -exp3BDerived_ik * exp3B_ij * potentia3B_factor; + forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; - f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_x - directorCos_ij_x); - f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_y - directorCos_ij_y); - f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_z - directorCos_ij_z); + f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); + f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); + f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); - forceModCoord += (forceMod3B_factor2 * - (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); + forceModCoord += (forceMod3B_factor2 * (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); - f[j].x += f_ij[0]; - f[j].y += f_ij[1]; - f[j].z += f_ij[2]; + f[j].x += f_ij[0]; + f[j].y += f_ij[1]; + f[j].z += f_ij[2]; - f[k].x += f_ik[0]; - f[k].y += f_ik[1]; - f[k].z += f_ik[2]; + f[k].x += f_ik[0]; + f[k].y += f_ik[1]; + f[k].z += f_ik[2]; - f[i].x -= f_ij[0] + f_ik[0]; - f[i].y -= f_ij[1] + f_ik[1]; - f[i].z -= f_ij[2] + f_ik[2]; + f[i].x -= f_ij[0] + f_ik[0]; + f[i].y -= f_ij[1] + f_ik[1]; + f[i].z -= f_ij[2] + f_ik[2]; - // potential energy + // potential energy - evdwl = (exp3B_ij * exp3B_ik * potentia3B_factor); + evdwl = (exp3B_ij * exp3B_ik * potentia3B_factor); - if (evflag) ev_tally3_thr(this,i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik,thr); + if (EVFLAG) ev_tally3_thr(this, i, j, k, evdwl, 0.0, f_ij, f_ik, dr_ij, dr_ik, thr); } } // forces due to environment coordination f(Z) for (int idx = 0; idx < numForceCoordPairs; idx++) { - double dr_ij[3], f_ij[3]; + double dr_ij[3], f_ij[3]; - preForceCoord_counter = idx * 5; - zeta_iDerivedInvR_ij=pre_thrForceCoord[preForceCoord_counter+0]; - dr_ij[0]=pre_thrForceCoord[preForceCoord_counter+1]; - dr_ij[1]=pre_thrForceCoord[preForceCoord_counter+2]; - dr_ij[2]=pre_thrForceCoord[preForceCoord_counter+3]; - j = static_cast (pre_thrForceCoord[preForceCoord_counter+4]); + preForceCoord_counter = idx * 5; + zeta_iDerivedInvR_ij = pre_thrForceCoord[preForceCoord_counter + 0]; + dr_ij[0] = pre_thrForceCoord[preForceCoord_counter + 1]; + dr_ij[1] = pre_thrForceCoord[preForceCoord_counter + 2]; + dr_ij[2] = pre_thrForceCoord[preForceCoord_counter + 3]; + j = static_cast(pre_thrForceCoord[preForceCoord_counter + 4]); - forceModCoord_ij = forceModCoord * zeta_iDerivedInvR_ij; + forceModCoord_ij = forceModCoord * zeta_iDerivedInvR_ij; - f_ij[0] = forceModCoord_ij * dr_ij[0]; - f_ij[1] = forceModCoord_ij * dr_ij[1]; - f_ij[2] = forceModCoord_ij * dr_ij[2]; + f_ij[0] = forceModCoord_ij * dr_ij[0]; + f_ij[1] = forceModCoord_ij * dr_ij[1]; + f_ij[2] = forceModCoord_ij * dr_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]; + f[j].x += f_ij[0]; + f[j].y += f_ij[1]; + f[j].z += f_ij[2]; - // potential energy - - evdwl = 0.0; - if (EVFLAG) ev_tally_thr(this,i, j, nlocal, /* newton_pair */ 1, 0.0, 0.0, - forceModCoord_ij, dr_ij[0], dr_ij[1], dr_ij[2],thr); + if (EVFLAG) + ev_tally_thr(this, i, j, nlocal, /* newton_pair */ 1, 0.0, 0.0, -forceModCoord_ij, dr_ij[0], + dr_ij[1], dr_ij[2], thr); } } } diff --git a/src/OPENMP/pair_tersoff_mod_c_omp.cpp b/src/OPENMP/pair_tersoff_mod_c_omp.cpp index 0d66e83b00..cd59e06673 100644 --- a/src/OPENMP/pair_tersoff_mod_c_omp.cpp +++ b/src/OPENMP/pair_tersoff_mod_c_omp.cpp @@ -60,20 +60,20 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); } else { if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -84,7 +84,7 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum; @@ -282,7 +282,7 @@ void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_mod_omp.cpp b/src/OPENMP/pair_tersoff_mod_omp.cpp index b3e3b89283..de19aa3872 100644 --- a/src/OPENMP/pair_tersoff_mod_omp.cpp +++ b/src/OPENMP/pair_tersoff_mod_omp.cpp @@ -60,20 +60,20 @@ void PairTersoffMODOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); } else { if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -84,7 +84,7 @@ void PairTersoffMODOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum; @@ -282,7 +282,7 @@ void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_omp.cpp b/src/OPENMP/pair_tersoff_omp.cpp index 6bf48792eb..2a3cfb5d54 100644 --- a/src/OPENMP/pair_tersoff_omp.cpp +++ b/src/OPENMP/pair_tersoff_omp.cpp @@ -61,10 +61,10 @@ void PairTersoffOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); @@ -73,10 +73,10 @@ void PairTersoffOMP::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -87,7 +87,7 @@ void PairTersoffOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum,maxshort_thr; @@ -293,7 +293,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_table_omp.cpp b/src/OPENMP/pair_tersoff_table_omp.cpp index a91f361f79..47ee2feca8 100644 --- a/src/OPENMP/pair_tersoff_table_omp.cpp +++ b/src/OPENMP/pair_tersoff_table_omp.cpp @@ -80,7 +80,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag) ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); if (evflag) - if (vflag_atom) eval<1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1>(ifrom, ito, thr); else eval<1,0>(ifrom, ito, thr); else eval<0,0>(ifrom, ito, thr); @@ -89,7 +89,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jnum; @@ -427,7 +427,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); } // second loop over neighbors of atom i except j, forces and virial only - part 2/2 @@ -493,8 +493,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); - + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); } } // loop on J f[i].x += fxtmp; diff --git a/src/OPENMP/thr_omp.cpp b/src/OPENMP/thr_omp.cpp index 4685e116be..545e3bbe88 100644 --- a/src/OPENMP/thr_omp.cpp +++ b/src/OPENMP/thr_omp.cpp @@ -725,18 +725,27 @@ void ThrOMP::ev_tally4_thr(Pair * const pair, const int i, const int j, } } - if (pair->vflag_atom) { - v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); - v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); - v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); - v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); - v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); - v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + if (pair->vflag_either) { + v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + if (pair->vflag_global) v_tally(thr->virial_pair,v); - v_tally(thr->vatom_pair[i],v); - v_tally(thr->vatom_pair[j],v); - v_tally(thr->vatom_pair[k],v); - v_tally(thr->vatom_pair[m],v); + if (pair->vflag_atom) { + v[0] *= 0.25; + v[1] *= 0.25; + v[2] *= 0.25; + v[3] *= 0.25; + v[4] *= 0.25; + v[5] *= 0.25; + v_tally(thr->vatom_pair[i],v); + v_tally(thr->vatom_pair[j],v); + v_tally(thr->vatom_pair[k],v); + v_tally(thr->vatom_pair[m],v); + } } } @@ -1465,20 +1474,29 @@ void ThrOMP::ev_tally_thr(Improper * const imprp, const int i1, const int i2, fpair is magnitude of force on atom I ------------------------------------------------------------------------- */ -void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair, +void ThrOMP::v_tally2_thr(Pair *const pair, const int i, const int j, const double fpair, const double * const drij, ThrData * const thr) { double v[6]; - v[0] = 0.5 * drij[0]*drij[0]*fpair; - v[1] = 0.5 * drij[1]*drij[1]*fpair; - v[2] = 0.5 * drij[2]*drij[2]*fpair; - v[3] = 0.5 * drij[0]*drij[1]*fpair; - v[4] = 0.5 * drij[0]*drij[2]*fpair; - v[5] = 0.5 * drij[1]*drij[2]*fpair; + v[0] = drij[0]*drij[0]*fpair; + v[1] = drij[1]*drij[1]*fpair; + v[2] = drij[2]*drij[2]*fpair; + v[3] = drij[0]*drij[1]*fpair; + v[4] = drij[0]*drij[2]*fpair; + v[5] = drij[1]*drij[2]*fpair; + if (pair->vflag_global) v_tally(thr->virial_pair,v); - v_tally(thr->vatom_pair[i],v); - v_tally(thr->vatom_pair[j],v); + if (pair->vflag_atom) { + v[0] *= 0.5; + v[1] *= 0.5; + v[2] *= 0.5; + v[3] *= 0.5; + v[4] *= 0.5; + v[5] *= 0.5; + v_tally(thr->vatom_pair[i],v); + v_tally(thr->vatom_pair[j],v); + } } /* ---------------------------------------------------------------------- @@ -1486,23 +1504,32 @@ void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair, called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ -void ThrOMP::v_tally3_thr(const int i, const int j, const int k, +void ThrOMP::v_tally3_thr(Pair *const pair, const int i, const int j, const int k, const double * const fi, const double * const fj, const double * const drik, const double * const drjk, ThrData * const thr) { double v[6]; - v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); - v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); - v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); - v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); - v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); - v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); + v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]); + if (pair->vflag_global) v_tally(thr->virial_pair,v); - v_tally(thr->vatom_pair[i],v); - v_tally(thr->vatom_pair[j],v); - v_tally(thr->vatom_pair[k],v); + if (pair->vflag_atom) { + v[0] *= THIRD; + v[1] *= THIRD; + v[2] *= THIRD; + v[3] *= THIRD; + v[4] *= THIRD; + v[5] *= THIRD; + v_tally(thr->vatom_pair[i],v); + v_tally(thr->vatom_pair[j],v); + v_tally(thr->vatom_pair[k],v); + } } /* ---------------------------------------------------------------------- @@ -1510,7 +1537,7 @@ void ThrOMP::v_tally3_thr(const int i, const int j, const int k, called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ -void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m, +void ThrOMP::v_tally4_thr(Pair *const pair, const int i, const int j, const int k, const int m, const double * const fi, const double * const fj, const double * const fk, const double * const drim, const double * const drjm, const double * const drkm, @@ -1518,17 +1545,26 @@ void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m, { double v[6]; - v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); - v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); - v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); - v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); - v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); - v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + if (pair->vflag_global) v_tally(thr->virial_pair,v); - v_tally(thr->vatom_pair[i],v); - v_tally(thr->vatom_pair[j],v); - v_tally(thr->vatom_pair[k],v); - v_tally(thr->vatom_pair[m],v); + if (pair->vflag_atom) { + v[0] *= 0.25; + v[1] *= 0.25; + v[2] *= 0.25; + v[3] *= 0.25; + v[4] *= 0.25; + v[5] *= 0.25; + v_tally(thr->vatom_pair[i],v); + v_tally(thr->vatom_pair[j],v); + v_tally(thr->vatom_pair[k],v); + v_tally(thr->vatom_pair[m],v); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/OPENMP/thr_omp.h b/src/OPENMP/thr_omp.h index 77835122e0..35cd43be5d 100644 --- a/src/OPENMP/thr_omp.h +++ b/src/OPENMP/thr_omp.h @@ -135,12 +135,19 @@ class ThrOMP { void ev_tally_xyz_full_thr(Pair *const, const int, const double, const double, const double, const double, const double, const double, const double, const double, ThrData *const); + void v_tally2_thr(Pair *const, const int, const int, const double, const double *const, + ThrData *const); void ev_tally3_thr(Pair *const, const int, const int, const int, const double, const double, const double *const, const double *const, const double *const, const double *const, ThrData *const); + void v_tally3_thr(Pair *const, const int, const int, const int, const double *const, + const double *const, const double *const, const double *const, ThrData *const); void ev_tally4_thr(Pair *const, const int, const int, const int, const int, const double, const double *const, const double *const, const double *const, const double *const, const double *const, const double *const, ThrData *const); + void v_tally4_thr(Pair *const, const int, const int, const int, const int, const double *const, + const double *const, const double *const, const double *const, + const double *const, const double *const, ThrData *const); // Bond void ev_tally_thr(Bond *const, const int, const int, const int, const int, const double, @@ -169,12 +176,6 @@ class ThrOMP { ThrData *const); // style independent versions - void v_tally2_thr(const int, const int, const double, const double *const, ThrData *const); - void v_tally3_thr(const int, const int, const int, const double *const, const double *const, - const double *const, const double *const, ThrData *const); - void v_tally4_thr(const int, const int, const int, const int, const double *const, - const double *const, const double *const, const double *const, - const double *const, const double *const, ThrData *const); void ev_tally_list_thr(Pair *const, const int, const int *const, const double *const, const double, const double, ThrData *const); }; diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp index 3f80609853..1624d7c908 100644 --- a/src/USER-MISC/pair_edip.cpp +++ b/src/USER-MISC/pair_edip.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -34,8 +33,8 @@ #include "neighbor.h" #include -#include #include +#include using namespace LAMMPS_NS; @@ -52,12 +51,11 @@ static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ PairEDIP::PairEDIP(LAMMPS *lmp) : - Pair(lmp), preInvR_ij(nullptr), preExp3B_ij(nullptr), preExp3BDerived_ij(nullptr), - preExp2B_ij(nullptr), preExp2BDerived_ij(nullptr), prePow2B_ij(nullptr), - preForceCoord(nullptr), cutoffFunction(nullptr), cutoffFunctionDerived(nullptr), - pow2B(nullptr), exp2B(nullptr), exp3B(nullptr), qFunctionGrid(nullptr), - expMinusBetaZeta_iZeta_iGrid(nullptr), tauFunctionGrid(nullptr), - tauFunctionDerivedGrid(nullptr) + Pair(lmp), preInvR_ij(nullptr), preExp3B_ij(nullptr), preExp3BDerived_ij(nullptr), + preExp2B_ij(nullptr), preExp2BDerived_ij(nullptr), prePow2B_ij(nullptr), preForceCoord(nullptr), + cutoffFunction(nullptr), cutoffFunctionDerived(nullptr), pow2B(nullptr), exp2B(nullptr), + exp3B(nullptr), qFunctionGrid(nullptr), expMinusBetaZeta_iZeta_iGrid(nullptr), + tauFunctionGrid(nullptr), tauFunctionDerivedGrid(nullptr) { single_enable = 0; restartinfo = 0; @@ -90,10 +88,10 @@ PairEDIP::~PairEDIP() void PairEDIP::compute(int eflag, int vflag) { - int i,j,k,ii,inum,jnum; - int itype,jtype,ktype,ijparam,ikparam; - double xtmp,ytmp,ztmp,evdwl; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, k, ii, inum, jnum; + int itype, jtype, ktype, ijparam, ikparam; + double xtmp, ytmp, ztmp, evdwl; + int *ilist, *jlist, *numneigh, **firstneigh; int preForceCoord_counter; double invR_ij; @@ -149,7 +147,7 @@ void PairEDIP::compute(int eflag, int vflag) double potential2B_factor; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -180,86 +178,82 @@ void PairEDIP::compute(int eflag, int vflag) // pre-loop to compute environment coordination f(Z) for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { - j = jlist[neighbor_j]; - j &= NEIGHMASK; + j = jlist[neighbor_j]; + j &= NEIGHMASK; - double dr_ij[3], r_ij; + double dr_ij[3], r_ij; - dr_ij[0] = xtmp - x[j][0]; - dr_ij[1] = ytmp - x[j][1]; - dr_ij[2] = ztmp - x[j][2]; - r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + dr_ij[0] = xtmp - x[j][0]; + dr_ij[1] = ytmp - x[j][1]; + dr_ij[2] = ztmp - x[j][2]; + r_ij = dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2]; - jtype = map[type[j]]; - ijparam = elem3param[itype][jtype][jtype]; - if (r_ij > params[ijparam].cutsq) continue; + jtype = map[type[j]]; + ijparam = elem3param[itype][jtype][jtype]; + if (r_ij > params[ijparam].cutsq) continue; - r_ij = sqrt(r_ij); + r_ij = sqrt(r_ij); - invR_ij = 1.0 / r_ij; - preInvR_ij[neighbor_j] = invR_ij; + invR_ij = 1.0 / r_ij; + preInvR_ij[neighbor_j] = invR_ij; - invRMinusCutoffA = 1.0 / (r_ij - cutoffA); - sigmaInvRMinusCutoffA = sigma * invRMinusCutoffA; - gammInvRMinusCutoffA = gamm * invRMinusCutoffA; + invRMinusCutoffA = 1.0 / (r_ij - cutoffA); + sigmaInvRMinusCutoffA = sigma * invRMinusCutoffA; + gammInvRMinusCutoffA = gamm * invRMinusCutoffA; - interpolDeltaX = r_ij - GRIDSTART; - interpolTMP = (interpolDeltaX * GRIDDENSITY); - interpolIDX = (int) interpolTMP; + interpolDeltaX = r_ij - GRIDSTART; + interpolTMP = (interpolDeltaX * GRIDDENSITY); + interpolIDX = (int) interpolTMP; - interpolY1 = exp3B[interpolIDX]; - interpolY2 = exp3B[interpolIDX+1]; - exp3B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = exp3B[interpolIDX]; + interpolY2 = exp3B[interpolIDX + 1]; + exp3B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - exp3BDerived_ij = - exp3B_ij * gammInvRMinusCutoffA * invRMinusCutoffA; + exp3BDerived_ij = -exp3B_ij * gammInvRMinusCutoffA * invRMinusCutoffA; - preExp3B_ij[neighbor_j] = exp3B_ij; - preExp3BDerived_ij[neighbor_j] = exp3BDerived_ij; + preExp3B_ij[neighbor_j] = exp3B_ij; + preExp3BDerived_ij[neighbor_j] = exp3BDerived_ij; - interpolY1 = exp2B[interpolIDX]; - interpolY2 = exp2B[interpolIDX+1]; - exp2B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = exp2B[interpolIDX]; + interpolY2 = exp2B[interpolIDX + 1]; + exp2B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - exp2BDerived_ij = - exp2B_ij * sigmaInvRMinusCutoffA * invRMinusCutoffA; + exp2BDerived_ij = -exp2B_ij * sigmaInvRMinusCutoffA * invRMinusCutoffA; - preExp2B_ij[neighbor_j] = exp2B_ij; - preExp2BDerived_ij[neighbor_j] = exp2BDerived_ij; + preExp2B_ij[neighbor_j] = exp2B_ij; + preExp2BDerived_ij[neighbor_j] = exp2BDerived_ij; - interpolY1 = pow2B[interpolIDX]; - interpolY2 = pow2B[interpolIDX+1]; - pow2B_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = pow2B[interpolIDX]; + interpolY2 = pow2B[interpolIDX + 1]; + pow2B_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - prePow2B_ij[neighbor_j] = pow2B_ij; + prePow2B_ij[neighbor_j] = pow2B_ij; - // zeta and its derivative + // zeta and its derivative - if (r_ij < cutoffC) zeta_i += 1.0; - else { - interpolY1 = cutoffFunction[interpolIDX]; - interpolY2 = cutoffFunction[interpolIDX+1]; - cutoffFunction_ij = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + if (r_ij < cutoffC) + zeta_i += 1.0; + else { + interpolY1 = cutoffFunction[interpolIDX]; + interpolY2 = cutoffFunction[interpolIDX + 1]; + cutoffFunction_ij = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - zeta_i += cutoffFunction_ij; + zeta_i += cutoffFunction_ij; - interpolY1 = cutoffFunctionDerived[interpolIDX]; - interpolY2 = cutoffFunctionDerived[interpolIDX+1]; - zeta_iDerived = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY1 = cutoffFunctionDerived[interpolIDX]; + interpolY2 = cutoffFunctionDerived[interpolIDX + 1]; + zeta_iDerived = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); - zeta_iDerivedInvR_ij = zeta_iDerived * invR_ij; + zeta_iDerivedInvR_ij = zeta_iDerived * invR_ij; - preForceCoord_counter=numForceCoordPairs*5; - preForceCoord[preForceCoord_counter+0]=zeta_iDerivedInvR_ij; - preForceCoord[preForceCoord_counter+1]=dr_ij[0]; - preForceCoord[preForceCoord_counter+2]=dr_ij[1]; - preForceCoord[preForceCoord_counter+3]=dr_ij[2]; - preForceCoord[preForceCoord_counter+4]=j; - numForceCoordPairs++; - } + preForceCoord_counter = numForceCoordPairs * 5; + preForceCoord[preForceCoord_counter + 0] = zeta_iDerivedInvR_ij; + preForceCoord[preForceCoord_counter + 1] = dr_ij[0]; + preForceCoord[preForceCoord_counter + 2] = dr_ij[1]; + preForceCoord[preForceCoord_counter + 3] = dr_ij[2]; + preForceCoord[preForceCoord_counter + 4] = j; + numForceCoordPairs++; + } } // quantities depending on zeta_i @@ -269,24 +263,20 @@ void PairEDIP::compute(int eflag, int vflag) interpolIDX = (int) interpolTMP; interpolY1 = expMinusBetaZeta_iZeta_iGrid[interpolIDX]; - interpolY2 = expMinusBetaZeta_iZeta_iGrid[interpolIDX+1]; - expMinusBetaZeta_iZeta_i = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = expMinusBetaZeta_iZeta_iGrid[interpolIDX + 1]; + expMinusBetaZeta_iZeta_i = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = qFunctionGrid[interpolIDX]; - interpolY2 = qFunctionGrid[interpolIDX+1]; - qFunction = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = qFunctionGrid[interpolIDX + 1]; + qFunction = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = tauFunctionGrid[interpolIDX]; - interpolY2 = tauFunctionGrid[interpolIDX+1]; - tauFunction = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = tauFunctionGrid[interpolIDX + 1]; + tauFunction = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); interpolY1 = tauFunctionDerivedGrid[interpolIDX]; - interpolY2 = tauFunctionDerivedGrid[interpolIDX+1]; - tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) * - (interpolTMP-interpolIDX); + interpolY2 = tauFunctionDerivedGrid[interpolIDX + 1]; + tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); forceModCoord_factor = 2.0 * beta * zeta_i * expMinusBetaZeta_iZeta_i; @@ -303,7 +293,7 @@ void PairEDIP::compute(int eflag, int vflag) dr_ij[0] = x[j][0] - xtmp; dr_ij[1] = x[j][1] - ytmp; dr_ij[2] = x[j][2] - ztmp; - r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + r_ij = dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2]; jtype = map[type[j]]; ijparam = elem3param[itype][jtype][jtype]; @@ -318,13 +308,12 @@ void PairEDIP::compute(int eflag, int vflag) exp2B_ij = preExp2B_ij[neighbor_j]; - pow2BDerived_ij = - rho * invR_ij * pow2B_ij; + pow2BDerived_ij = -rho * invR_ij * pow2B_ij; - forceModCoord += (forceModCoord_factor*exp2B_ij); + forceModCoord += (forceModCoord_factor * exp2B_ij); exp2BDerived_ij = preExp2BDerived_ij[neighbor_j]; - forceMod2B = exp2BDerived_ij * potential2B_factor + - exp2B_ij * pow2BDerived_ij; + forceMod2B = exp2BDerived_ij * potential2B_factor + exp2B_ij * pow2BDerived_ij; directorCos_ij_x = invR_ij * dr_ij[0]; directorCos_ij_y = invR_ij * dr_ij[1]; @@ -349,136 +338,123 @@ void PairEDIP::compute(int eflag, int vflag) evdwl = (exp2B_ij * potential2B_factor); - if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, - -forceMod2B*invR_ij, dr_ij[0], dr_ij[1], dr_ij[2]); + if (evflag) + ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, -forceMod2B * invR_ij, dr_ij[0], dr_ij[1], + dr_ij[2]); // three-body Forces for (int neighbor_k = neighbor_j + 1; neighbor_k < jnum; neighbor_k++) { - double dr_ik[3], r_ik, f_ik[3]; + double dr_ik[3], r_ik, f_ik[3]; - k = jlist[neighbor_k]; - k &= NEIGHMASK; - ktype = map[type[k]]; - ikparam = elem3param[itype][ktype][ktype]; + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem3param[itype][ktype][ktype]; - dr_ik[0] = x[k][0] - xtmp; - dr_ik[1] = x[k][1] - ytmp; - dr_ik[2] = x[k][2] - ztmp; - r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + dr_ik[0] = x[k][0] - xtmp; + dr_ik[1] = x[k][1] - ytmp; + dr_ik[2] = x[k][2] - ztmp; + r_ik = dr_ik[0] * dr_ik[0] + dr_ik[1] * dr_ik[1] + dr_ik[2] * dr_ik[2]; - if (r_ik > params[ikparam].cutsq) continue; + if (r_ik > params[ikparam].cutsq) continue; - r_ik = sqrt(r_ik); + r_ik = sqrt(r_ik); - invR_ik = preInvR_ij[neighbor_k]; + invR_ik = preInvR_ij[neighbor_k]; - directorCos_ik_x = invR_ik * dr_ik[0]; - directorCos_ik_y = invR_ik * dr_ik[1]; - directorCos_ik_z = invR_ik * dr_ik[2]; + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; - cosTeta = directorCos_ij_x * directorCos_ik_x + - directorCos_ij_y * directorCos_ik_y + + cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + directorCos_ij_z * directorCos_ik_z; - cosTetaDiff = cosTeta + tauFunction; - cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; - qFunctionCosTetaDiffCosTetaDiff = cosTetaDiffCosTetaDiff * qFunction; - expMinusQFunctionCosTetaDiffCosTetaDiff = - exp(-qFunctionCosTetaDiffCosTetaDiff); + cosTetaDiff = cosTeta + tauFunction; + cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; + qFunctionCosTetaDiffCosTetaDiff = cosTetaDiffCosTetaDiff * qFunction; + expMinusQFunctionCosTetaDiffCosTetaDiff = exp(-qFunctionCosTetaDiffCosTetaDiff); - potentia3B_factor = lambda * + potentia3B_factor = lambda * ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + eta * qFunctionCosTetaDiffCosTetaDiff); - exp3B_ik = preExp3B_ij[neighbor_k]; - exp3BDerived_ik = preExp3BDerived_ij[neighbor_k]; + exp3B_ik = preExp3B_ij[neighbor_k]; + exp3BDerived_ik = preExp3BDerived_ij[neighbor_k]; - forceMod3B_factor1_ij = - exp3BDerived_ij * exp3B_ik * - potentia3B_factor; - forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * - qFunction * cosTetaDiff * + forceMod3B_factor1_ij = -exp3BDerived_ij * exp3B_ik * potentia3B_factor; + forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff * (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); - forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; + forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; - f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_x - directorCos_ik_x); - f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_y - directorCos_ik_y); - f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + - forceMod3B_factor2_ij * - (cosTeta * directorCos_ij_z - directorCos_ik_z); + f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); + f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); + f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); - forceMod3B_factor1_ik = - exp3BDerived_ik * exp3B_ij * - potentia3B_factor; - forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; + forceMod3B_factor1_ik = -exp3BDerived_ik * exp3B_ij * potentia3B_factor; + forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; - f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_x - directorCos_ij_x); - f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_y - directorCos_ij_y); - f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + - forceMod3B_factor2_ik * - (cosTeta * directorCos_ik_z - directorCos_ij_z); + f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); + f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); + f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); - forceModCoord += (forceMod3B_factor2 * - (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); + forceModCoord += (forceMod3B_factor2 * (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); - f[j][0] += f_ij[0]; - f[j][1] += f_ij[1]; - f[j][2] += f_ij[2]; + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; - f[k][0] += f_ik[0]; - f[k][1] += f_ik[1]; - f[k][2] += f_ik[2]; + f[k][0] += f_ik[0]; + f[k][1] += f_ik[1]; + f[k][2] += f_ik[2]; - f[i][0] -= f_ij[0] + f_ik[0]; - f[i][1] -= f_ij[1] + f_ik[1]; - f[i][2] -= f_ij[2] + f_ik[2]; + f[i][0] -= f_ij[0] + f_ik[0]; + f[i][1] -= f_ij[1] + f_ik[1]; + f[i][2] -= f_ij[2] + f_ik[2]; - // potential energy + // potential energy - evdwl = (exp3B_ij * exp3B_ik * potentia3B_factor); + evdwl = (exp3B_ij * exp3B_ik * potentia3B_factor); - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + if (evflag) ev_tally3(i, j, k, evdwl, 0.0, f_ij, f_ik, dr_ij, dr_ik); } } // forces due to environment coordination f(Z) for (int idx = 0; idx < numForceCoordPairs; idx++) { - double dr_ij[3],f_ij[3]; + double dr_ij[3], f_ij[3]; - preForceCoord_counter = idx * 5; - zeta_iDerivedInvR_ij=preForceCoord[preForceCoord_counter+0]; - dr_ij[0]=preForceCoord[preForceCoord_counter+1]; - dr_ij[1]=preForceCoord[preForceCoord_counter+2]; - dr_ij[2]=preForceCoord[preForceCoord_counter+3]; - j = static_cast (preForceCoord[preForceCoord_counter+4]); + preForceCoord_counter = idx * 5; + zeta_iDerivedInvR_ij = preForceCoord[preForceCoord_counter + 0]; + dr_ij[0] = preForceCoord[preForceCoord_counter + 1]; + dr_ij[1] = preForceCoord[preForceCoord_counter + 2]; + dr_ij[2] = preForceCoord[preForceCoord_counter + 3]; + j = static_cast(preForceCoord[preForceCoord_counter + 4]); - forceModCoord_ij = forceModCoord * zeta_iDerivedInvR_ij; + forceModCoord_ij = forceModCoord * zeta_iDerivedInvR_ij; - f_ij[0] = forceModCoord_ij * dr_ij[0]; - f_ij[1] = forceModCoord_ij * dr_ij[1]; - f_ij[2] = forceModCoord_ij * dr_ij[2]; + f_ij[0] = forceModCoord_ij * dr_ij[0]; + f_ij[1] = forceModCoord_ij * dr_ij[1]; + f_ij[2] = forceModCoord_ij * dr_ij[2]; - f[i][0] -= f_ij[0]; - f[i][1] -= f_ij[1]; - f[i][2] -= f_ij[2]; + f[i][0] -= f_ij[0]; + f[i][1] -= f_ij[1]; + f[i][2] -= f_ij[2]; - f[j][0] += f_ij[0]; - f[j][1] += f_ij[1]; - f[j][2] += f_ij[2]; + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; - // potential energy - - evdwl = 0.0; - if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, - -forceModCoord_ij, dr_ij[0], dr_ij[1], dr_ij[2]); + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, -forceModCoord_ij, dr_ij[0], dr_ij[1], + dr_ij[2]); } } @@ -500,58 +476,51 @@ void PairEDIP::allocateGrids(void) double maxArgumentTauFunctionGrid; double maxArgumentQFunctionGrid; double maxArgumentExpMinusBetaZeta_iZeta_i; - double const leftLimitToZero = -DBL_MIN * 1000.0; + double const leftLimitToZero = -std::numeric_limits::min() * 1000.0; deallocateGrids(); // tauFunctionGrid maxArgumentTauFunctionGrid = leadDimInteractionList; - numGridPointsTauFunctionGrid = (int) - ((maxArgumentTauFunctionGrid) * GRIDDENSITY) + 2; + numGridPointsTauFunctionGrid = (int) ((maxArgumentTauFunctionGrid) *GRIDDENSITY) + 2; - memory->create(tauFunctionGrid,numGridPointsTauFunctionGrid, - "edip:tauFunctionGrid"); - memory->create(tauFunctionDerivedGrid,numGridPointsTauFunctionGrid, + memory->create(tauFunctionGrid, numGridPointsTauFunctionGrid, "edip:tauFunctionGrid"); + memory->create(tauFunctionDerivedGrid, numGridPointsTauFunctionGrid, "edip:tauFunctionDerivedGrid"); // expMinusBetaZeta_iZeta_iGrid maxArgumentExpMinusBetaZeta_iZeta_i = leadDimInteractionList; - numGridPointsExpMinusBetaZeta_iZeta_i = (int) - ((maxArgumentExpMinusBetaZeta_iZeta_i) * GRIDDENSITY) + 2; - memory->create(expMinusBetaZeta_iZeta_iGrid, - numGridPointsExpMinusBetaZeta_iZeta_i, + numGridPointsExpMinusBetaZeta_iZeta_i = + (int) ((maxArgumentExpMinusBetaZeta_iZeta_i) *GRIDDENSITY) + 2; + memory->create(expMinusBetaZeta_iZeta_iGrid, numGridPointsExpMinusBetaZeta_iZeta_i, "edip:expMinusBetaZeta_iZeta_iGrid"); // qFunctionGrid maxArgumentQFunctionGrid = leadDimInteractionList; - numGridPointsQFunctionGrid = (int) - ((maxArgumentQFunctionGrid) * GRIDDENSITY) + 2; - memory->create(qFunctionGrid,numGridPointsQFunctionGrid,"edip:qFunctionGrid"); + numGridPointsQFunctionGrid = (int) ((maxArgumentQFunctionGrid) *GRIDDENSITY) + 2; + memory->create(qFunctionGrid, numGridPointsQFunctionGrid, "edip:qFunctionGrid"); // cutoffFunction numGridPointsOneCutoffFunction = (int) ((cutoffC - GRIDSTART) * GRIDDENSITY); - numGridPointsNotOneCutoffFunction = (int) ((cutoffA-cutoffC) * GRIDDENSITY); - numGridPointsCutoffFunction = numGridPointsOneCutoffFunction + - numGridPointsNotOneCutoffFunction+2; + numGridPointsNotOneCutoffFunction = (int) ((cutoffA - cutoffC) * GRIDDENSITY); + numGridPointsCutoffFunction = + numGridPointsOneCutoffFunction + numGridPointsNotOneCutoffFunction + 2; - memory->create(cutoffFunction,numGridPointsCutoffFunction, - "edip:cutoffFunction"); - memory->create(cutoffFunctionDerived,numGridPointsCutoffFunction, - "edip:cutoffFunctionDerived"); + memory->create(cutoffFunction, numGridPointsCutoffFunction, "edip:cutoffFunction"); + memory->create(cutoffFunctionDerived, numGridPointsCutoffFunction, "edip:cutoffFunctionDerived"); // pow2B - numGridPointsR = (int) - ((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY); + numGridPointsR = (int) ((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY); numGridPointsRTotal = numGridPointsR + 2; - memory->create(pow2B,numGridPointsRTotal,"edip:pow2B"); - memory->create(exp2B,numGridPointsRTotal,"edip:exp2B"); - memory->create(exp3B,numGridPointsRTotal,"edip:exp3B"); + memory->create(pow2B, numGridPointsRTotal, "edip:pow2B"); + memory->create(exp2B, numGridPointsRTotal, "edip:exp2B"); + memory->create(exp3B, numGridPointsRTotal, "edip:exp3B"); } /* ---------------------------------------------------------------------- @@ -563,15 +532,13 @@ 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, - "edip:preExp3BDerived_ij"); - memory->create(preExp2B_ij,nthreads*leadDimInteractionList,"edip:preExp2B_ij"); - memory->create(preExp2BDerived_ij,nthreads*leadDimInteractionList, - "edip:preExp2BDerived_ij"); - memory->create(prePow2B_ij,nthreads*leadDimInteractionList,"edip:prePow2B_ij"); - memory->create(preForceCoord,5*nthreads*leadDimInteractionList,"edip:preForceCoord"); + memory->create(preInvR_ij, nthreads * leadDimInteractionList, "edip:preInvR_ij"); + memory->create(preExp3B_ij, nthreads * leadDimInteractionList, "edip:preExp3B_ij"); + memory->create(preExp3BDerived_ij, nthreads * leadDimInteractionList, "edip:preExp3BDerived_ij"); + memory->create(preExp2B_ij, nthreads * leadDimInteractionList, "edip:preExp2B_ij"); + memory->create(preExp2BDerived_ij, nthreads * leadDimInteractionList, "edip:preExp2BDerived_ij"); + memory->create(prePow2B_ij, nthreads * leadDimInteractionList, "edip:prePow2B_ij"); + memory->create(preForceCoord, 5 * nthreads * leadDimInteractionList, "edip:preForceCoord"); } /* ---------------------------------------------------------------------- @@ -613,19 +580,19 @@ void PairEDIP::allocate() allocated = 1; int n = atom->ntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); - map = new int[n+1]; + map = new int[n + 1]; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ -void PairEDIP::settings(int narg, char **/*arg*/) +void PairEDIP::settings(int narg, char ** /*arg*/) { - if (narg != 0) error->all(FLERR,"Illegal pair_style command"); + if (narg != 0) error->all(FLERR, "Illegal pair_style command"); } /* ---------------------------------------------------------------------- */ @@ -652,105 +619,97 @@ void PairEDIP::initGrids(void) double deltaArgumentQFunctionGrid; double deltaArgumentTauFunctionGrid; double deltaArgumentExpMinusBetaZeta_iZeta_i; - double const leftLimitToZero = -DBL_MIN * 1000.0; + double const leftLimitToZero = -std::numeric_limits::min() * 1000.0; // tauFunctionGrid maxArgumentTauFunctionGrid = leadDimInteractionList; - numGridPointsTauFunctionGrid = (int) - ((maxArgumentTauFunctionGrid) * GRIDDENSITY) + 2; + numGridPointsTauFunctionGrid = (int) ((maxArgumentTauFunctionGrid) *GRIDDENSITY) + 2; r = 0.0; deltaArgumentTauFunctionGrid = 1.0 / GRIDDENSITY; for (l = 0; l < numGridPointsTauFunctionGrid; l++) { - tauFunctionGrid[l] = u1 + u2 * u3 * exp(-u4 * r) - - u2 * exp(-2.0 * u4 * r); - tauFunctionDerivedGrid[l] = - u2 * u3 * u4 * exp(-u4 * r) + - 2.0 * u2 * u4 * exp(-2.0 * u4 * r); - r += deltaArgumentTauFunctionGrid; + tauFunctionGrid[l] = u1 + u2 * u3 * exp(-u4 * r) - u2 * exp(-2.0 * u4 * r); + tauFunctionDerivedGrid[l] = -u2 * u3 * u4 * exp(-u4 * r) + 2.0 * u2 * u4 * exp(-2.0 * u4 * r); + r += deltaArgumentTauFunctionGrid; } // expMinusBetaZeta_iZeta_iGrid maxArgumentExpMinusBetaZeta_iZeta_i = leadDimInteractionList; - numGridPointsExpMinusBetaZeta_iZeta_i = (int) - ((maxArgumentExpMinusBetaZeta_iZeta_i) * GRIDDENSITY) + 2; + numGridPointsExpMinusBetaZeta_iZeta_i = + (int) ((maxArgumentExpMinusBetaZeta_iZeta_i) *GRIDDENSITY) + 2; r = 0.0; deltaArgumentExpMinusBetaZeta_iZeta_i = 1.0 / GRIDDENSITY; for (l = 0; l < numGridPointsExpMinusBetaZeta_iZeta_i; l++) { - expMinusBetaZeta_iZeta_iGrid[l] = exp(-beta * r * r); - r += deltaArgumentExpMinusBetaZeta_iZeta_i; + expMinusBetaZeta_iZeta_iGrid[l] = exp(-beta * r * r); + r += deltaArgumentExpMinusBetaZeta_iZeta_i; } // qFunctionGrid maxArgumentQFunctionGrid = leadDimInteractionList; - numGridPointsQFunctionGrid = - (int) ((maxArgumentQFunctionGrid) * GRIDDENSITY) + 2; + numGridPointsQFunctionGrid = (int) ((maxArgumentQFunctionGrid) *GRIDDENSITY) + 2; r = 0.0; deltaArgumentQFunctionGrid = 1.0 / GRIDDENSITY; for (l = 0; l < numGridPointsQFunctionGrid; l++) { - qFunctionGrid[l] = Q0 * exp(-mu * r); - r += deltaArgumentQFunctionGrid; + qFunctionGrid[l] = Q0 * exp(-mu * r); + r += deltaArgumentQFunctionGrid; } // cutoffFunction - numGridPointsOneCutoffFunction = - (int) ((cutoffC - GRIDSTART) * GRIDDENSITY); - numGridPointsNotOneCutoffFunction = - (int) ((cutoffA-cutoffC) * GRIDDENSITY); + numGridPointsOneCutoffFunction = (int) ((cutoffC - GRIDSTART) * GRIDDENSITY); + numGridPointsNotOneCutoffFunction = (int) ((cutoffA - cutoffC) * GRIDDENSITY); numGridPointsCutoffFunction = - numGridPointsOneCutoffFunction+numGridPointsNotOneCutoffFunction+2; + numGridPointsOneCutoffFunction + numGridPointsNotOneCutoffFunction + 2; r = GRIDSTART; deltaArgumentCutoffFunction = 1.0 / GRIDDENSITY; for (l = 0; l < numGridPointsOneCutoffFunction; l++) { - cutoffFunction[l] = 1.0; - cutoffFunctionDerived[l] = 0.0; - r += deltaArgumentCutoffFunction; + cutoffFunction[l] = 1.0; + cutoffFunctionDerived[l] = 0.0; + r += deltaArgumentCutoffFunction; } - for (l = numGridPointsOneCutoffFunction; - l < numGridPointsCutoffFunction; l++) { - temp = (cutoffA - cutoffC)/(r - cutoffC); - temp3 = temp * temp * temp; - temp4 = temp3 * temp; - cutoffFunction[l] = exp(alpha/(1.0-temp3)); - cutoffFunctionDerived[l] = (-3*alpha/(cutoffA-cutoffC)) * - (temp4/((1-temp3)*(1-temp3)))*exp(alpha/(1.0-temp3)); - r += deltaArgumentCutoffFunction; + for (l = numGridPointsOneCutoffFunction; l < numGridPointsCutoffFunction; l++) { + temp = (cutoffA - cutoffC) / (r - cutoffC); + temp3 = temp * temp * temp; + temp4 = temp3 * temp; + cutoffFunction[l] = exp(alpha / (1.0 - temp3)); + cutoffFunctionDerived[l] = (-3 * alpha / (cutoffA - cutoffC)) * + (temp4 / ((1 - temp3) * (1 - temp3))) * exp(alpha / (1.0 - temp3)); + r += deltaArgumentCutoffFunction; } // pow2B - numGridPointsR = (int) - ((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY); + numGridPointsR = (int) ((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY); r = GRIDSTART; deltaArgumentR = 1.0 / GRIDDENSITY; for (l = 0; l < numGridPointsR; l++) { - pow2B[l] = pow((B/r),rho); - exp2B[l] = A * exp(sigma/(r-cutoffA)); - exp3B[l] = exp(gamm/(r-cutoffA)); - r += deltaArgumentR; + pow2B[l] = pow((B / r), rho); + exp2B[l] = A * exp(sigma / (r - cutoffA)); + exp3B[l] = exp(gamm / (r - cutoffA)); + r += deltaArgumentR; } - pow2B[numGridPointsR] = pow((B/r),rho); - exp2B[numGridPointsR]=0; - exp3B[numGridPointsR]=0; + pow2B[numGridPointsR] = pow((B / r), rho); + exp2B[numGridPointsR] = 0; + exp3B[numGridPointsR] = 0; r += deltaArgumentR; - pow2B[numGridPointsR+1] = pow((B/r),rho); - exp2B[numGridPointsR+1]=0; - exp3B[numGridPointsR+1]=0; + pow2B[numGridPointsR + 1] = pow((B / r), rho); + exp2B[numGridPointsR + 1] = 0; + exp3B[numGridPointsR + 1] = 0; } /* ---------------------------------------------------------------------- @@ -761,9 +720,8 @@ void PairEDIP::coeff(int narg, char **arg) { if (!allocated) allocate(); - map_element2type(narg-3,arg+3); - if (nelements != 1) - error->all(FLERR,"Pair style edip only supports single element potentials"); + map_element2type(narg - 3, arg + 3); + if (nelements != 1) error->all(FLERR, "Pair style edip only supports single element potentials"); // read potential file and initialize potential parameters @@ -783,12 +741,11 @@ void PairEDIP::coeff(int narg, char **arg) void PairEDIP::init_style() { - if (force->newton_pair == 0) - error->all(FLERR,"Pair style edip requires newton pair on"); + if (force->newton_pair == 0) error->all(FLERR, "Pair style edip requires newton pair on"); // need a full neighbor list - int irequest = neighbor->request(this,instance_me); + int irequest = neighbor->request(this, instance_me); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; } @@ -799,7 +756,7 @@ void PairEDIP::init_style() double PairEDIP::init_one(int i, int j) { - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); return cutmax; } @@ -809,7 +766,7 @@ double PairEDIP::init_one(int i, int j) void PairEDIP::read_file(char *file) { int params_per_line = 20; - char **words = new char*[params_per_line+1]; + char **words = new char *[params_per_line + 1]; memory->sfree(params); params = nullptr; @@ -819,11 +776,11 @@ void PairEDIP::read_file(char *file) FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file,lmp,nullptr); + fp = utils::open_potential(file, lmp, nullptr); if (fp == nullptr) { char str[128]; - snprintf(str,128,"Cannot open EDIP potential file %s",file); - error->one(FLERR,str); + snprintf(str, 128, "Cannot open EDIP potential file %s", file); + error->one(FLERR, str); } } @@ -831,26 +788,27 @@ void PairEDIP::read_file(char *file) // one set of params can span multiple lines // store params if all 3 element tags are in element list - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; + int n, nwords, ielement, jelement, kelement; + char line[MAXLINE], *ptr; int eof = 0; while (1) { if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); + ptr = fgets(line, MAXLINE, fp); if (ptr == nullptr) { eof = 1; fclose(fp); - } else n = strlen(line) + 1; + } else + n = strlen(line) + 1; } - MPI_Bcast(&eof,1,MPI_INT,0,world); + MPI_Bcast(&eof, 1, MPI_INT, 0, world); if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + MPI_Bcast(&n, 1, MPI_INT, 0, world); + MPI_Bcast(line, n, MPI_CHAR, 0, world); // strip comment, skip line if blank - if ((ptr = strchr(line,'#'))) *ptr = '\0'; + if ((ptr = strchr(line, '#'))) *ptr = '\0'; nwords = utils::count_words(line); if (nwords == 0) continue; @@ -859,54 +817,53 @@ void PairEDIP::read_file(char *file) while (nwords < params_per_line) { n = strlen(line); if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); + ptr = fgets(&line[n], MAXLINE - n, fp); if (ptr == nullptr) { eof = 1; fclose(fp); - } else n = strlen(line) + 1; + } else + n = strlen(line) + 1; } - MPI_Bcast(&eof,1,MPI_INT,0,world); + MPI_Bcast(&eof, 1, MPI_INT, 0, world); if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - if ((ptr = strchr(line,'#'))) *ptr = '\0'; + MPI_Bcast(&n, 1, MPI_INT, 0, world); + MPI_Bcast(line, n, MPI_CHAR, 0, world); + if ((ptr = strchr(line, '#'))) *ptr = '\0'; nwords = utils::count_words(line); } - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in EDIP potential file"); + if (nwords != params_per_line) error->all(FLERR, "Incorrect format in EDIP potential file"); // words = ptrs to all words in line nwords = 0; - words[nwords++] = strtok(line," \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue; + words[nwords++] = strtok(line, " \t\n\r\f"); + while ((words[nwords++] = strtok(nullptr, " \t\n\r\f"))) continue; // ielement,jelement,kelement = 1st args // if all 3 args are in element list, then parse this line // else skip to next entry in file for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0],elements[ielement]) == 0) break; + if (strcmp(words[0], elements[ielement]) == 0) break; if (ielement == nelements) continue; for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1],elements[jelement]) == 0) break; + if (strcmp(words[1], elements[jelement]) == 0) break; if (jelement == nelements) continue; for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2],elements[kelement]) == 0) break; + if (strcmp(words[2], elements[kelement]) == 0) break; if (kelement == nelements) continue; // load up parameter settings and error check their values if (nparams == maxparam) { maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); + params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params"); // make certain all addional allocated storage is initialized // to avoid false positives when checking with valgrind - memset(params + nparams, 0, DELTA*sizeof(Param)); + memset(params + nparams, 0, DELTA * sizeof(Param)); } params[nparams].ielement = ielement; @@ -930,25 +887,24 @@ void PairEDIP::read_file(char *file) params[nparams].u3 = atof(words[18]); params[nparams].u4 = atof(words[19]); - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || - params[nparams].cutoffA < 0.0 || params[nparams].cutoffC < 0.0 || - params[nparams].alpha < 0.0 || params[nparams].beta < 0.0 || - params[nparams].eta < 0.0 || params[nparams].gamm < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || - params[nparams].rho < 0.0 || params[nparams].sigma < 0.0) - error->all(FLERR,"Illegal EDIP parameter"); + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamm < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); nparams++; } - delete [] words; + delete[] words; } /* ---------------------------------------------------------------------- */ void PairEDIP::setup_params() { - int i,j,k,m,n; + int i, j, k, m, n; double rtmp; // set elem3param for all triplet combinations @@ -956,28 +912,25 @@ void PairEDIP::setup_params() // do not allow for ACB in place of ABC memory->destroy(elem3param); - memory->create(elem3param,nelements,nelements,nelements,"pair:elem3param"); + memory->create(elem3param, nelements, nelements, nelements, "pair:elem3param"); for (i = 0; i < nelements; i++) for (j = 0; j < nelements; j++) for (k = 0; k < nelements; k++) { n = -1; for (m = 0; m < nparams; m++) { - if (i == params[m].ielement && j == params[m].jelement && - k == params[m].kelement) { - if (n >= 0) error->all(FLERR,"Potential file has duplicate entry"); + if (i == params[m].ielement && j == params[m].jelement && k == params[m].kelement) { + if (n >= 0) error->all(FLERR, "Potential file has duplicate entry"); n = m; } } - if (n < 0) error->all(FLERR,"Potential file is missing an entry"); + if (n < 0) error->all(FLERR, "Potential file is missing an entry"); elem3param[i][j][k] = n; } // set cutoff square - for (m = 0; m < nparams; m++) { - params[m].cutsq = params[m].cutoffA*params[m].cutoffA; - } + for (m = 0; m < nparams; m++) { params[m].cutsq = params[m].cutoffA * params[m].cutoffA; } // set cutmax to max of all params diff --git a/src/USER-MISC/pair_extep.cpp b/src/USER-MISC/pair_extep.cpp index 554aafe8ee..efb5b0e7fb 100644 --- a/src/USER-MISC/pair_extep.cpp +++ b/src/USER-MISC/pair_extep.cpp @@ -354,6 +354,7 @@ void PairExTeP::compute(int eflag, int vflag) f[k][0] -= fc_prefac_ik * delr2[0]; f[k][1] -= fc_prefac_ik * delr2[1]; f[k][2] -= fc_prefac_ik * delr2[2]; + if (vflag_either) v_tally2(i,k,-fc_prefac_ik,delr2); if (itype != ktype) { fc_prefac_ik = dFc_dNdij * fc_prefac_ik_0; f[i][0] += fc_prefac_ik * delr2[0]; @@ -362,6 +363,7 @@ void PairExTeP::compute(int eflag, int vflag) f[k][0] -= fc_prefac_ik * delr2[0]; f[k][1] -= fc_prefac_ik * delr2[1]; f[k][2] -= fc_prefac_ik * delr2[2]; + if (vflag_either) v_tally2(i,k,-fc_prefac_ik,delr2); } /* END F_IJ (2) */ @@ -410,7 +412,7 @@ void PairExTeP::compute(int eflag, int vflag) f[k][1] += fk[1]; f[k][2] += fk[2]; - if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2); } } } diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp index 016213e505..e64e5bce02 100644 --- a/src/USER-MISC/pair_tersoff_table.cpp +++ b/src/USER-MISC/pair_tersoff_table.cpp @@ -145,13 +145,10 @@ void PairTersoffTable::compute(int eflag, int vflag) jlist = firstneigh[i]; jnum = numneigh[i]; - if (jnum > leadingDimensionInteractionList) { - char errmsg[256]; - sprintf(errmsg,"Too many neighbors for interaction list: %d vs %d.\n" - "Check your system or increase 'leadingDimensionInteractionList'", - jnum, leadingDimensionInteractionList); - error->one(FLERR,errmsg); - } + if (jnum > leadingDimensionInteractionList) + error->one(FLERR,"Too many neighbors for interaction list: {} vs {}.\n" + "Check your system or increase 'leadingDimensionInteractionList'", + jnum, leadingDimensionInteractionList); // Pre-calculate gteta and cutoff function for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { @@ -434,10 +431,7 @@ void PairTersoffTable::compute(int eflag, int vflag) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - // potential energy - evdwl = 0.0; - - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + if (vflag_either) v_tally3(i,j,k,f_ij,f_ik,dr_ij,dr_ik); } // second loop over neighbors of atom i except j, forces and virial only - part 2/2 @@ -499,11 +493,7 @@ void PairTersoffTable::compute(int eflag, int vflag) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - // potential energy - evdwl = 0.0; - - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); - + if (vflag_either) v_tally3(i,j,k,f_ij,f_ik,dr_ij,dr_ik); } } // loop on J f[i][0] += fxtmp; @@ -755,7 +745,7 @@ void PairTersoffTable::coeff(int narg, char **arg) void PairTersoffTable::init_style() { if (force->newton_pair == 0) - error->all(FLERR,"Pair style Tersoff requires newton pair on"); + error->all(FLERR,"Pair style tersoff/table requires newton pair on"); // need a full neighbor list @@ -888,9 +878,8 @@ void PairTersoffTable::read_file(char *file) MPI_Bcast(&nparams, 1, MPI_INT, 0, world); MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - if (comm->me != 0) { + if (comm->me != 0) params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); - } MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } diff --git a/src/pair.cpp b/src/pair.cpp index 5448a42e1c..5e8c666743 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1363,22 +1363,40 @@ void Pair::ev_tally4(int i, int j, int k, int m, double evdwl, } } - if (vflag_atom) { - v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); - v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); - v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); - v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); - v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); - v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + if (vflag_either) { + v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; - vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; - vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; - vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; - vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; - vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; - vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] *= 0.25; + v[1] *= 0.25; + v[2] *= 0.25; + v[3] *= 0.25; + v[4] *= 0.25; + v[5] *= 0.25; + + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; + vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + } } } @@ -1499,17 +1517,34 @@ void Pair::v_tally2(int i, int j, double fpair, double *drij) { double v[6]; - v[0] = 0.5 * drij[0]*drij[0]*fpair; - v[1] = 0.5 * drij[1]*drij[1]*fpair; - v[2] = 0.5 * drij[2]*drij[2]*fpair; - v[3] = 0.5 * drij[0]*drij[1]*fpair; - v[4] = 0.5 * drij[0]*drij[2]*fpair; - v[5] = 0.5 * drij[1]*drij[2]*fpair; + v[0] = drij[0]*drij[0]*fpair; + v[1] = drij[1]*drij[1]*fpair; + v[2] = drij[2]*drij[2]*fpair; + v[3] = drij[0]*drij[1]*fpair; + v[4] = drij[0]*drij[2]*fpair; + v[5] = drij[1]*drij[2]*fpair; - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; - vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; - vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] *= 0.5; + v[1] *= 0.5; + v[2] *= 0.5; + v[3] *= 0.5; + v[4] *= 0.5; + v[5] *= 0.5; + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + } } /* ---------------------------------------------------------------------- @@ -1517,24 +1552,40 @@ void Pair::v_tally2(int i, int j, double fpair, double *drij) called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ -void Pair::v_tally3(int i, int j, int k, - double *fi, double *fj, double *drik, double *drjk) +void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, double *drjk) { double v[6]; - v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); - v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); - v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); - v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); - v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); - v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); + v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]); - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; - vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; - vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; - vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; - vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] *= THIRD; + v[1] *= THIRD; + v[2] *= THIRD; + v[3] *= THIRD; + v[4] *= THIRD; + v[5] *= THIRD; + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + } } /* ---------------------------------------------------------------------- @@ -1548,21 +1599,38 @@ void Pair::v_tally4(int i, int j, int k, int m, { double v[6]; - v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); - v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); - v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); - v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); - v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); - v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; - vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; - vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; - vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; - vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; - vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; - vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] *= 0.25; + v[1] *= 0.25; + v[2] *= 0.25; + v[3] *= 0.25; + v[4] *= 0.25; + v[5] *= 0.25; + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; + vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + } } /* ---------------------------------------------------------------------- diff --git a/src/utils.h b/src/utils.h index 58b6df7725..276ae0a315 100644 --- a/src/utils.h +++ b/src/utils.h @@ -33,19 +33,19 @@ class LAMMPS; namespace utils { - /** Match text against a simplified regex pattern - * - * \param text the text to be matched against the pattern - * \param pattern the search pattern, which may contain regexp markers - * \return true if the pattern matches, false if not */ + /*! Match text against a simplified regex pattern + * + * \param text the text to be matched against the pattern + * \param pattern the search pattern, which may contain regexp markers + * \return true if the pattern matches, false if not */ bool strmatch(const std::string &text, const std::string &pattern); - /** Find sub-string that matches a simplified regex pattern - * - * \param text the text to be matched against the pattern - * \param pattern the search pattern, which may contain regexp markers - * \return the string that matches the patters or an empty one */ + /*! Find sub-string that matches a simplified regex pattern + * + * \param text the text to be matched against the pattern + * \param pattern the search pattern, which may contain regexp markers + * \return the string that matches the patters or an empty one */ std::string strfind(const std::string &text, const std::string &pattern); @@ -53,261 +53,261 @@ namespace utils { void fmtargs_logmesg(LAMMPS *lmp, fmt::string_view format, fmt::format_args args); - /** Send formatted message to screen and logfile, if available - * - * This function simplifies the repetitive task of outputting some - * message to both the screen and/or the log file. The template - * wrapper with fmtlib format and argument processing allows - * this function to work similar to ``fmt::print()``. - * - * \param lmp pointer to LAMMPS class instance - * \param format format string of message to be printed - * \param args arguments to format string */ + /*! Send formatted message to screen and logfile, if available + * + * This function simplifies the repetitive task of outputting some + * message to both the screen and/or the log file. The template + * wrapper with fmtlib format and argument processing allows + * this function to work similar to ``fmt::print()``. + * + * \param lmp pointer to LAMMPS class instance + * \param format format string of message to be printed + * \param ... arguments to format string */ template void logmesg(LAMMPS *lmp, const S &format, Args &&...args) { fmtargs_logmesg(lmp, format, fmt::make_args_checked(format, args...)); } - /** \overload - * - * \param lmp pointer to LAMMPS class instance - * \param mesg string with message to be printed */ + /*! \overload + * + * \param lmp pointer to LAMMPS class instance + * \param mesg string with message to be printed */ void logmesg(LAMMPS *lmp, const std::string &mesg); - /** Return a string representing the current system error status - * - * This is a wrapper around calling strerror(errno). - * - * \return error string */ + /*! Return a string representing the current system error status + * + * This is a wrapper around calling strerror(errno). + * + * \return error string */ std::string getsyserror(); - /** Wrapper around fgets() which reads whole lines but truncates the - * data to the buffer size and ensures a newline char at the end. - * - * This function is useful for reading line based text files with - * possible comments that should be parsed later. This applies to - * data files, potential files, atomfile variable files and so on. - * It is used instead of fgets() by utils::read_lines_from_file(). - * - * \param s buffer for storing the result of fgets() - * \param size size of buffer s (max number of bytes returned) - * \param fp file pointer used by fgets() */ + /*! Wrapper around fgets() which reads whole lines but truncates the + * data to the buffer size and ensures a newline char at the end. + * + * This function is useful for reading line based text files with + * possible comments that should be parsed later. This applies to + * data files, potential files, atomfile variable files and so on. + * It is used instead of fgets() by utils::read_lines_from_file(). + * + * \param s buffer for storing the result of fgets() + * \param size size of buffer s (max number of bytes returned) + * \param fp file pointer used by fgets() */ char *fgets_trunc(char *s, int size, FILE *fp); - /** Safe wrapper around fgets() which aborts on errors - * or EOF and prints a suitable error message to help debugging. - * - * Use nullptr as the error parameter to avoid the abort on EOF or error. - * - * \param srcname name of the calling source file (from FLERR macro) - * \param srcline line in the calling source file (from FLERR macro) - * \param s buffer for storing the result of fgets() - * \param size size of buffer s (max number of bytes read by fgets()) - * \param fp file pointer used by fgets() - * \param filename file name associated with fp (may be a null pointer; then LAMMPS will try to detect) - * \param error pointer to Error class instance (for abort) or nullptr */ + /*! Safe wrapper around fgets() which aborts on errors + * or EOF and prints a suitable error message to help debugging. + * + * Use nullptr as the error parameter to avoid the abort on EOF or error. + * + * \param srcname name of the calling source file (from FLERR macro) + * \param srcline line in the calling source file (from FLERR macro) + * \param s buffer for storing the result of fgets() + * \param size size of buffer s (max number of bytes read by fgets()) + * \param fp file pointer used by fgets() + * \param filename file name associated with fp (may be a null pointer; then LAMMPS will try to detect) + * \param error pointer to Error class instance (for abort) or nullptr */ void sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp, const char *filename, Error *error); - /** Safe wrapper around fread() which aborts on errors - * or EOF and prints a suitable error message to help debugging. - * - * Use nullptr as the error parameter to avoid the abort on EOF or error. - * - * \param srcname name of the calling source file (from FLERR macro) - * \param srcline line in the calling source file (from FLERR macro) - * \param s buffer for storing the result of fread() - * \param size size of data elements read by fread() - * \param num number of data elements read by fread() - * \param fp file pointer used by fread() - * \param filename file name associated with fp (may be a null pointer; then LAMMPS will try to detect) - * \param error pointer to Error class instance (for abort) or nullptr */ + /*! Safe wrapper around fread() which aborts on errors + * or EOF and prints a suitable error message to help debugging. + * + * Use nullptr as the error parameter to avoid the abort on EOF or error. + * + * \param srcname name of the calling source file (from FLERR macro) + * \param srcline line in the calling source file (from FLERR macro) + * \param s buffer for storing the result of fread() + * \param size size of data elements read by fread() + * \param num number of data elements read by fread() + * \param fp file pointer used by fread() + * \param filename file name associated with fp (may be a null pointer; then LAMMPS will try to detect) + * \param error pointer to Error class instance (for abort) or nullptr */ void sfread(const char *srcname, int srcline, void *s, size_t size, size_t num, FILE *fp, const char *filename, Error *error); - /** Read N lines of text from file into buffer and broadcast them - * - * This function uses repeated calls to fread() to fill a buffer with - * newline terminated text. If a line does not end in a newline (e.g. - * at the end of a file), it is added. The caller has to allocate an - * nlines by nmax sized buffer for storing the text data. - * Reading is done by MPI rank 0 of the given communicator only, and - * thus only MPI rank 0 needs to provide a valid file pointer. - * - * \param fp file pointer used by fread - * \param nlines number of lines to be read - * \param nmax maximum length of a single line - * \param buffer buffer for storing the data. - * \param me MPI rank of calling process in MPI communicator - * \param comm MPI communicator for broadcast - * \return 1 if the read was short, 0 if read was successful */ + /*! Read N lines of text from file into buffer and broadcast them + * + * This function uses repeated calls to fread() to fill a buffer with + * newline terminated text. If a line does not end in a newline (e.g. + * at the end of a file), it is added. The caller has to allocate an + * nlines by nmax sized buffer for storing the text data. + * Reading is done by MPI rank 0 of the given communicator only, and + * thus only MPI rank 0 needs to provide a valid file pointer. + * + * \param fp file pointer used by fread + * \param nlines number of lines to be read + * \param nmax maximum length of a single line + * \param buffer buffer for storing the data. + * \param me MPI rank of calling process in MPI communicator + * \param comm MPI communicator for broadcast + * \return 1 if the read was short, 0 if read was successful */ int read_lines_from_file(FILE *fp, int nlines, int nmax, char *buffer, int me, MPI_Comm comm); - /** Report if a requested style is in a package or may have a typo - * - * \param style type of style that is to be checked for - * \param name name of style that was not found - * \param lmp pointer to top-level LAMMPS class instance - * \return string usable for error messages */ + /*! Report if a requested style is in a package or may have a typo + * + * \param style type of style that is to be checked for + * \param name name of style that was not found + * \param lmp pointer to top-level LAMMPS class instance + * \return string usable for error messages */ std::string check_packages_for_style(const std::string &style, const std::string &name, LAMMPS *lmp); - /** Convert a string to a floating point number while checking - * if it is a valid floating point or integer number - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param str string to be converted to number - * \param do_abort determines whether to call Error::one() or Error::all() - * \param lmp pointer to top-level LAMMPS class instance - * \return double precision floating point number - */ + /*! Convert a string to a floating point number while checking + * if it is a valid floating point or integer number + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param str string to be converted to number + * \param do_abort determines whether to call Error::one() or Error::all() + * \param lmp pointer to top-level LAMMPS class instance + * \return double precision floating point number */ + double numeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** Convert a string to an integer number while checking - * if it is a valid integer number (regular int) - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param str string to be converted to number - * \param do_abort determines whether to call Error::one() or Error::all() - * \param lmp pointer to top-level LAMMPS class instance - * \return integer number (regular int) */ + /*! Convert a string to an integer number while checking + * if it is a valid integer number (regular int) + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param str string to be converted to number + * \param do_abort determines whether to call Error::one() or Error::all() + * \param lmp pointer to top-level LAMMPS class instance + * \return integer number (regular int) */ int inumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** Convert a string to an integer number while checking - * if it is a valid integer number (bigint) - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param str string to be converted to number - * \param do_abort determines whether to call Error::one() or Error::all() - * \param lmp pointer to top-level LAMMPS class instance - * \return integer number (bigint) */ + /*! Convert a string to an integer number while checking + * if it is a valid integer number (bigint) + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param str string to be converted to number + * \param do_abort determines whether to call Error::one() or Error::all() + * \param lmp pointer to top-level LAMMPS class instance + * \return integer number (bigint) */ bigint bnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** Convert a string to an integer number while checking - * if it is a valid integer number (tagint) - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param str string to be converted to number - * \param do_abort determines whether to call Error::one() or Error::all() - * \param lmp pointer to top-level LAMMPS class instance - * \return integer number (tagint) */ + /*! Convert a string to an integer number while checking + * if it is a valid integer number (tagint) + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param str string to be converted to number + * \param do_abort determines whether to call Error::one() or Error::all() + * \param lmp pointer to top-level LAMMPS class instance + * \return integer number (tagint) */ tagint tnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** Compute index bounds derived from a string with a possible wildcard - * - * This functions processes the string in *str* and set the values of *nlo* - * and *nhi* according to the following five cases: - * - * - a single number, i: nlo = i; nhi = i; - * - a single asterisk, \*: nlo = nmin; nhi = nmax; - * - a single number followed by an asterisk, i\*: nlo = i; nhi = nmax; - * - a single asterisk followed by a number, \*i: nlo = nmin; nhi = i; - * - two numbers with an asterisk in between. i\*j: nlo = i; nhi = j; - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param str string to be processed - * \param nmin smallest possible lower bound - * \param nmax largest allowed upper bound - * \param nlo lower bound - * \param nhi upper bound - * \param error pointer to Error class for out-of-bounds messages */ + /*! Compute index bounds derived from a string with a possible wildcard + * + * This functions processes the string in *str* and set the values of *nlo* + * and *nhi* according to the following five cases: + * + * - a single number, i: nlo = i; nhi = i; + * - a single asterisk, \*: nlo = nmin; nhi = nmax; + * - a single number followed by an asterisk, i\*: nlo = i; nhi = nmax; + * - a single asterisk followed by a number, \*i: nlo = nmin; nhi = i; + * - two numbers with an asterisk in between. i\*j: nlo = i; nhi = j; + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param str string to be processed + * \param nmin smallest possible lower bound + * \param nmax largest allowed upper bound + * \param nlo lower bound + * \param nhi upper bound + * \param error pointer to Error class for out-of-bounds messages */ template void bounds(const char *file, int line, const std::string &str, bigint nmin, bigint nmax, TYPE &nlo, TYPE &nhi, Error *error); - /** Expand list of arguments when containing fix/compute wildcards - * - * This function searches the list of arguments in *arg* for strings - * of the kind c_ID[*] or f_ID[*] referring to computes or fixes. - * Any such strings are replaced by one or more strings with the - * '*' character replaced by the corresponding possible numbers as - * determined from the fix or compute instance. Other strings are - * just copied. If the *mode* parameter is set to 0, expand global - * vectors, but not global arrays; if it is set to 1, expand global - * arrays (by column) but not global vectors. - * - * If any expansion happens, the earg list and all its - * strings are new allocations and must be freed explicitly by the - * caller. Otherwise arg and earg will point to the same address - * and no explicit de-allocation is needed by the caller. - * - * \param file name of source file for error message - * \param line line number in source file for error message - * \param narg number of arguments in current list - * \param arg argument list, possibly containing wildcards - * \param mode select between global vectors(=0) and arrays (=1) - * \param earg new argument list with wildcards expanded - * \param lmp pointer to top-level LAMMPS class instance - * \return number of arguments in expanded list */ + /*! Expand list of arguments when containing fix/compute wildcards + * + * This function searches the list of arguments in *arg* for strings + * of the kind c_ID[*] or f_ID[*] referring to computes or fixes. + * Any such strings are replaced by one or more strings with the + * '*' character replaced by the corresponding possible numbers as + * determined from the fix or compute instance. Other strings are + * just copied. If the *mode* parameter is set to 0, expand global + * vectors, but not global arrays; if it is set to 1, expand global + * arrays (by column) but not global vectors. + * + * If any expansion happens, the earg list and all its + * strings are new allocations and must be freed explicitly by the + * caller. Otherwise arg and earg will point to the same address + * and no explicit de-allocation is needed by the caller. + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param narg number of arguments in current list + * \param arg argument list, possibly containing wildcards + * \param mode select between global vectors(=0) and arrays (=1) + * \param earg new argument list with wildcards expanded + * \param lmp pointer to top-level LAMMPS class instance + * \return number of arguments in expanded list */ int expand_args(const char *file, int line, int narg, char **arg, int mode, char **&earg, LAMMPS *lmp); - /** Make C-style copy of string in new storage - * - * This allocates a storage buffer and copies the C-style or - * C++ style string into it. The buffer is allocated with "new" - * and thus needs to be deallocated with "delete[]". - * - * \param text string that should be copied - * \return new buffer with copy of string */ + /*! Make C-style copy of string in new storage + * + * This allocates a storage buffer and copies the C-style or + * C++ style string into it. The buffer is allocated with "new" + * and thus needs to be deallocated with "delete[]". + * + * \param text string that should be copied + * \return new buffer with copy of string */ char *strdup(const std::string &text); - /** Trim leading and trailing whitespace. Like TRIM() in Fortran. - * - * \param line string that should be trimmed - * \return new string without whitespace (string) */ + /*! Trim leading and trailing whitespace. Like TRIM() in Fortran. + * + * \param line string that should be trimmed + * \return new string without whitespace (string) */ std::string trim(const std::string &line); - /** Return string with anything from '#' onward removed - * - * \param line string that should be trimmed - * \return new string without comment (string) */ + /*! Return string with anything from '#' onward removed + * + * \param line string that should be trimmed + * \return new string without comment (string) */ std::string trim_comment(const std::string &line); - /** Check if a string will likely have UTF-8 encoded characters - * - * UTF-8 uses the 7-bit standard ASCII table for the first 127 characters and - * all other characters are encoded as multiple bytes. For the multi-byte - * characters the first byte has either the highest two, three, or four bits - * set followed by a zero bit and followed by one, two, or three more bytes, - * respectively, where the highest bit is set and the second highest bit set - * to 0. The remaining bits combined are the character code, which is thus - * limited to 21-bits. - * - * For the sake of efficiency this test only checks if a character in the string - * has the highest bit set and thus is very likely an UTF-8 character. It will - * not be able to tell this this is a valid UTF-8 character or whether it is a - * 2-byte, 3-byte, or 4-byte character. - * + /*! Check if a string will likely have UTF-8 encoded characters + * + * UTF-8 uses the 7-bit standard ASCII table for the first 127 characters and + * all other characters are encoded as multiple bytes. For the multi-byte + * characters the first byte has either the highest two, three, or four bits + * set followed by a zero bit and followed by one, two, or three more bytes, + * respectively, where the highest bit is set and the second highest bit set + * to 0. The remaining bits combined are the character code, which is thus + * limited to 21-bits. + * + * For the sake of efficiency this test only checks if a character in the string + * has the highest bit set and thus is very likely an UTF-8 character. It will + * not be able to tell this this is a valid UTF-8 character or whether it is a + * 2-byte, 3-byte, or 4-byte character. + * \verbatim embed:rst *See also* :cpp:func:`utils::utf8_subst` \endverbatim - * \param line string that should be checked - * \return true if string contains UTF-8 encoded characters (bool) */ + * \param line string that should be checked + * \return true if string contains UTF-8 encoded characters (bool) */ inline bool has_utf8(const std::string &line) { @@ -316,243 +316,245 @@ namespace utils { return false; } - /** Replace known UTF-8 characters with ASCII equivalents - * + /*! Replace known UTF-8 characters with ASCII equivalents + * \verbatim embed:rst *See also* :cpp:func:`utils::has_utf8` \endverbatim - * \param line string that should be converted - * \return new string with ascii replacements (string) */ + * \param line string that should be converted + * \return new string with ascii replacements (string) */ std::string utf8_subst(const std::string &line); - /** Count words in string with custom choice of separating characters - * - * \param text string that should be searched - * \param separators string containing characters that will be treated as whitespace - * \return number of words found */ + /*! Count words in string with custom choice of separating characters + * + * \param text string that should be searched + * \param separators string containing characters that will be treated as whitespace + * \return number of words found */ size_t count_words(const std::string &text, const std::string &separators); - /** Count words in string, ignore any whitespace matching " \t\r\n\f" - * - * \param text string that should be searched - * \return number of words found */ + /*! Count words in string, ignore any whitespace matching " \t\r\n\f" + * + * \param text string that should be searched + * \return number of words found */ size_t count_words(const std::string &text); - /** Count words in C-string, ignore any whitespace matching " \t\r\n\f" - * - * \param text string that should be searched - * \return number of words found */ + /*! Count words in C-string, ignore any whitespace matching " \t\r\n\f" + * + * \param text string that should be searched + * \return number of words found */ size_t count_words(const char *text); - /** Count words in a single line, trim anything from '#' onward - * - * \param text string that should be trimmed and searched - * \param separators string containing characters that will be treated as whitespace - * \return number of words found */ + /*! Count words in a single line, trim anything from '#' onward + * + * \param text string that should be trimmed and searched + * \param separators string containing characters that will be treated as whitespace + * \return number of words found */ size_t trim_and_count_words(const std::string &text, const std::string &separators = " \t\r\n\f"); - /** Take text and split into non-whitespace words. - * - * This can handle strings with single and double quotes, escaped quotes, - * and escaped codes within quotes, but due to using an STL container and - * STL strings is rather slow because of making copies. Designed for - * parsing command lines and similar text and not for time critical - * processing. Use a tokenizer class if performance matters. - * + /*! Take text and split into non-whitespace words. + * + * This can handle strings with single and double quotes, escaped quotes, + * and escaped codes within quotes, but due to using an STL container and + * STL strings is rather slow because of making copies. Designed for + * parsing command lines and similar text and not for time critical + * processing. Use a tokenizer class if performance matters. + * \verbatim embed:rst *See also* :cpp:class:`Tokenizer`, :cpp:class:`ValueTokenizer` \endverbatim - * \param text string that should be split - * \return STL vector with the words */ + * \param text string that should be split + * \return STL vector with the words */ std::vector split_words(const std::string &text); - /** Take multi-line text and split into lines - * - * \param text string that should be split - * \return STL vector with the lines */ + /*! Take multi-line text and split into lines + * + * \param text string that should be split + * \return STL vector with the lines */ std::vector split_lines(const std::string &text); - /** Check if string can be converted to valid integer - * - * \param str string that should be checked - * \return true, if string contains valid a integer, false otherwise */ + /*! Check if string can be converted to valid integer + * + * \param str string that should be checked + * \return true, if string contains valid a integer, false otherwise */ bool is_integer(const std::string &str); - /** Check if string can be converted to valid floating-point number - * - * \param str string that should be checked - * \return true, if string contains valid number, false otherwise */ + /*! Check if string can be converted to valid floating-point number + * + * \param str string that should be checked + * \return true, if string contains valid number, false otherwise */ bool is_double(const std::string &str); - /** Check if string is a valid ID - * ID strings may contain only letters, numbers, and underscores. - * - * \param str string that should be checked - * \return true, if string contains valid id, false otherwise */ + /*! Check if string is a valid ID + * ID strings may contain only letters, numbers, and underscores. + * + * \param str string that should be checked + * \return true, if string contains valid id, false otherwise */ bool is_id(const std::string &str); - /** Try to detect pathname from FILE pointer. - * - * Currently only supported on Linux, otherwise will report "(unknown)". - * - * \param buf storage buffer for pathname. output will be truncated if not large enough - * \param len size of storage buffer. output will be truncated to this length - 1 - * \param fp FILE pointer struct from STDIO library for which we want to detect the name - * \return pointer to the storage buffer, i.e. buf */ + /*! Try to detect pathname from FILE pointer. + * + * Currently only supported on Linux, otherwise will report "(unknown)". + * + * \param buf storage buffer for pathname. output will be truncated if not large enough + * \param len size of storage buffer. output will be truncated to this length - 1 + * \param fp FILE pointer struct from STDIO library for which we want to detect the name + * \return pointer to the storage buffer, i.e. buf */ const char *guesspath(char *buf, int len, FILE *fp); - /** Strip off leading part of path, return just the filename - * - * \param path file path - * \return file name */ + /*! Strip off leading part of path, return just the filename + * + * \param path file path + * \return file name */ std::string path_basename(const std::string &path); - /** Return the directory part of a path. Return "." if empty - * - * \param path file path - * \return directory name */ + /*! Return the directory part of a path. Return "." if empty + * + * \param path file path + * \return directory name */ std::string path_dirname(const std::string &path); - /** Join two pathname segments - * - * This uses the forward slash '/' character unless LAMMPS is compiled - * for Windows where it used the equivalent backward slash '\\'. - * - * \param a first path - * \param b second path - * \return combined path */ + /*! Join two pathname segments + * + * This uses the forward slash '/' character unless LAMMPS is compiled + * for Windows where it used the equivalent backward slash '\\'. + * + * \param a first path + * \param b second path + * \return combined path */ std::string path_join(const std::string &a, const std::string &b); - /** Check if file exists and is readable - * - * \param path file path - * \return true if file exists and is readable */ + /*! Check if file exists and is readable + * + * \param path file path + * \return true if file exists and is readable */ bool file_is_readable(const std::string &path); - /** Determine full path of potential file. If file is not found in current directory, - * search directories listed in LAMMPS_POTENTIALS environment variable - * - * \param path file path - * \return full path to potential file */ + /*! Determine full path of potential file. If file is not found in current directory, + * search directories listed in LAMMPS_POTENTIALS environment variable + * + * \param path file path + * \return full path to potential file */ std::string get_potential_file_path(const std::string &path); - /** Read potential file and return DATE field if it is present - * - * \param path file path - * \param potential_name name of potential that is being read - * \return DATE field if present */ + /*! Read potential file and return DATE field if it is present + * + * \param path file path + * \param potential_name name of potential that is being read + * \return DATE field if present */ std::string get_potential_date(const std::string &path, const std::string &potential_name); - /** Read potential file and return UNITS field if it is present - * - * \param path file path - * \param potential_name name of potential that is being read - * \return UNITS field if present */ + /*! Read potential file and return UNITS field if it is present + * + * \param path file path + * \param potential_name name of potential that is being read + * \return UNITS field if present */ std::string get_potential_units(const std::string &path, const std::string &potential_name); enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1 << 1 }; enum { UNKNOWN = 0, ENERGY }; - /** Return bitmask of available conversion factors for a given property - * - * \param property property to be converted - * \return bitmask indicating available conversions */ + /*! Return bitmask of available conversion factors for a given property + * + * \param property property to be converted + * \return bitmask indicating available conversions */ + int get_supported_conversions(const int property); - /** Return unit conversion factor for given property and selected from/to units - * - * \param property property to be converted - * \param conversion constant indicating the conversion - * \return conversion factor */ + /*! Return unit conversion factor for given property and selected from/to units + * + * \param property property to be converted + * \param conversion constant indicating the conversion + * \return conversion factor */ double get_conversion_factor(const int property, const int conversion); - /** Open a potential file as specified by *name* - * - * If opening the file directly fails, the function will search for - * it in the list of folder pointed to by the environment variable - * ``LAMMPS_POTENTIALS`` (if it is set). - * - * If the potential file has a ``UNITS`` tag in the first line, the - * tag's value is compared to the current unit style setting. - * The behavior of the function then depends on the value of the - * *auto_convert* parameter. If it is a null pointer, then the unit - * values must match or else the open will fail with an error. Otherwise - * the bitmask that *auto_convert* points to is used check for - * compatibility with possible automatic conversions by the calling - * function. If compatible, the bitmask is set to the required - * conversion or ``utils::NOCONVERT``. - * - * \param name file- or pathname of the potential file - * \param lmp pointer to top-level LAMMPS class instance - * \param auto_convert pointer to unit conversion bitmask or ``nullptr`` - * \return FILE pointer of the opened potential file or ``nullptr`` */ + /*! Open a potential file as specified by *name* + * + * If opening the file directly fails, the function will search for + * it in the list of folder pointed to by the environment variable + * ``LAMMPS_POTENTIALS`` (if it is set). + * + * If the potential file has a ``UNITS`` tag in the first line, the + * tag's value is compared to the current unit style setting. + * The behavior of the function then depends on the value of the + * *auto_convert* parameter. If it is a null pointer, then the unit + * values must match or else the open will fail with an error. Otherwise + * the bitmask that *auto_convert* points to is used check for + * compatibility with possible automatic conversions by the calling + * function. If compatible, the bitmask is set to the required + * conversion or ``utils::NOCONVERT``. + * + * \param name file- or pathname of the potential file + * \param lmp pointer to top-level LAMMPS class instance + * \param auto_convert pointer to unit conversion bitmask or ``nullptr`` + * \return FILE pointer of the opened potential file or ``nullptr`` */ FILE *open_potential(const std::string &name, LAMMPS *lmp, int *auto_convert); - /** Convert a time string to seconds - * - * The strings "off" and "unlimited" result in -1 - * - * \param timespec a string in the following format: ([[HH:]MM:]SS) - * \return total in seconds */ + /*! Convert a time string to seconds + * + * The strings "off" and "unlimited" result in -1 + * + * \param timespec a string in the following format: ([[HH:]MM:]SS) + * \return total in seconds */ double timespec2seconds(const std::string ×pec); - /** Convert a LAMMPS version date to a number - * - * This will generate a number YYYYMMDD from a date string - * (with or without blanks) that is suitable for numerical - * comparisons, i.e. later dates will generate a larger number. - * - * The day may or may not have a leading zero, the month - * is identified by the first 3 letters (so there may be more) - * and the year may be 2 or 4 digits (the missing 2 digits will - * be assumed as 20. That is 04 corresponds to 2004). - * - * No check is made whether the date is valid. - * - * \param date string in the format (Day Month Year) - * \return date code */ + /*! Convert a LAMMPS version date to a number + * + * This will generate a number YYYYMMDD from a date string + * (with or without blanks) that is suitable for numerical + * comparisons, i.e. later dates will generate a larger number. + * + * The day may or may not have a leading zero, the month + * is identified by the first 3 letters (so there may be more) + * and the year may be 2 or 4 digits (the missing 2 digits will + * be assumed as 20. That is 04 corresponds to 2004). + * + * No check is made whether the date is valid. + * + * \param date string in the format (Day Month Year) + * \return date code */ + int date2num(const std::string &date); - /** Custom merge sort implementation - * - * This function provides a custom upward hybrid merge sort - * implementation with support to pass an opaque pointer to - * the comparison function, e.g. for access to class members. - * This avoids having to use global variables. For improved - * performance, it uses an in-place insertion sort on initial - * chunks of up to 64 elements and switches to merge sort from - * then on. - * - * \param index Array with indices to be sorted - * \param num Length of the index array - * \param ptr Pointer to opaque object passed to comparison function - * \param comp Pointer to comparison function */ + /*! Custom merge sort implementation + * + * This function provides a custom upward hybrid merge sort + * implementation with support to pass an opaque pointer to + * the comparison function, e.g. for access to class members. + * This avoids having to use global variables. For improved + * performance, it uses an in-place insertion sort on initial + * chunks of up to 64 elements and switches to merge sort from + * then on. + * + * \param index Array with indices to be sorted + * \param num Length of the index array + * \param ptr Pointer to opaque object passed to comparison function + * \param comp Pointer to comparison function */ void merge_sort(int *index, int num, void *ptr, int (*comp)(int, int, void *)); } // namespace utils diff --git a/unittest/force-styles/test_pair_style.cpp b/unittest/force-styles/test_pair_style.cpp index 939bf5dbb5..b0221cf56e 100644 --- a/unittest/force-styles/test_pair_style.cpp +++ b/unittest/force-styles/test_pair_style.cpp @@ -158,7 +158,7 @@ void run_lammps(LAMMPS *lmp) command("run 4 post no"); } -void restart_lammps(LAMMPS *lmp, const TestConfig &cfg) +void restart_lammps(LAMMPS *lmp, const TestConfig &cfg, bool nofdotr = false, bool newton = true) { // utility lambda to improve readability auto command = [&](const std::string &line) { @@ -166,6 +166,10 @@ void restart_lammps(LAMMPS *lmp, const TestConfig &cfg) }; command("clear"); + if (newton) + command("newton on"); + else + command("newton off"); command("read_restart " + cfg.basename + ".restart"); if (!lmp->force->pair) { @@ -180,6 +184,7 @@ void restart_lammps(LAMMPS *lmp, const TestConfig &cfg) for (auto &post_command : cfg.post_commands) { command(post_command); } + if (nofdotr) command("pair_modify nofdotr"); command("run 0 post no"); } @@ -518,6 +523,41 @@ TEST(PairStyle, plain) EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon); if (print_stats) std::cerr << "restart_energy stats:" << stats << std::endl; + // pair style rann does not support pair_modify nofdotr + // temporarily disable testing pair style reax/c until we merge the other fixes + if ((test_config.pair_style != "rann") && !(test_config.pair_style.find("reax") != std::string::npos)) { + if (!verbose) ::testing::internal::CaptureStdout(); + restart_lammps(lmp, test_config, true); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + f = lmp->atom->f; + tag = lmp->atom->tag; + stats.reset(); + ASSERT_EQ(nlocal + 1, f_ref.size()); + for (int i = 0; i < nlocal; ++i) { + EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon); + } + if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl; + + pair = lmp->force->pair; + stress = pair->virial; + stats.reset(); + EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon); + EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon); + EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon); + EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon); + EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon); + EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon); + if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl; + + stats.reset(); + EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon); + EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon); + if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl; + } + if (!verbose) ::testing::internal::CaptureStdout(); data_lammps(lmp, test_config); if (!verbose) ::testing::internal::GetCapturedStdout(); @@ -704,14 +744,14 @@ TEST(PairStyle, omp) EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon); if (print_stats) std::cerr << "run_energy stats, newton on: " << stats << std::endl; - if (!verbose) ::testing::internal::CaptureStdout(); - cleanup_lammps(lmp, test_config); - lmp = init_lammps(argc, argv, test_config, false); - if (!verbose) ::testing::internal::GetCapturedStdout(); - // skip over these tests if newton pair is forced to be on if (lmp->force->newton_pair == 0) { + if (!verbose) ::testing::internal::CaptureStdout(); + cleanup_lammps(lmp, test_config); + lmp = init_lammps(argc, argv, test_config, false); + if (!verbose) ::testing::internal::GetCapturedStdout(); + f = lmp->atom->f; tag = lmp->atom->tag; stats.reset(); @@ -770,6 +810,38 @@ TEST(PairStyle, omp) EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon); if (print_stats) std::cerr << "run_energy stats, newton off:" << stats << std::endl; } + + if (!verbose) ::testing::internal::CaptureStdout(); + restart_lammps(lmp, test_config, true); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + f = lmp->atom->f; + tag = lmp->atom->tag; + stats.reset(); + ASSERT_EQ(nlocal + 1, f_ref.size()); + for (int i = 0; i < nlocal; ++i) { + EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon); + } + if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl; + + pair = lmp->force->pair; + stress = pair->virial; + stats.reset(); + EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon); + if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl; + + stats.reset(); + EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon); + if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl; + if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); if (!verbose) ::testing::internal::GetCapturedStdout(); @@ -781,8 +853,10 @@ TEST(PairStyle, gpu) if (!Info::has_gpu_device()) GTEST_SKIP(); if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP(); - const char *args_neigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "gpu"}; - const char *args_noneigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "gpu", "-pk", "gpu", "0", "neigh", "no"}; + const char *args_neigh[] = {"PairStyle", "-log", "none", "-echo", + "screen", "-nocite", "-sf", "gpu"}; + const char *args_noneigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", + "gpu", "-pk", "gpu", "0", "neigh", "no"}; char **argv = (char **)args_neigh; int argc = sizeof(args_neigh) / sizeof(char *); @@ -790,7 +864,7 @@ TEST(PairStyle, gpu) // cannot use GPU neighbor list with hybrid pair style (yet) if (test_config.pair_style.substr(0, 6) == "hybrid") { argv = (char **)args_noneigh; - argc = sizeof(args_noneigh) / sizeof(char *); + argc = sizeof(args_noneigh) / sizeof(char *); } ::testing::internal::CaptureStdout(); @@ -817,14 +891,17 @@ TEST(PairStyle, gpu) // relax error for GPU package depending on precision setting double epsilon = test_config.epsilon; - if (Info::has_accelerator_feature("GPU","precision","double")) + if (Info::has_accelerator_feature("GPU", "precision", "double")) epsilon *= 7.5; - else if (Info::has_accelerator_feature("GPU","precision","mixed")) + else if (Info::has_accelerator_feature("GPU", "precision", "mixed")) epsilon *= 5.0e8; - else epsilon *= 1.0e10; - // relax test precision when using pppm and single precision FFTs, but only when also running with double precision + else + epsilon *= 1.0e10; + // relax test precision when using pppm and single precision FFTs, but only when also + // running with double precision #if defined(FFT_SINGLE) - if (lmp->force->kspace && lmp->force->kspace->compute_flag && Info::has_accelerator_feature("GPU","precision","double")) + if (lmp->force->kspace && lmp->force->kspace->compute_flag && + Info::has_accelerator_feature("GPU", "precision", "double")) if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8; #endif const std::vector &f_ref = test_config.init_forces; @@ -910,7 +987,7 @@ TEST(PairStyle, intel) int argc = sizeof(args) / sizeof(char *); ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(argc, argv, test_config); + LAMMPS *lmp = init_lammps(argc, argv, test_config, true); std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -1124,6 +1201,37 @@ TEST(PairStyle, opt) EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon); if (print_stats) std::cerr << "run_energy stats:" << stats << std::endl; + if (!verbose) ::testing::internal::CaptureStdout(); + restart_lammps(lmp, test_config, true); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + f = lmp->atom->f; + tag = lmp->atom->tag; + stats.reset(); + ASSERT_EQ(nlocal + 1, f_ref.size()); + for (int i = 0; i < nlocal; ++i) { + EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon); + } + if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl; + + pair = lmp->force->pair; + stress = pair->virial; + stats.reset(); + EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon); + EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon); + if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl; + + stats.reset(); + EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon); + EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon); + if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl; + if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); if (!verbose) ::testing::internal::GetCapturedStdout(); diff --git a/unittest/force-styles/tests/atomic-pair-polymorphic_eam.yaml b/unittest/force-styles/tests/atomic-pair-polymorphic_eam.yaml index f3ad72b9cd..6a2c16a80d 100644 --- a/unittest/force-styles/tests/atomic-pair-polymorphic_eam.yaml +++ b/unittest/force-styles/tests/atomic-pair-polymorphic_eam.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:09:04 2021 -epsilon: 1e-14 +epsilon: 5e-14 prerequisites: ! | pair polymorphic pre_commands: ! | diff --git a/unittest/force-styles/tests/atomic-pair-reax_c.yaml b/unittest/force-styles/tests/atomic-pair-reax_c.yaml index a0dabff323..ed0908beae 100644 --- a/unittest/force-styles/tests/atomic-pair-reax_c.yaml +++ b/unittest/force-styles/tests/atomic-pair-reax_c.yaml @@ -2,6 +2,7 @@ lammps_version: 24 Aug 2020 date_generated: Tue Sep 15 09:44:24 202 epsilon: 5e-11 +skip_tests: omp prerequisites: ! | pair reax/c fix qeq/reax diff --git a/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml b/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml index 567bddf5ce..bf71811de9 100644 --- a/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml +++ b/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml @@ -2,6 +2,7 @@ lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:09:06 2021 epsilon: 2e-10 +skip_tests: omp prerequisites: ! | pair reax/c fix qeq/reax diff --git a/unittest/force-styles/tests/atomic-pair-reax_c_noqeq.yaml b/unittest/force-styles/tests/atomic-pair-reax_c_noqeq.yaml index efdf3ff5de..8f4b8b81c3 100644 --- a/unittest/force-styles/tests/atomic-pair-reax_c_noqeq.yaml +++ b/unittest/force-styles/tests/atomic-pair-reax_c_noqeq.yaml @@ -2,6 +2,7 @@ lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:09:08 2021 epsilon: 5e-13 +skip_tests: omp prerequisites: ! | pair reax/c pre_commands: ! | diff --git a/unittest/force-styles/tests/manybody-pair-bop.yaml b/unittest/force-styles/tests/manybody-pair-bop.yaml index 3855d2d263..bc5d198aaf 100644 --- a/unittest/force-styles/tests/manybody-pair-bop.yaml +++ b/unittest/force-styles/tests/manybody-pair-bop.yaml @@ -1,7 +1,7 @@ --- lammps_version: 8 Apr 2021 date_generated: Wed May 5 11:50:15 2021 -epsilon: 1e-14 +epsilon: 2e-13 prerequisites: ! | pair bop pre_commands: ! | diff --git a/unittest/force-styles/tests/manybody-pair-bop_save.yaml b/unittest/force-styles/tests/manybody-pair-bop_save.yaml index 5f13b370a4..9c8ca13067 100644 --- a/unittest/force-styles/tests/manybody-pair-bop_save.yaml +++ b/unittest/force-styles/tests/manybody-pair-bop_save.yaml @@ -1,7 +1,7 @@ --- lammps_version: 8 Apr 2021 date_generated: Wed May 5 11:50:24 2021 -epsilon: 1e-14 +epsilon: 2e-14 prerequisites: ! | pair bop pre_commands: ! | diff --git a/unittest/force-styles/tests/manybody-pair-edip_multi.yaml b/unittest/force-styles/tests/manybody-pair-edip_multi.yaml index ad2dc9018e..d86aa02b3f 100644 --- a/unittest/force-styles/tests/manybody-pair-edip_multi.yaml +++ b/unittest/force-styles/tests/manybody-pair-edip_multi.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:09:15 2021 -epsilon: 5e-14 +epsilon: 7.5e-14 prerequisites: ! | pair edip/multi pre_commands: ! | diff --git a/unittest/force-styles/tests/manybody-pair-snap.yaml b/unittest/force-styles/tests/manybody-pair-snap.yaml index 139d1b94d0..57e0e9d599 100644 --- a/unittest/force-styles/tests/manybody-pair-snap.yaml +++ b/unittest/force-styles/tests/manybody-pair-snap.yaml @@ -14,7 +14,8 @@ pair_style: hybrid/overlay zbl 4.0 4.8 snap pair_coeff: ! | 1*8 1*8 zbl 73 73 * * snap Ta06A.snapcoeff Ta06A.snapparam Ta Ta Ta Ta Ta Ta Ta Ta -extract: ! "" +extract: ! | + scale 2 natoms: 64 init_vdwl: -473.569864629026 init_coul: 0 diff --git a/unittest/force-styles/tests/manybody-pair-snap_chem.yaml b/unittest/force-styles/tests/manybody-pair-snap_chem.yaml index 39d2abd804..b97709863e 100644 --- a/unittest/force-styles/tests/manybody-pair-snap_chem.yaml +++ b/unittest/force-styles/tests/manybody-pair-snap_chem.yaml @@ -16,7 +16,8 @@ pair_coeff: ! | 1*4 5*8 zbl 49 15 5*8 5*8 zbl 15 15 * * snap InP_JCPA2020.snapcoeff InP_JCPA2020.snapparam In In In In P P P P -extract: ! "" +extract: ! | + scale 2 natoms: 64 init_vdwl: -185.3871232982 init_coul: 0 diff --git a/unittest/force-styles/tests/mol-pair-lj_charmm_coul_charmm_implicit.yaml b/unittest/force-styles/tests/mol-pair-lj_charmm_coul_charmm_implicit.yaml index 91a73b31e3..d9bd7786ba 100644 --- a/unittest/force-styles/tests/mol-pair-lj_charmm_coul_charmm_implicit.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_charmm_coul_charmm_implicit.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:08:45 2021 -epsilon: 2e-13 +epsilon: 2.5e-13 prerequisites: ! | atom full pair lj/charmm/coul/charmm/implicit diff --git a/unittest/force-styles/tests/mol-pair-lj_charmm_coul_long.yaml b/unittest/force-styles/tests/mol-pair-lj_charmm_coul_long.yaml index c58a4849da..eaa7ca4fa5 100644 --- a/unittest/force-styles/tests/mol-pair-lj_charmm_coul_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_charmm_coul_long.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:08:45 2021 -epsilon: 7.5e-14 +epsilon: 1e-13 prerequisites: ! | atom full pair lj/charmm/coul/long diff --git a/unittest/force-styles/tests/mol-pair-lj_cut_coul_long_cs.yaml b/unittest/force-styles/tests/mol-pair-lj_cut_coul_long_cs.yaml index 863b41f597..0dc745ca26 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cut_coul_long_cs.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cut_coul_long_cs.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Feb 2021 date_generated: Fri Feb 26 23:08:48 2021 -epsilon: 7.5e-14 +epsilon: 1e-13 prerequisites: ! | atom full pair lj/cut/coul/long/cs diff --git a/unittest/formats/test_file_operations.cpp b/unittest/formats/test_file_operations.cpp index eb30b69cfd..794222fe5a 100644 --- a/unittest/formats/test_file_operations.cpp +++ b/unittest/formats/test_file_operations.cpp @@ -208,11 +208,11 @@ TEST_F(FileOperationsTest, read_lines_from_file) MPI_Comm world = MPI_COMM_WORLD; int me, rv; memset(buf, 0, MAX_BUF_SIZE); + MPI_Comm_rank(world, &me); rv = utils::read_lines_from_file(nullptr, 1, MAX_BUF_SIZE, buf, me, world); ASSERT_EQ(rv, 1); - MPI_Comm_rank(world, &me); if (me == 0) { fp = fopen("safe_file_read_test.txt", "r"); ASSERT_NE(fp, nullptr);