diff --git a/cmake/Modules/Packages/GPU.cmake b/cmake/Modules/Packages/GPU.cmake index fe15917f47..aec8887c30 100644 --- a/cmake/Modules/Packages/GPU.cmake +++ b/cmake/Modules/Packages/GPU.cmake @@ -30,7 +30,15 @@ file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/[^.]*.cpp) file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu) if(GPU_API STREQUAL "CUDA") - find_package(CUDA REQUIRED) + find_package(CUDA QUIET) + # augment search path for CUDA toolkit libraries to include the stub versions. Needed to find libcuda.so on machines without a CUDA driver installation + if(CUDA_FOUND) + set(CMAKE_LIBRARY_PATH "${CUDA_TOOLKIT_ROOT_DIR}/lib64/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib64;${CUDA_TOOLKIT_ROOT_DIR}/lib;${CMAKE_LIBRARY_PATH}") + find_package(CUDA REQUIRED) + else() + message(FATAL_ERROR "CUDA Toolkit not found") + endif() + find_program(BIN2C bin2c) if(NOT BIN2C) message(FATAL_ERROR "Could not find bin2c, use -DBIN2C=/path/to/bin2c to help cmake finding it.") diff --git a/doc/src/Build_settings.rst b/doc/src/Build_settings.rst index b68313aaed..6053bdbf8a 100644 --- a/doc/src/Build_settings.rst +++ b/doc/src/Build_settings.rst @@ -4,15 +4,15 @@ Optional build settings LAMMPS can be built with several optional settings. Each sub-section explain how to do this for building both with CMake and make. -* :ref:`C++11 standard compliance ` when building all of LAMMPS -* :ref:`FFT library ` for use with the :doc:`kspace_style pppm ` command -* :ref:`Size of LAMMPS integer types ` -* :ref:`Read or write compressed files ` -* :ref:`Output of JPG and PNG files ` via the :doc:`dump image ` command -* :ref:`Output of movie files ` via the :doc:`dump_movie ` command -* :ref:`Memory allocation alignment ` -* :ref:`Workaround for long long integers ` -* :ref:`Error handling exceptions ` when using LAMMPS as a library +* `C++11 standard compliance`_ when building all of LAMMPS +* `FFT library`_ for use with the :doc:`kspace_style pppm ` command +* `Size of LAMMPS integer types and size limits`_ +* `Read or write compressed files`_ +* `Output of JPG, PNG, and move files` via the :doc:`dump image ` or :doc:`dump movie ` commands +* `Memory allocation alignment`_ +* `Workaround for long long integers`_ +* `Exception handling when using LAMMPS as a library`_ to capture errors +* `Trigger selected floating-point exceptions`_ ---------- diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 590b3d2ea8..6ad5ed4435 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -142,11 +142,11 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`sph/t/atom ` * :doc:`spin ` * :doc:`stress/atom ` - * :doc:`stress/cartesian ` - * :doc:`stress/cylinder ` + * :doc:`stress/cartesian ` + * :doc:`stress/cylinder ` * :doc:`stress/mop ` * :doc:`stress/mop/profile ` - * :doc:`stress/spherical ` + * :doc:`stress/spherical ` * :doc:`stress/tally ` * :doc:`tdpd/cc/atom ` * :doc:`temp (k) ` diff --git a/doc/src/Developer_notes.rst b/doc/src/Developer_notes.rst index cd3c51fed6..e58dec94ac 100644 --- a/doc/src/Developer_notes.rst +++ b/doc/src/Developer_notes.rst @@ -25,18 +25,17 @@ following benefits: - transparent support for translating unsupported UTF-8 characters to their ASCII equivalents (the text-to-value conversion functions **only** accept ASCII characters) -In most cases (e.g. potential files) the same data is needed on all -MPI ranks. Then it is best to do the reading and parsing only on MPI -rank 0, and communicate the data later with one or more -``MPI_Bcast()`` calls. For reading generic text and potential -parameter files the custom classes :cpp:class:`TextFileReader -` and :cpp:class:`PotentialFileReader -` are available. They allow reading -the file as individual lines for which they can return a tokenizer -class (see below) for parsing the line. Or they can return blocks of -numbers as a vector directly. The documentation on `File reader -classes `_ contains an example for a typical -case. +In most cases (e.g. potential files) the same data is needed on all MPI +ranks. Then it is best to do the reading and parsing only on MPI rank +0, and communicate the data later with one or more ``MPI_Bcast()`` +calls. For reading generic text and potential parameter files the +custom classes :cpp:class:`TextFileReader ` +and :cpp:class:`PotentialFileReader ` +are available. They allow reading the file as individual lines for which +they can return a tokenizer class (see below) for parsing the line. Or +they can return blocks of numbers as a vector directly. The +documentation on :ref:`File reader classes ` +contains an example for a typical case. When reading per-atom data, the data on each line of the file usually needs to include an atom ID so it can be associated with a particular diff --git a/doc/src/Developer_org.rst b/doc/src/Developer_org.rst index ce36d9a590..63ff30abfe 100644 --- a/doc/src/Developer_org.rst +++ b/doc/src/Developer_org.rst @@ -252,12 +252,6 @@ follows: - The Timer class logs timing information, output at the end of a run. -.. TODO section on "Spatial decomposition and parallel operations" -.. diagram of 3d processor grid, brick vs. tiled. local vs. ghost -.. atoms, 6-way communication with pack/unpack functions, -.. PBC as part of the communication, forward and reverse communication -.. rendezvous communication, ring communication. - .. TODO section on "Fixes, Computes, and Variables" .. how and when data is computed and provided and how it is .. referenced. flags in Fix/Compute/Variable classes tell diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index a9df85c899..720eececcc 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -396,7 +396,7 @@ A typical code segment would look like this: ---------- -.. file-reader-classes: +.. _file-reader-classes: File reader classes ------------------- diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index 1b4b82a5e7..396ecd6442 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -9,7 +9,7 @@ gives links to documentation, example scripts, and pictures/movies (if available) that illustrate use of the package. The majority of packages can be included in a LAMMPS build with a -single setting (``-D PGK_=on`` for CMake) or command +single setting (``-D PKG_=on`` for CMake) or command (``make yes-`` for make). See the :doc:`Build package ` page for more info. A few packages may require additional steps; this is indicated in the descriptions below. The :doc:`Build extras ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6973559d16..2768543179 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -288,11 +288,11 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`sph/t/atom ` - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms * :doc:`spin ` - magnetic quantities for a system of atoms having spins * :doc:`stress/atom ` - stress tensor for each atom -* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates -* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates +* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates +* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates * :doc:`stress/mop ` - normal components of the local stress tensor using the method of planes * :doc:`stress/mop/profile ` - profile of the normal components of the local stress tensor using the method of planes -* :doc:`stress/spherical ` - stress tensor in spherical coordinates +* :doc:`stress/spherical ` - stress tensor in spherical coordinates * :doc:`stress/tally ` - stress between two groups of atoms via the tally callback mechanism * :doc:`tdpd/cc/atom ` - per-atom chemical concentration of a specified species for each tDPD particle * :doc:`temp ` - temperature of group of atoms diff --git a/doc/src/compute_momentum.rst b/doc/src/compute_momentum.rst index c86901933b..40d2f527f5 100644 --- a/doc/src/compute_momentum.rst +++ b/doc/src/compute_momentum.rst @@ -23,11 +23,10 @@ Examples Description """"""""""" -Define a computation that calculates the translational momentum -of a group of particles. - -The momentum of each particles is computed as m v, where m and v are -the mass and velocity of the particle. +Define a computation that calculates the translational momentum *p* +of a group of particles. It is computed as the sum :math:`\vec{p} = \sum_i m_i \cdot \vec{v}_i` +over all particles in the compute group, where *m* and *v* are +the mass and velocity vector of the particle, respectively. Output info """"""""""" diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 18b9e703c5..d439e18d34 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -73,7 +73,7 @@ command, since those are contributions to the global system pressure. NOTE 3: The local stress profile generated by compute *stress/mop/profile* is similar to that obtained by compute -:doc:`stress/cartesian `. +:doc:`stress/cartesian `. A key difference is that compute *stress/mop/profile* considers particles crossing a set of planes, @@ -122,7 +122,7 @@ intra-molecular interactions, and long range (kspace) interactions. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` +:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` Default """"""" diff --git a/doc/src/compute_stress_cartesian.rst b/doc/src/compute_stress_profile.rst similarity index 100% rename from doc/src/compute_stress_cartesian.rst rename to doc/src/compute_stress_profile.rst diff --git a/doc/src/compute_tally.rst b/doc/src/compute_tally.rst index 313b54645c..45ebb42843 100644 --- a/doc/src/compute_tally.rst +++ b/doc/src/compute_tally.rst @@ -31,7 +31,7 @@ Syntax compute ID group-ID style group2-ID * ID, group-ID are documented in :doc:`compute ` command -* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or * or *pe/tally* or *pe/mol/tally* or *stress/tally* +* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or *pe/tally* or *pe/mol/tally* or *stress/tally* * group2-ID = group ID of second (or same) group Examples @@ -61,7 +61,7 @@ mechanism. Compute *pe/mol/tally* is one such style, that can - through using this mechanism - separately tally intermolecular and intramolecular energies. Something that would otherwise be impossible without integrating this as a core functionality into -the based classes of LAMMPS. +the base classes of LAMMPS. ---------- @@ -148,30 +148,38 @@ pairwise property computations. Output info """"""""""" -Compute *pe/tally* calculates a global scalar (the energy) and a per -atom scalar (the contributions of the single atom to the global -scalar). Compute *pe/mol/tally* calculates a global 4-element vector -containing (in this order): *evdwl* and *ecoul* for intramolecular pairs -and *evdwl* and *ecoul* for intermolecular pairs. Since molecules are -identified by their molecule IDs, the partitioning does not have to be -related to molecules, but the energies are tallied into the respective -slots depending on whether the molecule IDs of a pair are the same or -different. Compute *force/tally* calculates a global scalar (the force -magnitude) and a per atom 3-element vector (force contribution from -each atom). Compute *stress/tally* calculates a global scalar -(average of the diagonal elements of the stress tensor) and a per atom -vector (the 6 elements of stress tensor contributions from the -individual atom). As in :doc:`compute heat/flux `, -compute *heat/flux/tally* calculates a global vector of length 6, -where the first 3 components are the :math:`x`, :math:`y`, :math:`z` -components of the full heat flow vector, -and the next 3 components are the corresponding components -of just the convective portion of the flow, i.e. the -first term in the equation for :math:`\mathbf{Q}`. -Compute *heat/flux/virial/tally* calculates a global scalar (heat flow) -and a per atom 3-element vector -(contribution to the force acting over atoms in the first group -from individual atoms in both groups). +- Compute *pe/tally* calculates a global scalar (the energy) and a per + atom scalar (the contributions of the single atom to the global + scalar). + +- Compute *pe/mol/tally* calculates a global 4-element vector containing + (in this order): *evdwl* and *ecoul* for intramolecular pairs and + *evdwl* and *ecoul* for intermolecular pairs. Since molecules are + identified by their molecule IDs, the partitioning does not have to be + related to molecules, but the energies are tallied into the respective + slots depending on whether the molecule IDs of a pair are the same or + different. + +- Compute *force/tally* calculates a global scalar (the force magnitude) + and a per atom 3-element vector (force contribution from each atom). + +- Compute *stress/tally* calculates a global scalar + (average of the diagonal elements of the stress tensor) and a per atom + vector (the 6 elements of stress tensor contributions from the + individual atom). + +- As in :doc:`compute heat/flux `, + compute *heat/flux/tally* calculates a global vector of length 6, + where the first 3 components are the :math:`x`, :math:`y`, :math:`z` + components of the full heat flow vector, + and the next 3 components are the corresponding components + of just the convective portion of the flow, i.e. the + first term in the equation for :math:`\mathbf{Q}`. + +- Compute *heat/flux/virial/tally* calculates a global scalar (heat flow) + and a per atom 3-element vector + (contribution to the force acting over atoms in the first group + from individual atoms in both groups). Both the scalar and vector values calculated by this compute are "extensive". diff --git a/src/.gitignore b/src/.gitignore index 8b60afb189..8803d8a7e3 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1,6 +1,7 @@ /Makefile.package /Makefile.package.settings /MAKE/MINE +/ADIOS/Makefile.lammps /Make.py.last /lmp_* @@ -267,6 +268,8 @@ /fix_drag.h /fix_numdiff.cpp /fix_numdiff.h +/fix_numdiff_virial.cpp +/fix_numdiff_virial.h /fix_spring_rg.cpp /fix_spring_rg.h /fix_temp_csld.cpp @@ -934,6 +937,8 @@ /msm_cg.h /neb.cpp /neb.h +/netcdf_units.cpp +/netcdf_units.h /pair_adp.cpp /pair_adp.h /pair_agni.cpp @@ -1069,6 +1074,8 @@ /pair_hdnnp.h /pair_ilp_graphene_hbn.cpp /pair_ilp_graphene_hbn.h +/pair_ilp_tmd.cpp +/pair_ilp_tmd.h /pair_kolmogorov_crespi_full.cpp /pair_kolmogorov_crespi_full.h /pair_kolmogorov_crespi_z.cpp @@ -1175,8 +1182,8 @@ /pair_nm_cut_coul_cut.h /pair_nm_cut_coul_long.cpp /pair_nm_cut_coul_long.h -/pait_nm_cut_split.cpp -/pait_nm_cut_split.h +/pair_nm_cut_split.cpp +/pair_nm_cut_split.h /pair_oxdna_*.cpp /pair_oxdna_*.h /pair_oxdna2_*.cpp @@ -1202,6 +1209,8 @@ /pair_rebo.h /pair_resquared.cpp /pair_resquared.h +/pair_saip_metal.cpp +/pair_saip_metal.h /pair_sdpd_taitwater_isothermal.cpp /pair_sdpd_taitwater_isothermal.h /pair_sph_heatconduction.cpp diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp index 7762e261cd..066869e851 100644 --- a/src/CLASS2/pair_lj_class2.cpp +++ b/src/CLASS2/pair_lj_class2.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -65,13 +64,13 @@ PairLJClass2::~PairLJClass2() void PairLJClass2::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -104,34 +103,32 @@ void PairLJClass2::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -144,10 +141,10 @@ void PairLJClass2::compute(int eflag, int vflag) void PairLJClass2::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -165,8 +162,8 @@ void PairLJClass2::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -187,28 +184,28 @@ void PairLJClass2::compute_inner() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -219,10 +216,10 @@ void PairLJClass2::compute_inner() void PairLJClass2::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -243,10 +240,10 @@ void PairLJClass2::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -267,32 +264,32 @@ void PairLJClass2::compute_middle() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -303,13 +300,13 @@ void PairLJClass2::compute_middle() void PairLJClass2::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -327,8 +324,8 @@ void PairLJClass2::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -349,55 +346,54 @@ void PairLJClass2::compute_outer(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } if (eflag) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (vflag) { if (rsq <= cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; } else if (rsq < cut_in_on_sq) - fpair = factor_lj*forcelj*r2inv; + fpair = factor_lj * forcelj * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -409,23 +405,21 @@ void PairLJClass2::compute_outer(int eflag, int vflag) void PairLJClass2::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(cut, np1, np1, "pair:cut"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -434,14 +428,14 @@ void PairLJClass2::allocate() void PairLJClass2::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; @@ -454,22 +448,22 @@ void PairLJClass2::settings(int narg, char **arg) void PairLJClass2::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); double cut_one = cut_global; - if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp); + if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; @@ -478,7 +472,7 @@ void PairLJClass2::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -492,12 +486,12 @@ void PairLJClass2::init_style() int irequest; int respa = 0; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; } - irequest = neighbor->request(this,instance_me); + irequest = neighbor->request(this, instance_me); if (respa >= 1) { neighbor->requests[irequest]->respaouter = 1; @@ -507,10 +501,11 @@ void PairLJClass2::init_style() // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; } /* ---------------------------------------------------------------------- @@ -523,23 +518,22 @@ double PairLJClass2::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; @@ -550,7 +544,7 @@ double PairLJClass2::init_one(int i, int j) // check interior rRESPA cutoff if (cut_respa && cut[i][j] < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -559,22 +553,21 @@ double PairLJClass2::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut[i][j] * cut[i][j] * cut[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut[i][j]; @@ -588,14 +581,14 @@ void PairLJClass2::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -609,21 +602,21 @@ void PairLJClass2::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -634,10 +627,10 @@ void PairLJClass2::read_restart(FILE *fp) void PairLJClass2::write_restart_settings(FILE *fp) { - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -648,26 +641,24 @@ void PairLJClass2::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } - /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairLJClass2::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -678,27 +669,25 @@ void PairLJClass2::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut[i][j]); } /* ---------------------------------------------------------------------- */ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, - double /*factor_coul*/, double factor_lj, - double &fforce) + double /*factor_coul*/, double factor_lj, double &fforce) { - double r2inv,rinv,r3inv,r6inv,forcelj,philj; + double r2inv, rinv, r3inv, r6inv, forcelj, philj; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fforce = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fforce = factor_lj * forcelj * r2inv; - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - return factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + return factor_lj * philj; } /* ---------------------------------------------------------------------- */ @@ -706,7 +695,7 @@ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double r void *PairLJClass2::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index c4ac5ae7c8..844c979994 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,17 +13,17 @@ #include "pair_lj_class2_coul_cut.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -65,14 +64,14 @@ PairLJClass2CoulCut::~PairLJClass2CoulCut() void PairLJClass2CoulCut::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj; - double factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj; + double factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -110,47 +109,49 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - else forcecoul = 0.0; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + fpair = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq[itype][jtype]) - ecoul = factor_coul * qqrd2e * qtmp*q[j]*sqrt(r2inv); - else ecoul = 0.0; + ecoul = factor_coul * qqrd2e * qtmp * q[j] * sqrt(r2inv); + else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -165,26 +166,25 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) void PairLJClass2CoulCut::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(cut_coul,n+1,n+1,"pair:cut_coul"); - memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(cut_coul, np1, np1, "pair:cut_coul"); + memory->create(cut_coulsq, np1, np1, "pair:cut_coulsq"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -193,16 +193,18 @@ void PairLJClass2CoulCut::allocate() void PairLJClass2CoulCut::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul_global = cut_lj_global; - else cut_coul_global = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul_global = cut_lj_global; + else + cut_coul_global = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) { @@ -218,26 +220,25 @@ void PairLJClass2CoulCut::settings(int narg, char **arg) void PairLJClass2CoulCut::coeff(int narg, char **arg) { - if (narg < 4 || narg > 6) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 6) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; - if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp); - if (narg == 6) cut_coul_one = utils::numeric(FLERR,arg[5],false,lmp); + if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR, arg[4], false, lmp); + if (narg == 6) cut_coul_one = utils::numeric(FLERR, arg[5], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; @@ -247,7 +248,7 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -256,10 +257,9 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg) void PairLJClass2CoulCut::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style lj/class2/coul/cut requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/cut requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); } /* ---------------------------------------------------------------------- @@ -272,28 +272,27 @@ double PairLJClass2CoulCut::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); - cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]); + cut_coul[i][j] = mix_distance(cut_coul[i][i], cut_coul[j][j]); } - double cut = MAX(cut_lj[i][j],cut_coul[i][j]); + double cut = MAX(cut_lj[i][j], cut_coul[i][j]); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; cut_coulsq[j][i] = cut_coulsq[i][j]; @@ -310,22 +309,21 @@ double PairLJClass2CoulCut::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut; @@ -339,15 +337,15 @@ void PairLJClass2CoulCut::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); - fwrite(&cut_coul[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); + fwrite(&cut_coul[i][j], sizeof(double), 1, fp); } } } @@ -361,23 +359,23 @@ void PairLJClass2CoulCut::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -388,11 +386,11 @@ void PairLJClass2CoulCut::read_restart(FILE *fp) void PairLJClass2CoulCut::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -402,17 +400,17 @@ void PairLJClass2CoulCut::write_restart_settings(FILE *fp) void PairLJClass2CoulCut::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- @@ -421,8 +419,7 @@ void PairLJClass2CoulCut::read_restart_settings(FILE *fp) void PairLJClass2CoulCut::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -433,39 +430,38 @@ void PairLJClass2CoulCut::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { - double r2inv,rinv,r3inv,r6inv,forcecoul,forcelj,phicoul,philj; + double r2inv, rinv, r3inv, r6inv, forcecoul, forcelj, phicoul, philj; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - else forcecoul = 0.0; + forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; - fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; + fforce = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; if (rsq < cut_coulsq[itype][jtype]) { - phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - eng += factor_coul*phicoul; + phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + eng += factor_coul * phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + eng += factor_lj * philj; } return eng; @@ -476,9 +472,8 @@ double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, void *PairLJClass2CoulCut::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } - diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index 63265509c9..79773d6408 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,32 +13,32 @@ #include "pair_lj_class2_coul_long.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "kspace.h" -#include "update.h" -#include "respa.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "respa.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 +static constexpr double EWALD_F = 1.12837917; +static constexpr double EWALD_P = 0.3275911; +static constexpr double A1 = 0.254829592; +static constexpr double A2 = -0.284496736; +static constexpr double A3 = 1.421413741; +static constexpr double A4 = -1.453152027; +static constexpr double A5 = 1.061405429; /* ---------------------------------------------------------------------- */ @@ -79,16 +78,16 @@ PairLJClass2CoulLong::~PairLJClass2CoulLong() void PairLJClass2CoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itable,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj; - double grij,expm2,prefactor,t,erfc; - double factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itable, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double rsq, r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj; + double grij, expm2, prefactor, t, erfc; + double factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -126,75 +125,77 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - ecoul = prefactor*erfc; + ecoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; } - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -206,11 +207,11 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) void PairLJClass2CoulLong::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -231,8 +232,8 @@ void PairLJClass2CoulLong::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -255,34 +256,35 @@ void PairLJClass2CoulLong::compute_inner() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -293,11 +295,11 @@ void PairLJClass2CoulLong::compute_inner() void PairLJClass2CoulLong::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -321,10 +323,10 @@ void PairLJClass2CoulLong::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -347,38 +349,39 @@ void PairLJClass2CoulLong::compute_middle() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -389,17 +392,17 @@ void PairLJClass2CoulLong::compute_middle() void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; + int i, j, ii, jj, inum, jnum, itype, jtype, itable; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; + double grij, expm2, prefactor, t, erfc; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double rsq; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -420,8 +423,8 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -444,32 +447,30 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0); + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2 - 1.0); if (rsq > cut_in_off_sq) { if (rsq < cut_in_on_sq) { - rsw = (r - cut_in_off)/cut_in_diff; - forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + rsw = (r - cut_in_off) / cut_in_diff; + forcecoul += prefactor * rsw * rsw * (3.0 - 2.0 * rsw); if (factor_coul < 1.0) - forcecoul -= - (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + forcecoul -= (1.0 - factor_coul) * prefactor * rsw * rsw * (3.0 - 2.0 * rsw); } else { forcecoul += prefactor; - if (factor_coul < 1.0) - forcecoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } } } else { @@ -478,96 +479,99 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + forcelj *= rsw * rsw * (3.0 - 2.0 * rsw); } - } else forcelj = 0.0; + } else + forcelj = 0.0; fpair = (forcecoul + forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - ecoul = prefactor*erfc; - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + ecoul = prefactor * erfc; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - ecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + ecoul -= (1.0 - factor_coul) * prefactor; } } - } else ecoul = 0.0; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } if (vflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { - table = vtable[itable] + fraction*dvtable[itable]; - forcecoul = qtmp*q[j] * table; + table = vtable[itable] + fraction * dvtable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq <= cut_in_off_sq) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } else if (rsq <= cut_in_on_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -580,24 +584,22 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) void PairLJClass2CoulLong::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -606,16 +608,18 @@ void PairLJClass2CoulLong::allocate() void PairLJClass2CoulLong::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul = cut_lj_global; - else cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul = cut_lj_global; + else + cut_coul = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; @@ -628,24 +632,23 @@ void PairLJClass2CoulLong::settings(int narg, char **arg) void PairLJClass2CoulLong::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); - double cut_lj_one = cut_lj_global; - if (narg == 5) cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp); + double cut_lj_one = cut_lj_global; + if (narg == 5) cut_lj_one = utils::numeric(FLERR, arg[4], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; @@ -654,7 +657,7 @@ void PairLJClass2CoulLong::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -663,25 +666,23 @@ void PairLJClass2CoulLong::coeff(int narg, char **arg) void PairLJClass2CoulLong::init_style() { - if (!atom->q_flag) - error->all(FLERR, - "Pair style lj/class2/coul/long requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/long requires atom attribute q"); // request regular or rRESPA neighbor list int irequest; int respa = 0; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; } - irequest = neighbor->request(this,instance_me); + irequest = neighbor->request(this, instance_me); if (respa >= 1) { neighbor->requests[irequest]->respaouter = 1; - neighbor->requests[irequest]->respainner = 1; + neighbor->requests[irequest]->respainner = 1; } if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; @@ -689,19 +690,19 @@ void PairLJClass2CoulLong::init_style() // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; // insure use of KSpace long-range solver, set g_ewald - if (force->kspace == nullptr) - error->all(FLERR,"Pair style requires a KSpace style"); + if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables - if (ncoultablebits) init_tables(cut_coul,cut_respa); + if (ncoultablebits) init_tables(cut_coul, cut_respa); } /* ---------------------------------------------------------------------- @@ -714,26 +715,25 @@ double PairLJClass2CoulLong::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]); } - double cut = MAX(cut_lj[i][j],cut_coul); + double cut = MAX(cut_lj[i][j], cut_coul); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; lj1[j][i] = lj1[i][j]; @@ -744,8 +744,8 @@ double PairLJClass2CoulLong::init_one(int i, int j) // check interior rRESPA cutoff - if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + if (cut_respa && MIN(cut_lj[i][j], cut_coul) < cut_respa[3]) + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -754,22 +754,21 @@ double PairLJClass2CoulLong::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut; @@ -783,14 +782,14 @@ void PairLJClass2CoulLong::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); } } } @@ -804,21 +803,21 @@ void PairLJClass2CoulLong::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -829,13 +828,13 @@ void PairLJClass2CoulLong::read_restart(FILE *fp) void PairLJClass2CoulLong::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); - fwrite(&ncoultablebits,sizeof(int),1,fp); - fwrite(&tabinner,sizeof(double),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); + fwrite(&ncoultablebits, sizeof(int), 1, fp); + fwrite(&tabinner, sizeof(double), 1, fp); } /* ---------------------------------------------------------------------- @@ -845,21 +844,21 @@ void PairLJClass2CoulLong::write_restart_settings(FILE *fp) void PairLJClass2CoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); - MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); - MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world); + MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- @@ -868,8 +867,7 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp) void PairLJClass2CoulLong::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -880,69 +878,68 @@ void PairLJClass2CoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { - double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor; - double fraction,table,forcecoul,forcelj,phicoul,philj; + double r2inv, r, rinv, r3inv, r6inv, grij, expm2, t, erfc, prefactor; + double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = atom->q[i]*atom->q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = atom->q[i] * atom->q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = atom->q[i]*atom->q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = atom->q[i] * atom->q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; - fforce = (forcecoul + factor_lj*forcelj) * r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; + fforce = (forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor*erfc; + phicoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - phicoul = atom->q[i]*atom->q[j] * table; + table = etable[itable] + fraction * detable[itable]; + phicoul = atom->q[i] * atom->q[j] * table; } - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + eng += factor_lj * philj; } return eng; @@ -953,9 +950,9 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, void *PairLJClass2CoulLong::extract(const char *str, int &dim) { dim = 0; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } diff --git a/src/DPD-REACT/pair_multi_lucy.cpp b/src/DPD-REACT/pair_multi_lucy.cpp index 0f54a57870..41c9d9fb66 100644 --- a/src/DPD-REACT/pair_multi_lucy.cpp +++ b/src/DPD-REACT/pair_multi_lucy.cpp @@ -37,6 +37,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -104,7 +105,6 @@ void PairMultiLucy::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - double pi = MathConst::MY_PI; double A_i; double A_j; double fraction_i,fraction_j; @@ -198,7 +198,7 @@ void PairMultiLucy::compute(int eflag, int vflag) evdwl = tb->e[itable] + fraction_i*tb->de[itable]; } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy"); - evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0; + evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0; if (evflag) ev_tally(0,0,nlocal,newton_pair, evdwl,0.0,0.0,0.0,0.0,0.0); @@ -733,7 +733,6 @@ void PairMultiLucy::computeLocalDensity() numneigh = list->numneigh; firstneigh = list->firstneigh; - double pi = MathConst::MY_PI; double *rho = atom->rho; // zero out density @@ -766,7 +765,7 @@ void PairMultiLucy::computeLocalDensity() if (rsq < cutsq[itype][jtype]) { double rcut = sqrt(cutsq[itype][jtype]); - factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut); + factor= (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut); rho[i] += factor; if (newton_pair || j < nlocal) { rho[j] += factor; diff --git a/src/DPD-REACT/pair_multi_lucy_rx.cpp b/src/DPD-REACT/pair_multi_lucy_rx.cpp index 64b2b4c7d5..8b348810fd 100644 --- a/src/DPD-REACT/pair_multi_lucy_rx.cpp +++ b/src/DPD-REACT/pair_multi_lucy_rx.cpp @@ -39,6 +39,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -127,7 +128,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag) double *uCG = atom->uCG; double *uCGnew = atom->uCGnew; - double pi = MathConst::MY_PI; double A_i, A_j; double fraction_i,fraction_j; int jtable; @@ -276,7 +276,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag) evdwl = tb->e[itable] + fraction_i*tb->de[itable]; } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx"); - evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0; + evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0; evdwlOld = mixWtSite1old_i*evdwl; evdwl = mixWtSite1_i*evdwl; @@ -866,15 +866,13 @@ void PairMultiLucyRX::computeLocalDensity() const int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - const double pi = MathConst::MY_PI; - const bool newton_pair = force->newton_pair; const bool one_type = (atom->ntypes == 1); // Special cut-off values for when there's only one type. const double cutsq_type11 = cutsq[1][1]; const double rcut_type11 = sqrt(cutsq_type11); - const double factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11); + const double factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11); double *rho = atom->rho; @@ -925,7 +923,7 @@ void PairMultiLucyRX::computeLocalDensity() const double rcut = sqrt(cutsq[itype][jtype]); const double tmpFactor = 1.0-sqrt(rsq)/rcut; const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; - const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; + const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; rho_i += factor; if (newton_pair || j < nlocal) rho[j] += factor; diff --git a/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp b/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp index 8c558db32c..c2942288ff 100644 --- a/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp +++ b/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp @@ -23,20 +23,24 @@ ------------------------------------------------------------------------------------------- */ #include "pair_multi_lucy_rx_kokkos.h" + +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "kokkos.h" +#include "math_const.h" +#include "math_const.h" +#include "memory_kokkos.h" +#include "neigh_list.h" +#include "neigh_request.h" + #include #include -#include "math_const.h" -#include "atom_kokkos.h" -#include "force.h" -#include "comm.h" -#include "neigh_list.h" -#include "memory_kokkos.h" -#include "error.h" -#include "atom_masks.h" -#include "neigh_request.h" -#include "kokkos.h" using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -282,7 +286,6 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXCompute::operator()(TagPairMultiLucyRXCompute()() = 3; - evdwl *=(pi*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0; + evdwl *=(MY_PI*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0; evdwlOld = mixWtSite1old_i*evdwl; evdwl = mixWtSite1_i*evdwl; @@ -458,15 +461,13 @@ void PairMultiLucyRXKokkos::computeLocalDensity() d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; - const double pi = MathConst::MY_PI; - const bool newton_pair = force->newton_pair; const bool one_type = (atom->ntypes == 1); // Special cut-off values for when there's only one type. cutsq_type11 = cutsq[1][1]; rcut_type11 = sqrt(cutsq_type11); - factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11); + factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11); // zero out density int m = nlocal; @@ -548,8 +549,6 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXComputeLoca const int itype = type[i]; const int jnum = d_numneigh[i]; - const double pi = MathConst::MY_PI; - for (int jj = 0; jj < jnum; jj++) { const int j = (d_neighbors(i,jj) & NEIGHMASK); const int jtype = type[j]; @@ -574,7 +573,7 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXComputeLoca const double rcut = sqrt(d_cutsq(itype,jtype)); const double tmpFactor = 1.0-sqrt(rsq)/rcut; const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; - const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; + const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; rho_i_contrib += factor; if ((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR || j < nlocal)) a_rho[j] += factor; diff --git a/src/LATTE/fix_latte.cpp b/src/LATTE/fix_latte.cpp index e1d0224e62..63b24b1cdd 100644 --- a/src/LATTE/fix_latte.cpp +++ b/src/LATTE/fix_latte.cpp @@ -17,20 +17,18 @@ ------------------------------------------------------------------------- */ #include "fix_latte.h" -#include -#include + #include "atom.h" #include "comm.h" -#include "update.h" -#include "neighbor.h" -#include "domain.h" -#include "force.h" -#include "neigh_request.h" -#include "neigh_list.h" -#include "modify.h" #include "compute.h" -#include "memory.h" +#include "domain.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -82,13 +80,12 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[3],"NULL") != 0) { coulomb = 1; - error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation " - "of a Coulomb potential"); + error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation of a Coulomb potential"); id_pe = utils::strdup(arg[3]); - int ipe = modify->find_compute(id_pe); - if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID"); - if (modify->compute[ipe]->peatomflag == 0) + c_pe = modify->get_compute_by_id(id_pe); + if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); + if (c_pe->peatomflag == 0) error->all(FLERR,"Fix latte compute ID does not compute pe/atom"); } @@ -105,7 +102,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : FixLatte::~FixLatte() { - delete [] id_pe; + delete[] id_pe; memory->destroy(qpotential); memory->destroy(flatte); } @@ -134,9 +131,8 @@ void FixLatte::init() if (atom->q_flag == 0 || force->pair == nullptr || force->kspace == nullptr) error->all(FLERR,"Fix latte cannot compute Coulomb potential"); - int ipe = modify->find_compute(id_pe); - if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID"); - c_pe = modify->compute[ipe]; + c_pe = modify->get_compute_by_id(id_pe); + if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); } // must be fully periodic or fully non-periodic @@ -153,33 +149,6 @@ void FixLatte::init() memory->create(qpotential,atom->nlocal,"latte:qpotential"); memory->create(flatte,atom->nlocal,3,"latte:flatte"); } - - /* - // warn if any integrate fix comes after this one - // is it actually necessary for q(n) update to come after x,v update ?? - - int after = 0; - int flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(id,modify->fix[i]->id) == 0) after = 1; - else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1; - } - if (flag && comm->me == 0) - error->warning(FLERR,"Fix latte should come after all other " - "integration fixes"); - */ - - /* - // need a full neighbor list - // could we use a half list? - // perpetual list, built whenever re-neighboring occurs - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->fix = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; - */ } /* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index bac65a2464..239b125f93 100644 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -292,7 +292,7 @@ void PairSW::read_file(char *file) if (comm->me == 0) { PotentialFileReader reader(lmp, file, "sw", unit_convert_flag); - char * line; + char *line; // transparently convert units for supported conversions diff --git a/src/MANYBODY/pair_tersoff_table.cpp b/src/MANYBODY/pair_tersoff_table.cpp index 0c11760737..12733e1c60 100644 --- a/src/MANYBODY/pair_tersoff_table.cpp +++ b/src/MANYBODY/pair_tersoff_table.cpp @@ -27,6 +27,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "math_const.h" #include "memory.h" #include "neigh_list.h" #include "neigh_request.h" @@ -37,6 +38,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; #define MAXLINE 1024 #define DELTA 4 @@ -543,7 +545,6 @@ void PairTersoffTable::allocateGrids() double deltaArgumentCutoffFunction, deltaArgumentExponential, deltaArgumentBetaZetaPower; double deltaArgumentGtetaFunction; double r, minMu, maxLambda, maxCutoff; - double const PI=acos(-1.0); deallocateGrids(); @@ -652,8 +653,8 @@ void PairTersoffTable::allocateGrids() } for (l = numGridPointsOneCutoffFunction; l < numGridPointsCutoffFunction; l++) { - cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (PI * (r - cutoffR)/(cutoffS-cutoffR)) ; - cutoffFunctionDerived[i][j][l] = -0.5 * PI * sin (PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR) ; + cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) ; + cutoffFunctionDerived[i][j][l] = -0.5 * MY_PI * sin (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR); r += deltaArgumentCutoffFunction; } } diff --git a/src/MEAM/meam_setup_param.cpp b/src/MEAM/meam_setup_param.cpp index b15b70e263..5612936c20 100644 --- a/src/MEAM/meam_setup_param.cpp +++ b/src/MEAM/meam_setup_param.cpp @@ -233,7 +233,6 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3 // 21 = theta // see alloyparams(void) in meam_setup_done.cpp case 21: - // const double PI = 3.141592653589793238463; meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; diff --git a/src/OPENMP/fix_rigid_nh_omp.cpp b/src/OPENMP/fix_rigid_nh_omp.cpp index d95585e0ac..4e73eaa89b 100644 --- a/src/OPENMP/fix_rigid_nh_omp.cpp +++ b/src/OPENMP/fix_rigid_nh_omp.cpp @@ -246,12 +246,11 @@ void FixRigidNHOMP::compute_forces_and_torques() if (rstyle == SINGLE) { // we have just one rigid body. use OpenMP reduction to get sum[] double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; - int i; #if defined(_OPENMP) -#pragma omp parallel for LMP_DEFAULT_NONE private(i) reduction(+:s0,s1,s2,s3,s4,s5) +#pragma omp parallel for LMP_DEFAULT_NONE reduction(+:s0,s1,s2,s3,s4,s5) #endif - for (i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if (ibody < 0) continue; @@ -285,12 +284,11 @@ void FixRigidNHOMP::compute_forces_and_torques() for (int ib = 0; ib < nbody; ++ib) { double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; - int i; #if defined(_OPENMP) -#pragma omp parallel for LMP_DEFAULT_NONE private(i) LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5) +#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5) #endif - for (i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if (ibody != ib) continue; diff --git a/src/OPENMP/pair_comb_omp.cpp b/src/OPENMP/pair_comb_omp.cpp index 116db7194d..b3306d678c 100644 --- a/src/OPENMP/pair_comb_omp.cpp +++ b/src/OPENMP/pair_comb_omp.cpp @@ -384,7 +384,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) double PairCombOMP::yasu_char(double *qf_fix, int &igroup) { - int ii; double potal,fac11,fac11e; const double * const * const x = atom->x; @@ -401,7 +400,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) const int groupbit = group->bitmask[igroup]; qf = qf_fix; - for (ii = 0; ii < inum; ii++) { + for (int ii = 0; ii < inum; ii++) { const int i = ilist[ii]; if (mask[i] & groupbit) qf[i] = 0.0; @@ -417,9 +416,9 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) // loop over full neighbor list of my atoms #if defined(_OPENMP) -#pragma omp parallel for private(ii) LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e) +#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e) #endif - for (ii = 0; ii < inum; ii ++) { + for (int ii = 0; ii < inum; ii ++) { double fqi,fqij,fqji,fqjj,delr1[3]; double sr1,sr2,sr3; int mr1,mr2,mr3; @@ -533,7 +532,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) // sum charge force on each node and return it double eneg = 0.0; - for (ii = 0; ii < inum; ii++) { + for (int ii = 0; ii < inum; ii++) { const int i = ilist[ii]; if (mask[i] & groupbit) eneg += qf[i]; diff --git a/src/OPENMP/pair_reaxff_omp.cpp b/src/OPENMP/pair_reaxff_omp.cpp index 1f320f76fc..a329809c4e 100644 --- a/src/OPENMP/pair_reaxff_omp.cpp +++ b/src/OPENMP/pair_reaxff_omp.cpp @@ -402,8 +402,7 @@ int PairReaxFFOMP::estimate_reax_lists() int PairReaxFFOMP::write_reax_lists() { - int itr_i, itr_j, i, j, num_mynbrs; - int *jlist; + int num_mynbrs; double d_sqr, dist, cutoff_sqr; rvec dvec; @@ -425,19 +424,19 @@ int PairReaxFFOMP::write_reax_lists() num_nbrs = 0; - for (itr_i = 0; itr_i < numall; ++itr_i) { - i = ilist[itr_i]; + for (int itr_i = 0; itr_i < numall; ++itr_i) { + int i = ilist[itr_i]; num_nbrs_offset[i] = num_nbrs; num_nbrs += numneigh[i]; } #if defined(_OPENMP) #pragma omp parallel for schedule(dynamic,50) default(shared) \ - private(itr_i, itr_j, i, j, jlist, cutoff_sqr, num_mynbrs, d_sqr, dvec, dist) + private(cutoff_sqr, num_mynbrs, d_sqr, dvec, dist) #endif - for (itr_i = 0; itr_i < numall; ++itr_i) { - i = ilist[itr_i]; - jlist = firstneigh[i]; + for (int itr_i = 0; itr_i < numall; ++itr_i) { + int i = ilist[itr_i]; + auto jlist = firstneigh[i]; Set_Start_Index(i, num_nbrs_offset[i], far_nbrs); if (i < inum) @@ -447,8 +446,8 @@ int PairReaxFFOMP::write_reax_lists() num_mynbrs = 0; - for (itr_j = 0; itr_j < numneigh[i]; ++itr_j) { - j = jlist[itr_j]; + for (int itr_j = 0; itr_j < numneigh[i]; ++itr_j) { + int j = jlist[itr_j]; j &= NEIGHMASK; get_distance(x[j], x[i], &d_sqr, &dvec); diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index 105b24f83b..8b773a30ea 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -1192,13 +1192,7 @@ void FixRigidNHSmall::compute_dof() for (int k = 0; k < dimension; k++) if (fabs(b->inertia[k]) < EPSILON) nf_r--; } - } else if (dimension == 2) { - nf_r = nlocal_body; - for (int ibody = 0; ibody < nlocal_body; ibody++) { - Body *b = &body[ibody]; - if (fabs(b->inertia[2]) < EPSILON) nf_r--; - } - } + } else if (dimension == 2) nf_r = nlocal_body; double nf[2], nfall[2]; nf[0] = nf_t; diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 49116cb8cc..f012789414 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -51,12 +51,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : type_flag(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), - shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), - list(nullptr), b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), - b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), - a_count(nullptr), a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), - a_ave_all(nullptr), a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), - onemols(nullptr) + shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), + b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), b_max(nullptr), + b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), + a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr), + a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); @@ -172,8 +171,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : if (imol == -1) error->all(FLERR,"Molecule template ID for fix shake does not exist"); if (atom->molecules[imol]->nset > 1 && comm->me == 0) - error->warning(FLERR,"Molecule template for " - "fix shake has multiple molecules"); + error->warning(FLERR,"Molecule template for fix shake has multiple molecules"); onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; iarg += 2; @@ -199,6 +197,8 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : int nb = atom->nbondtypes + 1; b_count = new int[nb]; b_count_all = new int[nb]; + b_atom = new int[nb]; + b_atom_all = new int[nb]; b_ave = new double[nb]; b_ave_all = new double[nb]; b_max = new double[nb]; @@ -228,8 +228,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : find_clusters(); if (comm->me == 0) - utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n", - platform::walltime()-time1); + utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n",platform::walltime()-time1); // initialize list of SHAKE clusters to constrain @@ -281,32 +280,34 @@ FixShake::~FixShake() memory->destroy(vtmp); - delete [] bond_flag; - delete [] angle_flag; - delete [] type_flag; - delete [] mass_list; + delete[] bond_flag; + delete[] angle_flag; + delete[] type_flag; + delete[] mass_list; - delete [] bond_distance; - delete [] angle_distance; + delete[] bond_distance; + delete[] angle_distance; if (output_every) { - delete [] b_count; - delete [] b_count_all; - delete [] b_ave; - delete [] b_ave_all; - delete [] b_max; - delete [] b_max_all; - delete [] b_min; - delete [] b_min_all; + delete[] b_count; + delete[] b_count_all; + delete[] b_atom; + delete[] b_atom_all; + delete[] b_ave; + delete[] b_ave_all; + delete[] b_max; + delete[] b_max_all; + delete[] b_min; + delete[] b_min_all; - delete [] a_count; - delete [] a_count_all; - delete [] a_ave; - delete [] a_ave_all; - delete [] a_max; - delete [] a_max_all; - delete [] a_min; - delete [] a_min_all; + delete[] a_count; + delete[] a_count_all; + delete[] a_ave; + delete[] a_ave_all; + delete[] a_max; + delete[] a_max_all; + delete[] a_min; + delete[] a_min_all; } memory->destroy(list); @@ -2473,6 +2474,7 @@ void FixShake::stats() b_count[i] = 0; b_ave[i] = b_max[i] = 0.0; b_min[i] = BIG; + b_atom[i] = -1; } for (i = 0; i < na; i++) { a_count[i] = 0; @@ -2504,6 +2506,7 @@ void FixShake::stats() m = shake_type[i][j-1]; b_count[m]++; + b_atom[m] = n; b_ave[m] += r; b_max[m] = MAX(b_max[m],r); b_min[m] = MIN(b_min[m],r); @@ -2547,6 +2550,7 @@ void FixShake::stats() // sum across all procs MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world); + MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world); MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world); @@ -2563,16 +2567,16 @@ void FixShake::stats() auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n", update->ntimestep); for (i = 1; i < nb; i++) { - const auto bcnt = b_count_all[i]/2; + const auto bcnt = b_count_all[i]; if (bcnt) mesg += fmt::format("Bond: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, - b_ave_all[i]/bcnt/2.0,b_max_all[i]-b_min_all[i],bcnt); + b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]); } for (i = 1; i < na; i++) { - const auto acnt = a_count_all[i]/3; + const auto acnt = a_count_all[i]; if (acnt) mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, - a_ave_all[i]/acnt/3.0,a_max_all[i]-a_min_all[i],acnt); + a_ave_all[i]/acnt,a_max_all[i]-a_min_all[i],acnt/3); } utils::logmesg(lmp,mesg); } diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 5a956bd798..677cdfa942 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -108,10 +108,10 @@ class FixShake : public Fix { int nlist, maxlist; // size and max-size of list // stat quantities - int *b_count, *b_count_all; // counts for each bond type - double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type - double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays - int *a_count, *a_count_all; // ditto for angle types + int *b_count, *b_count_all, *b_atom, *b_atom_all; // counts for each bond type, atoms in bond cluster + double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type + double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays + int *a_count, *a_count_all; // ditto for angle types double *a_ave, *a_max, *a_min; double *a_ave_all, *a_max_all, *a_min_all; diff --git a/src/thermo.cpp b/src/thermo.cpp index 27d74c58b6..f95b665e4e 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -905,7 +905,8 @@ void Thermo::parse_fields(char *str) if (argi.get_type() == ArgInfo::COMPUTE) { auto icompute = modify->get_compute_by_id(argi.get_name()); - if (!icompute) error->all(FLERR,"Could not find thermo custom compute ID"); + if (!icompute) + error->all(FLERR,"Could not find thermo custom compute ID: {}", argi.get_name()); if (argindex1[nfield] == 0 && icompute->scalar_flag == 0) error->all(FLERR,"Thermo compute does not compute scalar"); if (argindex1[nfield] > 0 && argindex2[nfield] == 0) { @@ -935,7 +936,7 @@ void Thermo::parse_fields(char *str) } else if (argi.get_type() == ArgInfo::FIX) { auto ifix = modify->get_fix_by_id(argi.get_name()); - if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID"); + if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID: {}", argi.get_name()); if (argindex1[nfield] == 0 && ifix->scalar_flag == 0) error->all(FLERR,"Thermo fix does not compute scalar"); if (argindex1[nfield] > 0 && argindex2[nfield] == 0) { @@ -961,7 +962,7 @@ void Thermo::parse_fields(char *str) } else if (argi.get_type() == ArgInfo::VARIABLE) { int n = input->variable->find(argi.get_name()); if (n < 0) - error->all(FLERR,"Could not find thermo custom variable name"); + error->all(FLERR,"Could not find thermo custom variable name: {}", argi.get_name()); if (argindex1[nfield] == 0 && input->variable->equalstyle(n) == 0) error->all(FLERR,"Thermo custom variable is not equal-style variable"); if (argindex1[nfield] && input->variable->vectorstyle(n) == 0) diff --git a/unittest/utils/test_platform.cpp b/unittest/utils/test_platform.cpp index 37e749b9be..a809b61458 100644 --- a/unittest/utils/test_platform.cpp +++ b/unittest/utils/test_platform.cpp @@ -21,7 +21,7 @@ TEST(Platform, clock) // spend some time computing pi constexpr double known_pi = 3.141592653589793238462643; - constexpr int n = 10000000; + constexpr int n = 100000000; constexpr double h = 1.0 / (double)n; double my_pi = 0.0, x; for (int i = 0; i < n; ++i) { @@ -34,7 +34,8 @@ TEST(Platform, clock) ASSERT_NEAR(my_pi, known_pi, 1e-12); ASSERT_GT(wt_used, 1e-4); - ASSERT_GT(ct_used, 1e-4); + // windows sometimes incorrectly reports used CPU time as 0.0 + if (ct_used != 0.0) ASSERT_GT(ct_used, 1e-4); } TEST(Platform, putenv_unsetenv)