diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index 6847b8c030..4dd17d71f8 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -352,7 +352,7 @@ KIM Extra unit tests (CMake only) During development, testing, or debugging, if :doc:`unit testing ` is enabled in LAMMPS, one can also enable extra tests on :doc:`KIM commands ` by setting the -``KIM_EXTRA_UNITTESTS`` to *yes* (or *on*). +``KIM_EXTRA_UNITTESTS`` to *yes* (or *on*). Enabling the extra unit tests have some requirements, @@ -367,9 +367,9 @@ Enabling the extra unit tests have some requirements, *conda-forge* channel as ``conda install kim-property`` if LAMMPS is built in Conda. More detailed information is available at: `kim-property installation `_. -* It is also necessary to install - ``EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000``, and - ``EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005`` KIM models. +* It is also necessary to install + ``EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000``, and + ``EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005`` KIM models. See `Obtaining KIM Models `_ to learn how to install a pre-build binary of the OpenKIM Repository of Models or see diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 9a2ae60f1e..864c544fed 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -84,13 +84,15 @@ information is available, then also a heuristic based on that bond length is computed. It is used as communication cutoff, if there is no pair style present and no *comm_modify cutoff* command used. Otherwise a warning is printed, if this bond based estimate is larger than the -communication cutoff used. A +communication cutoff used. The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to communication mode *multi* instead. Since in this case the communication cutoffs are determined per atom type, a type specifier is needed and cutoff for one or multiple types can be extended. Also ranges of types -using the usual asterisk notation can be given. +using the usual asterisk notation can be given. For granular pair styles, +the default cutoff is set to the sum of the current maximum atomic radii +for each type. These are simulation scenarios in which it may be useful or even necessary to set a ghost cutoff > neighbor cutoff: diff --git a/doc/src/fix_temp_csvr.rst b/doc/src/fix_temp_csvr.rst index 5e6074bc66..440a264ed5 100644 --- a/doc/src/fix_temp_csvr.rst +++ b/doc/src/fix_temp_csvr.rst @@ -124,16 +124,16 @@ temperature is calculated taking the bias into account, bias is removed from each atom, thermostatting is performed on the remaining thermal degrees of freedom, and the bias is added back in. -An important feature of these thermostats is that they have an +An important feature of these thermostats is that they have an associated effective energy that is a constant of motion. -The effective energy is the total energy (kinetic + potential) plus -the accumulated kinetic energy changes due to the thermostat. The -latter quantity is the global scalar computed by these fixes. This -feature is useful to check the integration of the equations of motion -against discretization errors. In other words, the conservation of -the effective energy can be used to choose an appropriate integration -:doc:`timestep `. This is similar to the usual paradigm of -checking the conservation of the total energy in the microcanonical +The effective energy is the total energy (kinetic + potential) plus +the accumulated kinetic energy changes due to the thermostat. The +latter quantity is the global scalar computed by these fixes. This +feature is useful to check the integration of the equations of motion +against discretization errors. In other words, the conservation of +the effective energy can be used to choose an appropriate integration +:doc:`timestep `. This is similar to the usual paradigm of +checking the conservation of the total energy in the microcanonical ensemble. diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 1bb591587c..98ee0d0b6a 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -49,7 +49,9 @@ sometimes be faster. Either style should give the same answers. The *multi* style is a modified binning algorithm that is useful for systems with a wide range of cutoff distances, e.g. due to different -size particles. For the *bin* style, the bin size is set to 1/2 of +size particles. For granular pair styles, cutoffs are set to the +sum of the maximum atomic radii for each atom type. +For the *bin* style, the bin size is set to 1/2 of the largest cutoff distance between any pair of atom types and a single set of bins is defined to search over for all atom types. This can be inefficient if one pair of types has a very long cutoff, but @@ -57,8 +59,10 @@ other type pairs have a much shorter cutoff. For style *multi* the bin size is set to 1/2 of the shortest cutoff distance and multiple sets of bins are defined to search over for different atom types. This imposes some extra setup overhead, but the searches themselves -may be much faster for the short-cutoff cases. See the :doc:`comm_modify mode multi ` command for a communication option -that may also be beneficial for simulations of this kind. +may be much faster for the short-cutoff cases. +See the :doc:`comm_modify mode multi ` command for a +communication option that may also be beneficial for simulations of +this kind. The :doc:`neigh_modify ` command has additional options that control how often neighbor lists are built and which pairs are diff --git a/python/lammps/core.py b/python/lammps/core.py index 161583b78c..c0cbaac533 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -77,7 +77,7 @@ class lammps(object): modpath = dirname(abspath(getsourcefile(lambda:0))) # for windows installers the shared library is in a different folder - winpath = abspath(os.path.join(modpath,'..','bin')) + winpath = abspath(os.path.join(modpath,'..','..','bin')) self.lib = None self.lmp = None diff --git a/src/GPU/pair_lj_cubic_gpu.cpp b/src/GPU/pair_lj_cubic_gpu.cpp index a669d52a19..35062a5d71 100644 --- a/src/GPU/pair_lj_cubic_gpu.cpp +++ b/src/GPU/pair_lj_cubic_gpu.cpp @@ -16,10 +16,11 @@ ------------------------------------------------------------------------- */ #include "pair_lj_cubic_gpu.h" + #include #include - #include + #include "atom.h" #include "atom_vec.h" #include "comm.h" @@ -36,6 +37,8 @@ #include "gpu_extra.h" #include "suffix.h" +#include "pair_lj_cubic_const.h" + using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp index ad671de5b3..ee6f5f98a5 100644 --- a/src/KIM/pair_kim.cpp +++ b/src/KIM/pair_kim.cpp @@ -941,7 +941,7 @@ void PairKIM::set_argument_pointers() // Set KIM pointer appropriately for particleVirial if (KIM_SupportStatus_Equal(kim_model_support_for_particleVirial, - KIM_SUPPORT_STATUS_required) + KIM_SUPPORT_STATUS_required) && (vflag_atom == 0)) { // reallocate per-atom virial array if necessary if (atom->nmax > maxvatom) { diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp index 1323e589d3..3703769e90 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp @@ -21,6 +21,7 @@ #include "atom_masks.h" #include "kokkos.h" #include "math_const.h" +#include "math_special.h" #include "memory_kokkos.h" #include "neigh_list.h" #include "neigh_request.h" @@ -32,6 +33,7 @@ using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; using namespace std; #ifdef DBL_EPSILON diff --git a/src/KOKKOS/pair_zbl_kokkos.cpp b/src/KOKKOS/pair_zbl_kokkos.cpp index a3cdbf3b6e..780ae19cc8 100644 --- a/src/KOKKOS/pair_zbl_kokkos.cpp +++ b/src/KOKKOS/pair_zbl_kokkos.cpp @@ -30,6 +30,8 @@ #include "atom_masks.h" #include "kokkos.h" +#include "pair_zbl_const.h" + // From J.F. Zeigler, J. P. Biersack and U. Littmark, // "The Stopping and Range of Ions in Matter" volume 1, Pergamon, 1985. diff --git a/src/MANYBODY/pair_comb.h b/src/MANYBODY/pair_comb.h index 7a3d279033..f809dfaf0e 100644 --- a/src/MANYBODY/pair_comb.h +++ b/src/MANYBODY/pair_comb.h @@ -38,7 +38,7 @@ class PairComb : public Pair { virtual double yasu_char(double *, int &); double enegtot; - static const int NPARAMS_PER_LINE = 49; + static constexpr int NPARAMS_PER_LINE = 49; protected: struct Param { diff --git a/src/MANYBODY/pair_comb3.h b/src/MANYBODY/pair_comb3.h index 567859127b..b1c216fee6 100644 --- a/src/MANYBODY/pair_comb3.h +++ b/src/MANYBODY/pair_comb3.h @@ -37,7 +37,7 @@ class PairComb3 : public Pair { virtual double combqeq(double *, int &); double enegtot; - static const int NPARAMS_PER_LINE = 74; + static constexpr int NPARAMS_PER_LINE = 74; protected: // general potential parameters diff --git a/src/MANYBODY/pair_gw.h b/src/MANYBODY/pair_gw.h index 77ff5c0399..2dbf3dcf84 100644 --- a/src/MANYBODY/pair_gw.h +++ b/src/MANYBODY/pair_gw.h @@ -34,7 +34,7 @@ class PairGW : public Pair { void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: struct Param { diff --git a/src/MANYBODY/pair_gw_zbl.h b/src/MANYBODY/pair_gw_zbl.h index 95129c4eb8..de344e94f7 100644 --- a/src/MANYBODY/pair_gw_zbl.h +++ b/src/MANYBODY/pair_gw_zbl.h @@ -29,7 +29,7 @@ class PairGWZBL : public PairGW { PairGWZBL(class LAMMPS *); ~PairGWZBL() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; private: double global_a_0; // Bohr radius for Coulomb repulsion diff --git a/src/MANYBODY/pair_nb3b_harmonic.h b/src/MANYBODY/pair_nb3b_harmonic.h index 0267b6d0e2..b68974e15f 100644 --- a/src/MANYBODY/pair_nb3b_harmonic.h +++ b/src/MANYBODY/pair_nb3b_harmonic.h @@ -34,7 +34,7 @@ class PairNb3bHarmonic : public Pair { double init_one(int, int); void init_style(); - static const int NPARAMS_PER_LINE = 6; + static constexpr int NPARAMS_PER_LINE = 6; protected: struct Param { diff --git a/src/MANYBODY/pair_sw.h b/src/MANYBODY/pair_sw.h index 441b1c0303..4fdc538430 100644 --- a/src/MANYBODY/pair_sw.h +++ b/src/MANYBODY/pair_sw.h @@ -34,7 +34,7 @@ class PairSW : public Pair { virtual double init_one(int, int); virtual void init_style(); - static const int NPARAMS_PER_LINE = 14; + static constexpr int NPARAMS_PER_LINE = 14; struct Param { double epsilon,sigma; diff --git a/src/MANYBODY/pair_tersoff.h b/src/MANYBODY/pair_tersoff.h index 8ebbc83544..eb29cc227d 100644 --- a/src/MANYBODY/pair_tersoff.h +++ b/src/MANYBODY/pair_tersoff.h @@ -34,7 +34,7 @@ class PairTersoff : public Pair { virtual void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: diff --git a/src/MANYBODY/pair_tersoff_mod.h b/src/MANYBODY/pair_tersoff_mod.h index 2cdbae9518..4f1d7a12a6 100644 --- a/src/MANYBODY/pair_tersoff_mod.h +++ b/src/MANYBODY/pair_tersoff_mod.h @@ -30,7 +30,7 @@ class PairTersoffMOD : public PairTersoff { PairTersoffMOD(class LAMMPS *); ~PairTersoffMOD() {} - static const int NPARAMS_PER_LINE = 20; + static constexpr int NPARAMS_PER_LINE = 20; protected: virtual void read_file(char *); diff --git a/src/MANYBODY/pair_tersoff_mod_c.h b/src/MANYBODY/pair_tersoff_mod_c.h index bad5f988be..4949cb76db 100644 --- a/src/MANYBODY/pair_tersoff_mod_c.h +++ b/src/MANYBODY/pair_tersoff_mod_c.h @@ -29,7 +29,7 @@ class PairTersoffMODC : public PairTersoffMOD { PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {}; ~PairTersoffMODC() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; protected: void read_file(char *); diff --git a/src/MANYBODY/pair_tersoff_zbl.h b/src/MANYBODY/pair_tersoff_zbl.h index be3f496b62..42483b7d89 100644 --- a/src/MANYBODY/pair_tersoff_zbl.h +++ b/src/MANYBODY/pair_tersoff_zbl.h @@ -29,7 +29,7 @@ class PairTersoffZBL : public PairTersoff { PairTersoffZBL(class LAMMPS *); ~PairTersoffZBL() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; protected: double global_a_0; // Bohr radius for Coulomb repulsion diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index bd250d328a..19a8e01cd9 100644 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -34,7 +34,7 @@ class PairVashishta : public Pair { double init_one(int, int); void init_style(); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; struct Param { double bigb,gamma,r0,bigc,costheta; diff --git a/src/RIGID/rigid_const.h b/src/RIGID/rigid_const.h index 3aae988197..b345a6b9dc 100644 --- a/src/RIGID/rigid_const.h +++ b/src/RIGID/rigid_const.h @@ -33,21 +33,21 @@ namespace LAMMPS_NS { TORQUE = 1<<8 }; - static const double TOLERANCE = 1.0e-6; - static const double EPSILON = 1.0e-7; - static const double BIG = 1.0e20; + static constexpr double TOLERANCE = 1.0e-6; + static constexpr double EPSILON = 1.0e-7; + static constexpr double BIG = 1.0e20; // moment of inertia prefactor for sphere - static const double SINERTIA = 0.4; + static constexpr double SINERTIA = 0.4; // moment of inertia prefactor for ellipsoid - static const double EINERTIA = 0.2; + static constexpr double EINERTIA = 0.2; // moment of inertia prefactor for line segment - static const double LINERTIA = 1.0/12.0; + static constexpr double LINERTIA = 1.0/12.0; - static const int MAXLINE = 1024; - static const int CHUNK = 1024; - static const int DELTA_BODY = 10000; - static const int ATTRIBUTE_PERBODY = 20; + static constexpr int MAXLINE = 1024; + static constexpr int CHUNK = 1024; + static constexpr int DELTA_BODY = 10000; + static constexpr int ATTRIBUTE_PERBODY = 20; } } diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp index 6d4197b59b..3c7a40c0dc 100644 --- a/src/SNAP/sna.cpp +++ b/src/SNAP/sna.cpp @@ -18,6 +18,7 @@ #include "sna.h" #include #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" #include "comm.h" @@ -25,6 +26,7 @@ using namespace std; using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; /* ---------------------------------------------------------------------- @@ -1363,196 +1365,6 @@ void SNA::destroy_twojmax_arrays() } -/* ---------------------------------------------------------------------- - factorial n, wrapper for precomputed table -------------------------------------------------------------------------- */ - -double SNA::factorial(int n) -{ - if (n < 0 || n > nmaxfactorial) { - char str[128]; - sprintf(str, "Invalid argument to factorial %d", n); - error->all(FLERR, str); - } - - return nfac_table[n]; -} - -/* ---------------------------------------------------------------------- - factorial n table, size SNA::nmaxfactorial+1 -------------------------------------------------------------------------- */ - -const double SNA::nfac_table[] = { - 1, - 1, - 2, - 6, - 24, - 120, - 720, - 5040, - 40320, - 362880, - 3628800, - 39916800, - 479001600, - 6227020800, - 87178291200, - 1307674368000, - 20922789888000, - 355687428096000, - 6.402373705728e+15, - 1.21645100408832e+17, - 2.43290200817664e+18, - 5.10909421717094e+19, - 1.12400072777761e+21, - 2.5852016738885e+22, - 6.20448401733239e+23, - 1.5511210043331e+25, - 4.03291461126606e+26, - 1.08888694504184e+28, - 3.04888344611714e+29, - 8.8417619937397e+30, - 2.65252859812191e+32, - 8.22283865417792e+33, - 2.63130836933694e+35, - 8.68331761881189e+36, - 2.95232799039604e+38, - 1.03331479663861e+40, - 3.71993326789901e+41, - 1.37637530912263e+43, - 5.23022617466601e+44, - 2.03978820811974e+46, - 8.15915283247898e+47, - 3.34525266131638e+49, - 1.40500611775288e+51, - 6.04152630633738e+52, - 2.65827157478845e+54, - 1.1962222086548e+56, - 5.50262215981209e+57, - 2.58623241511168e+59, - 1.24139155925361e+61, - 6.08281864034268e+62, - 3.04140932017134e+64, - 1.55111875328738e+66, - 8.06581751709439e+67, - 4.27488328406003e+69, - 2.30843697339241e+71, - 1.26964033536583e+73, - 7.10998587804863e+74, - 4.05269195048772e+76, - 2.35056133128288e+78, - 1.3868311854569e+80, - 8.32098711274139e+81, - 5.07580213877225e+83, - 3.14699732603879e+85, - 1.98260831540444e+87, - 1.26886932185884e+89, - 8.24765059208247e+90, - 5.44344939077443e+92, - 3.64711109181887e+94, - 2.48003554243683e+96, - 1.71122452428141e+98, - 1.19785716699699e+100, - 8.50478588567862e+101, - 6.12344583768861e+103, - 4.47011546151268e+105, - 3.30788544151939e+107, - 2.48091408113954e+109, - 1.88549470166605e+111, - 1.45183092028286e+113, - 1.13242811782063e+115, - 8.94618213078297e+116, - 7.15694570462638e+118, - 5.79712602074737e+120, - 4.75364333701284e+122, - 3.94552396972066e+124, - 3.31424013456535e+126, - 2.81710411438055e+128, - 2.42270953836727e+130, - 2.10775729837953e+132, - 1.85482642257398e+134, - 1.65079551609085e+136, - 1.48571596448176e+138, - 1.3520015276784e+140, - 1.24384140546413e+142, - 1.15677250708164e+144, - 1.08736615665674e+146, - 1.03299784882391e+148, - 9.91677934870949e+149, - 9.61927596824821e+151, - 9.42689044888324e+153, - 9.33262154439441e+155, - 9.33262154439441e+157, - 9.42594775983835e+159, - 9.61446671503512e+161, - 9.90290071648618e+163, - 1.02990167451456e+166, - 1.08139675824029e+168, - 1.14628056373471e+170, - 1.22652020319614e+172, - 1.32464181945183e+174, - 1.44385958320249e+176, - 1.58824554152274e+178, - 1.76295255109024e+180, - 1.97450685722107e+182, - 2.23119274865981e+184, - 2.54355973347219e+186, - 2.92509369349301e+188, - 3.3931086844519e+190, - 3.96993716080872e+192, - 4.68452584975429e+194, - 5.5745857612076e+196, - 6.68950291344912e+198, - 8.09429852527344e+200, - 9.8750442008336e+202, - 1.21463043670253e+205, - 1.50614174151114e+207, - 1.88267717688893e+209, - 2.37217324288005e+211, - 3.01266001845766e+213, - 3.8562048236258e+215, - 4.97450422247729e+217, - 6.46685548922047e+219, - 8.47158069087882e+221, - 1.118248651196e+224, - 1.48727070609069e+226, - 1.99294274616152e+228, - 2.69047270731805e+230, - 3.65904288195255e+232, - 5.01288874827499e+234, - 6.91778647261949e+236, - 9.61572319694109e+238, - 1.34620124757175e+241, - 1.89814375907617e+243, - 2.69536413788816e+245, - 3.85437071718007e+247, - 5.5502938327393e+249, - 8.04792605747199e+251, - 1.17499720439091e+254, - 1.72724589045464e+256, - 2.55632391787286e+258, - 3.80892263763057e+260, - 5.71338395644585e+262, - 8.62720977423323e+264, - 1.31133588568345e+267, - 2.00634390509568e+269, - 3.08976961384735e+271, - 4.78914290146339e+273, - 7.47106292628289e+275, - 1.17295687942641e+278, - 1.85327186949373e+280, - 2.94670227249504e+282, - 4.71472363599206e+284, - 7.59070505394721e+286, - 1.22969421873945e+289, - 2.0044015765453e+291, - 3.28721858553429e+293, - 5.42391066613159e+295, - 9.00369170577843e+297, - 1.503616514865e+300, // nmaxfactorial = 167 -}; - /* ---------------------------------------------------------------------- the function delta given by VMK Eq. 8.2(1) ------------------------------------------------------------------------- */ diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h index 30c4aa6bbc..746b55fb70 100644 --- a/src/SNAP/sna.h +++ b/src/SNAP/sna.h @@ -100,10 +100,6 @@ private: double** dulist_r, ** dulist_i; int elem_duarray; // element of j in derivative - static const int nmaxfactorial = 167; - static const double nfac_table[]; - double factorial(int); - void create_twojmax_arrays(); void destroy_twojmax_arrays(); void init_clebsch_gordan(); diff --git a/src/USER-MISC/dihedral_quadratic.cpp b/src/USER-MISC/dihedral_quadratic.cpp index bfe13d13d1..ef5711215c 100644 --- a/src/USER-MISC/dihedral_quadratic.cpp +++ b/src/USER-MISC/dihedral_quadratic.cpp @@ -195,6 +195,8 @@ void DihedralQuadratic::compute(int eflag, int vflag) siinv = 1.0/si; double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; p = k[type]*dphi; pd = - 2.0 * p * siinv; p = p * dphi; diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp index b29f499de8..12390ba69d 100644 --- a/src/USER-MISC/pair_edip.cpp +++ b/src/USER-MISC/pair_edip.cpp @@ -46,7 +46,7 @@ using namespace LAMMPS_NS; // max number of interaction per atom for f(Z) environment potential -#define leadDimInteractionList 64 +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ diff --git a/src/USER-MISC/pair_edip_multi.cpp b/src/USER-MISC/pair_edip_multi.cpp index 811f80abfc..9d36db7041 100644 --- a/src/USER-MISC/pair_edip_multi.cpp +++ b/src/USER-MISC/pair_edip_multi.cpp @@ -32,13 +32,11 @@ #include "error.h" #include "citeme.h" - using namespace LAMMPS_NS; #define MAXLINE 1024 #define DELTA 4 - static const char cite_pair_edip[] = "@article{cjiang2012\n" " author = {Jian, Chao and Morgan, Dane, and Szlufarska, Izabella},\n" @@ -56,7 +54,9 @@ static const char cite_pair_edip[] = " year = {2010},\n" "}\n\n"; +// max number of interaction per atom for f(Z) environment potential +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ @@ -94,7 +94,6 @@ PairEDIPMulti::~PairEDIPMulti() memory->destroy(cutsq); delete [] map; -//XXX deallocateGrids(); deallocatePreLoops(); } } diff --git a/src/USER-MISC/pair_edip_multi.h b/src/USER-MISC/pair_edip_multi.h index 26433a66fa..5120ece5d2 100644 --- a/src/USER-MISC/pair_edip_multi.h +++ b/src/USER-MISC/pair_edip_multi.h @@ -36,17 +36,17 @@ class PairEDIPMulti : public Pair { protected: struct Param { - double A, B;//coefficients for pair interaction I-J - double cutoffA;//cut-off distance for pair interaction I-J - double cutoffC;//lower cut-off distance for calculating Z_I - double alpha;//coefficient for calculating Z_I - double beta;//attractive term for pair I-J - double sigma;//cut-off coefficient for pair I-J - double rho;//pair I-J - double gamma;//coefficient for three-body interaction I-J-K - double eta, lambda;//coefficients for function h(l,Z) - double mu, Q0;//coefficients for function Q(Z) - double u1, u2, u3, u4;//coefficients for function tau(Z) + double A, B; // coefficients for pair interaction I-J + double cutoffA; // cut-off distance for pair interaction I-J + double cutoffC; // lower cut-off distance for calculating Z_I + double alpha; // coefficient for calculating Z_I + double beta; // attractive term for pair I-J + double sigma; // cut-off coefficient for pair I-J + double rho; // pair I-J + double gamma; // coefficient for three-body interaction I-J-K + double eta, lambda; // coefficients for function h(l,Z) + double mu, Q0; // coefficients for function Q(Z) + double u1, u2, u3, u4; // coefficients for function tau(Z) double cutsq; int ielement,jelement,kelement; }; @@ -62,10 +62,6 @@ class PairEDIPMulti : public Pair { int maxparam; // max # of parameter sets Param *params; // parameter set for an I-J-K interaction - // max number of interaction per atom for f(Z) environment potential - - static const int leadDimInteractionList = 64; - void allocate(); void allocatePreLoops(void); void deallocatePreLoops(void); diff --git a/src/USER-MISC/pair_tersoff_table.h b/src/USER-MISC/pair_tersoff_table.h index faa704191c..d63a5be5a6 100644 --- a/src/USER-MISC/pair_tersoff_table.h +++ b/src/USER-MISC/pair_tersoff_table.h @@ -43,7 +43,7 @@ class PairTersoffTable : public Pair { void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: struct Param { diff --git a/src/USER-OMP/dihedral_quadratic_omp.cpp b/src/USER-OMP/dihedral_quadratic_omp.cpp index 2e25258b13..d183500182 100644 --- a/src/USER-OMP/dihedral_quadratic_omp.cpp +++ b/src/USER-OMP/dihedral_quadratic_omp.cpp @@ -23,11 +23,13 @@ #include "neighbor.h" #include "force.h" #include "update.h" +#include "math_const.h" #include "error.h" #include "suffix.h" using namespace LAMMPS_NS; +using namespace MathConst; #define TOLERANCE 0.05 #define SMALL 0.001 @@ -218,6 +220,8 @@ void DihedralQuadraticOMP::eval(int nfrom, int nto, ThrData * const thr) siinv = 1.0/si; double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; p = k[type]*dphi; pd = - 2.0 * p * siinv; p = p * dphi; diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp new file mode 100644 index 0000000000..33722eb9d3 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp @@ -0,0 +1,126 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newtoff_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.h b/src/USER-OMP/npair_half_size_multi_newtoff_omp.h new file mode 100644 index 0000000000..dd883d7f3c --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.h @@ -0,0 +1,44 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newtoff/omp, + NPairHalfSizeMultiNewtoffOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_OMP | + NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtoffOmp : public NPair { + public: + NPairHalfSizeMultiNewtoffOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtoffOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp new file mode 100644 index 0000000000..0be87457e3 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp @@ -0,0 +1,151 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newton_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.h b/src/USER-OMP/npair_half_size_multi_newton_omp.h new file mode 100644 index 0000000000..03d712145f --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton/omp, + NPairHalfSizeMultiNewtonOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_OMP | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonOmp : public NPair { + public: + NPairHalfSizeMultiNewtonOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtonOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp new file mode 100644 index 0000000000..ef26a5f2bd --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp @@ -0,0 +1,134 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newton_tri_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h new file mode 100644 index 0000000000..6e936c8da4 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton/tri/omp, + NPairHalfSizeMultiNewtonTriOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI | NP_OMP) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonTriOmp : public NPair { + public: + NPairHalfSizeMultiNewtonTriOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtonTriOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/USER-OMP/pair_edip_omp.cpp b/src/USER-OMP/pair_edip_omp.cpp index 2eb8384113..106d038743 100644 --- a/src/USER-OMP/pair_edip_omp.cpp +++ b/src/USER-OMP/pair_edip_omp.cpp @@ -28,7 +28,7 @@ using namespace LAMMPS_NS; // max number of interaction per atom for f(Z) environment potential -#define leadDimInteractionList 64 +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj_cubic_omp.cpp b/src/USER-OMP/pair_lj_cubic_omp.cpp index 872574a3b2..e96b3d0a79 100644 --- a/src/USER-OMP/pair_lj_cubic_omp.cpp +++ b/src/USER-OMP/pair_lj_cubic_omp.cpp @@ -23,6 +23,8 @@ #include #include "omp_compat.h" +#include "pair_lj_cubic_const.h" + using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/USER-OMP/pair_zbl_omp.cpp b/src/USER-OMP/pair_zbl_omp.cpp index 993e75f499..6794789f25 100644 --- a/src/USER-OMP/pair_zbl_omp.cpp +++ b/src/USER-OMP/pair_zbl_omp.cpp @@ -24,7 +24,6 @@ #include "omp_compat.h" using namespace LAMMPS_NS; -using namespace PairZBLConstants; /* ---------------------------------------------------------------------- */ diff --git a/src/USER-SMTBQ/pair_smtbq.cpp b/src/USER-SMTBQ/pair_smtbq.cpp index 94647cb552..aedf90da07 100644 --- a/src/USER-SMTBQ/pair_smtbq.cpp +++ b/src/USER-SMTBQ/pair_smtbq.cpp @@ -72,17 +72,6 @@ using namespace MathSpecial; #define PGDELTA 1 #define MAXNEIGH 24 -/* ------------------------------------------------------------------------------------ - - Calculates the factorial of an integer n via recursion. - - ------------------------------------------------------------------------------------ */ -static double factorial(int n) -{ - if (n <= 1) return 1.0; - else return static_cast(n)*factorial(n-1); -} - /* ---------------------------------------------------------------------- */ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp) diff --git a/src/balance.cpp b/src/balance.cpp index c47412cc9e..6294e023b3 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -924,66 +924,66 @@ int Balance::shift() i = 0; while (i < np) { if (split[i+1] - split[i] < close) { - j = i+1; + j = i+1; - // I,J = set of consecutive splits that are collectively too close - // if can expand set and not become too close to splits I-1 or J+1, do it - // else add split I-1 or J+1 to set and try again - // delta = size of expanded split set that will satisy criterion + // I,J = set of consecutive splits that are collectively too close + // if can expand set and not become too close to splits I-1 or J+1, do it + // else add split I-1 or J+1 to set and try again + // delta = size of expanded split set that will satisy criterion - while (1) { - delta = (j-i) * close; - midpt = 0.5 * (split[i]+split[j]); - start = midpt - 0.5*delta; - stop = midpt + 0.5*delta; + while (1) { + delta = (j-i) * close; + midpt = 0.5 * (split[i]+split[j]); + start = midpt - 0.5*delta; + stop = midpt + 0.5*delta; - if (i > 0) lbound = split[i-1] + close; - else lbound = 0.0; - if (j < np) ubound = split[j+1] - close; - else ubound = 1.0; + if (i > 0) lbound = split[i-1] + close; + else lbound = 0.0; + if (j < np) ubound = split[j+1] - close; + else ubound = 1.0; - // start/stop are within bounds, reset the splits + // start/stop are within bounds, reset the splits - if (start >= lbound && stop <= ubound) break; + if (start >= lbound && stop <= ubound) break; - // try a shift to either bound, reset the splits if delta fits - // these tests change start/stop + // try a shift to either bound, reset the splits if delta fits + // these tests change start/stop - if (start < lbound) { - start = lbound; - stop = start + delta; - if (stop <= ubound) break; - } else if (stop > ubound) { - stop = ubound; - start = stop - delta; - if (start >= lbound) break; - } + if (start < lbound) { + start = lbound; + stop = start + delta; + if (stop <= ubound) break; + } else if (stop > ubound) { + stop = ubound; + start = stop - delta; + if (start >= lbound) break; + } - // delta does not fit between lbound and ubound - // exit if can't expand set, else expand set - // if can expand in either direction, - // pick new split closest to current midpt of set + // delta does not fit between lbound and ubound + // exit if can't expand set, else expand set + // if can expand in either direction, + // pick new split closest to current midpt of set - if (i == 0 && j == np) { - start = 0.0; stop = 1.0; - break; - } - if (i == 0) j++; - else if (j == np) i--; - else if (midpt-lbound < ubound-midpt) i--; - else j++; - } + if (i == 0 && j == np) { + start = 0.0; stop = 1.0; + break; + } + if (i == 0) j++; + else if (j == np) i--; + else if (midpt-lbound < ubound-midpt) i--; + else j++; + } - // reset all splits between I,J inclusive to be equi-spaced + // reset all splits between I,J inclusive to be equi-spaced - spacing = (stop-start) / (j-i); - for (m = i; m <= j; m++) - split[m] = start + (m-i)*spacing; + spacing = (stop-start) / (j-i); + for (m = i; m <= j; m++) + split[m] = start + (m-i)*spacing; if (j == np) split[np] = 1.0; - // continue testing beyond the J split + // continue testing beyond the J split - i = j+1; + i = j+1; } else i++; } diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index db8a2f404b..b334654e82 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -23,6 +23,7 @@ #include "error.h" #include "force.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" @@ -36,6 +37,8 @@ using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; + #ifdef DBL_EPSILON #define MY_EPSILON (10.0*DBL_EPSILON) @@ -708,191 +711,3 @@ void ComputeOrientOrderAtom::init_clebsch_gordan() } } } - -/* ---------------------------------------------------------------------- - factorial n, wrapper for precomputed table -------------------------------------------------------------------------- */ - -double ComputeOrientOrderAtom::factorial(int n) -{ - if (n < 0 || n > nmaxfactorial) - error->all(FLERR,fmt::format("Invalid argument to factorial {}", n)); - - return nfac_table[n]; -} - -/* ---------------------------------------------------------------------- - factorial n table, size SNA::nmaxfactorial+1 -------------------------------------------------------------------------- */ - -const double ComputeOrientOrderAtom::nfac_table[] = { - 1, - 1, - 2, - 6, - 24, - 120, - 720, - 5040, - 40320, - 362880, - 3628800, - 39916800, - 479001600, - 6227020800, - 87178291200, - 1307674368000, - 20922789888000, - 355687428096000, - 6.402373705728e+15, - 1.21645100408832e+17, - 2.43290200817664e+18, - 5.10909421717094e+19, - 1.12400072777761e+21, - 2.5852016738885e+22, - 6.20448401733239e+23, - 1.5511210043331e+25, - 4.03291461126606e+26, - 1.08888694504184e+28, - 3.04888344611714e+29, - 8.8417619937397e+30, - 2.65252859812191e+32, - 8.22283865417792e+33, - 2.63130836933694e+35, - 8.68331761881189e+36, - 2.95232799039604e+38, - 1.03331479663861e+40, - 3.71993326789901e+41, - 1.37637530912263e+43, - 5.23022617466601e+44, - 2.03978820811974e+46, - 8.15915283247898e+47, - 3.34525266131638e+49, - 1.40500611775288e+51, - 6.04152630633738e+52, - 2.65827157478845e+54, - 1.1962222086548e+56, - 5.50262215981209e+57, - 2.58623241511168e+59, - 1.24139155925361e+61, - 6.08281864034268e+62, - 3.04140932017134e+64, - 1.55111875328738e+66, - 8.06581751709439e+67, - 4.27488328406003e+69, - 2.30843697339241e+71, - 1.26964033536583e+73, - 7.10998587804863e+74, - 4.05269195048772e+76, - 2.35056133128288e+78, - 1.3868311854569e+80, - 8.32098711274139e+81, - 5.07580213877225e+83, - 3.14699732603879e+85, - 1.98260831540444e+87, - 1.26886932185884e+89, - 8.24765059208247e+90, - 5.44344939077443e+92, - 3.64711109181887e+94, - 2.48003554243683e+96, - 1.71122452428141e+98, - 1.19785716699699e+100, - 8.50478588567862e+101, - 6.12344583768861e+103, - 4.47011546151268e+105, - 3.30788544151939e+107, - 2.48091408113954e+109, - 1.88549470166605e+111, - 1.45183092028286e+113, - 1.13242811782063e+115, - 8.94618213078297e+116, - 7.15694570462638e+118, - 5.79712602074737e+120, - 4.75364333701284e+122, - 3.94552396972066e+124, - 3.31424013456535e+126, - 2.81710411438055e+128, - 2.42270953836727e+130, - 2.10775729837953e+132, - 1.85482642257398e+134, - 1.65079551609085e+136, - 1.48571596448176e+138, - 1.3520015276784e+140, - 1.24384140546413e+142, - 1.15677250708164e+144, - 1.08736615665674e+146, - 1.03299784882391e+148, - 9.91677934870949e+149, - 9.61927596824821e+151, - 9.42689044888324e+153, - 9.33262154439441e+155, - 9.33262154439441e+157, - 9.42594775983835e+159, - 9.61446671503512e+161, - 9.90290071648618e+163, - 1.02990167451456e+166, - 1.08139675824029e+168, - 1.14628056373471e+170, - 1.22652020319614e+172, - 1.32464181945183e+174, - 1.44385958320249e+176, - 1.58824554152274e+178, - 1.76295255109024e+180, - 1.97450685722107e+182, - 2.23119274865981e+184, - 2.54355973347219e+186, - 2.92509369349301e+188, - 3.3931086844519e+190, - 3.96993716080872e+192, - 4.68452584975429e+194, - 5.5745857612076e+196, - 6.68950291344912e+198, - 8.09429852527344e+200, - 9.8750442008336e+202, - 1.21463043670253e+205, - 1.50614174151114e+207, - 1.88267717688893e+209, - 2.37217324288005e+211, - 3.01266001845766e+213, - 3.8562048236258e+215, - 4.97450422247729e+217, - 6.46685548922047e+219, - 8.47158069087882e+221, - 1.118248651196e+224, - 1.48727070609069e+226, - 1.99294274616152e+228, - 2.69047270731805e+230, - 3.65904288195255e+232, - 5.01288874827499e+234, - 6.91778647261949e+236, - 9.61572319694109e+238, - 1.34620124757175e+241, - 1.89814375907617e+243, - 2.69536413788816e+245, - 3.85437071718007e+247, - 5.5502938327393e+249, - 8.04792605747199e+251, - 1.17499720439091e+254, - 1.72724589045464e+256, - 2.55632391787286e+258, - 3.80892263763057e+260, - 5.71338395644585e+262, - 8.62720977423323e+264, - 1.31133588568345e+267, - 2.00634390509568e+269, - 3.08976961384735e+271, - 4.78914290146339e+273, - 7.47106292628289e+275, - 1.17295687942641e+278, - 1.85327186949373e+280, - 2.94670227249504e+282, - 4.71472363599206e+284, - 7.59070505394721e+286, - 1.22969421873945e+289, - 2.0044015765453e+291, - 3.28721858553429e+293, - 5.42391066613159e+295, - 9.00369170577843e+297, - 1.503616514865e+300, // nmaxfactorial = 167 -}; - diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h index 4e5084bfcd..4fc122c5c8 100644 --- a/src/compute_orientorder_atom.h +++ b/src/compute_orientorder_atom.h @@ -56,9 +56,6 @@ class ComputeOrientOrderAtom : public Compute { double polar_prefactor(int, int, double); double associated_legendre(int, int, double); - static const int nmaxfactorial = 167; - static const double nfac_table[]; - double factorial(int); virtual void init_clebsch_gordan(); double *cglist; // Clebsch-Gordan coeffs int idxcg_max; diff --git a/src/fix.h b/src/fix.h index 75c644e87f..ec951e3bb2 100644 --- a/src/fix.h +++ b/src/fix.h @@ -251,31 +251,32 @@ class Fix : protected Pointers { }; namespace FixConst { - static const int INITIAL_INTEGRATE = 1<<0; - static const int POST_INTEGRATE = 1<<1; - static const int PRE_EXCHANGE = 1<<2; - static const int PRE_NEIGHBOR = 1<<3; - static const int POST_NEIGHBOR = 1<<4; - static const int PRE_FORCE = 1<<5; - static const int PRE_REVERSE = 1<<6; - static const int POST_FORCE = 1<<7; - static const int FINAL_INTEGRATE = 1<<8; - static const int END_OF_STEP = 1<<9; - static const int POST_RUN = 1<<10; - static const int THERMO_ENERGY = 1<<11; - static const int INITIAL_INTEGRATE_RESPA = 1<<12; - static const int POST_INTEGRATE_RESPA = 1<<13; - static const int PRE_FORCE_RESPA = 1<<14; - static const int POST_FORCE_RESPA = 1<<15; - static const int FINAL_INTEGRATE_RESPA = 1<<16; - static const int MIN_PRE_EXCHANGE = 1<<17; - static const int MIN_PRE_NEIGHBOR = 1<<18; - static const int MIN_POST_NEIGHBOR = 1<<19; - static const int MIN_PRE_FORCE = 1<<20; - static const int MIN_PRE_REVERSE = 1<<21; - static const int MIN_POST_FORCE = 1<<22; - static const int MIN_ENERGY = 1<<23; - static const int FIX_CONST_LAST = 1<<24; + enum { + INITIAL_INTEGRATE = 1<<0, + POST_INTEGRATE = 1<<1, + PRE_EXCHANGE = 1<<2, + PRE_NEIGHBOR = 1<<3, + POST_NEIGHBOR = 1<<4, + PRE_FORCE = 1<<5, + PRE_REVERSE = 1<<6, + POST_FORCE = 1<<7, + FINAL_INTEGRATE = 1<<8, + END_OF_STEP = 1<<9, + POST_RUN = 1<<10, + THERMO_ENERGY = 1<<11, + INITIAL_INTEGRATE_RESPA = 1<<12, + POST_INTEGRATE_RESPA = 1<<13, + PRE_FORCE_RESPA = 1<<14, + POST_FORCE_RESPA = 1<<15, + FINAL_INTEGRATE_RESPA = 1<<16, + MIN_PRE_EXCHANGE = 1<<17, + MIN_PRE_NEIGHBOR = 1<<18, + MIN_POST_NEIGHBOR = 1<<19, + MIN_PRE_FORCE = 1<<20, + MIN_PRE_REVERSE = 1<<21, + MIN_POST_FORCE = 1<<22, + MIN_ENERGY = 1<<23 + }; } } diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp index f01de7c2fa..3b522c185f 100644 --- a/src/fix_temp_csld.cpp +++ b/src/fix_temp_csld.cpp @@ -305,16 +305,17 @@ double FixTempCSLD::compute_scalar() void FixTempCSLD::write_restart(FILE *fp) { - int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy + const int PRNGSIZE = 98+2+3; + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; list[0] = energy; list[1] = comm->nprocs; } - double state[103]; + double state[PRNGSIZE]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world); + MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp index 0130631172..674c436e9d 100644 --- a/src/fix_temp_csvr.cpp +++ b/src/fix_temp_csvr.cpp @@ -338,16 +338,17 @@ double FixTempCSVR::compute_scalar() void FixTempCSVR::write_restart(FILE *fp) { - int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy + const int PRNGSIZE = 98+2+3; + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; list[0] = energy; list[1] = comm->nprocs; } - double state[103]; + double state[PRNGSIZE]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world); + MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); diff --git a/src/math_const.h b/src/math_const.h index fefb442f68..53a485bf99 100644 --- a/src/math_const.h +++ b/src/math_const.h @@ -17,18 +17,18 @@ namespace LAMMPS_NS { namespace MathConst { - static const double THIRD = 1.0/3.0; - static const double TWOTHIRDS = 2.0/3.0; - static const double MY_PI = 3.14159265358979323846; // pi - static const double MY_2PI = 6.28318530717958647692; // 2pi - static const double MY_3PI = 9.42477796076937971538; // 3pi - static const double MY_4PI = 12.56637061435917295384; // 4pi - static const double MY_PI2 = 1.57079632679489661923; // pi/2 - static const double MY_PI4 = 0.78539816339744830962; // pi/4 - static const double MY_PIS = 1.77245385090551602729; // sqrt(pi) - static const double MY_ISPI4 = 1.12837916709551257389; // 1/sqrt(pi/4) - static const double MY_SQRT2 = 1.41421356237309504880; // sqrt(2) - static const double MY_CBRT2 = 1.25992104989487316476; // 2*(1/3) + static constexpr double THIRD = 1.0/3.0; + static constexpr double TWOTHIRDS = 2.0/3.0; + static constexpr double MY_PI = 3.14159265358979323846; // pi + static constexpr double MY_2PI = 6.28318530717958647692; // 2pi + static constexpr double MY_3PI = 9.42477796076937971538; // 3pi + static constexpr double MY_4PI = 12.56637061435917295384; // 4pi + static constexpr double MY_PI2 = 1.57079632679489661923; // pi/2 + static constexpr double MY_PI4 = 0.78539816339744830962; // pi/4 + static constexpr double MY_PIS = 1.77245385090551602729; // sqrt(pi) + static constexpr double MY_ISPI4 = 1.12837916709551257389; // 1/sqrt(pi/4) + static constexpr double MY_SQRT2 = 1.41421356237309504880; // sqrt(2) + static constexpr double MY_CBRT2 = 1.25992104989487316476; // 2*(1/3) } } diff --git a/src/math_special.cpp b/src/math_special.cpp index 5ec0fc3f61..243d5a05f3 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -1,9 +1,203 @@ #include "math_special.h" + #include #include // IWYU pragma: keep +#include + +#include "error.h" using namespace LAMMPS_NS; +static constexpr int nmaxfactorial = 167; + +/* ---------------------------------------------------------------------- + factorial n table, size nmaxfactorial+1 +------------------------------------------------------------------------- */ + +static const double nfac_table[] = { + 1, + 1, + 2, + 6, + 24, + 120, + 720, + 5040, + 40320, + 362880, + 3628800, + 39916800, + 479001600, + 6227020800, + 87178291200, + 1307674368000, + 20922789888000, + 355687428096000, + 6.402373705728e+15, + 1.21645100408832e+17, + 2.43290200817664e+18, + 5.10909421717094e+19, + 1.12400072777761e+21, + 2.5852016738885e+22, + 6.20448401733239e+23, + 1.5511210043331e+25, + 4.03291461126606e+26, + 1.08888694504184e+28, + 3.04888344611714e+29, + 8.8417619937397e+30, + 2.65252859812191e+32, + 8.22283865417792e+33, + 2.63130836933694e+35, + 8.68331761881189e+36, + 2.95232799039604e+38, + 1.03331479663861e+40, + 3.71993326789901e+41, + 1.37637530912263e+43, + 5.23022617466601e+44, + 2.03978820811974e+46, + 8.15915283247898e+47, + 3.34525266131638e+49, + 1.40500611775288e+51, + 6.04152630633738e+52, + 2.65827157478845e+54, + 1.1962222086548e+56, + 5.50262215981209e+57, + 2.58623241511168e+59, + 1.24139155925361e+61, + 6.08281864034268e+62, + 3.04140932017134e+64, + 1.55111875328738e+66, + 8.06581751709439e+67, + 4.27488328406003e+69, + 2.30843697339241e+71, + 1.26964033536583e+73, + 7.10998587804863e+74, + 4.05269195048772e+76, + 2.35056133128288e+78, + 1.3868311854569e+80, + 8.32098711274139e+81, + 5.07580213877225e+83, + 3.14699732603879e+85, + 1.98260831540444e+87, + 1.26886932185884e+89, + 8.24765059208247e+90, + 5.44344939077443e+92, + 3.64711109181887e+94, + 2.48003554243683e+96, + 1.71122452428141e+98, + 1.19785716699699e+100, + 8.50478588567862e+101, + 6.12344583768861e+103, + 4.47011546151268e+105, + 3.30788544151939e+107, + 2.48091408113954e+109, + 1.88549470166605e+111, + 1.45183092028286e+113, + 1.13242811782063e+115, + 8.94618213078297e+116, + 7.15694570462638e+118, + 5.79712602074737e+120, + 4.75364333701284e+122, + 3.94552396972066e+124, + 3.31424013456535e+126, + 2.81710411438055e+128, + 2.42270953836727e+130, + 2.10775729837953e+132, + 1.85482642257398e+134, + 1.65079551609085e+136, + 1.48571596448176e+138, + 1.3520015276784e+140, + 1.24384140546413e+142, + 1.15677250708164e+144, + 1.08736615665674e+146, + 1.03299784882391e+148, + 9.91677934870949e+149, + 9.61927596824821e+151, + 9.42689044888324e+153, + 9.33262154439441e+155, + 9.33262154439441e+157, + 9.42594775983835e+159, + 9.61446671503512e+161, + 9.90290071648618e+163, + 1.02990167451456e+166, + 1.08139675824029e+168, + 1.14628056373471e+170, + 1.22652020319614e+172, + 1.32464181945183e+174, + 1.44385958320249e+176, + 1.58824554152274e+178, + 1.76295255109024e+180, + 1.97450685722107e+182, + 2.23119274865981e+184, + 2.54355973347219e+186, + 2.92509369349301e+188, + 3.3931086844519e+190, + 3.96993716080872e+192, + 4.68452584975429e+194, + 5.5745857612076e+196, + 6.68950291344912e+198, + 8.09429852527344e+200, + 9.8750442008336e+202, + 1.21463043670253e+205, + 1.50614174151114e+207, + 1.88267717688893e+209, + 2.37217324288005e+211, + 3.01266001845766e+213, + 3.8562048236258e+215, + 4.97450422247729e+217, + 6.46685548922047e+219, + 8.47158069087882e+221, + 1.118248651196e+224, + 1.48727070609069e+226, + 1.99294274616152e+228, + 2.69047270731805e+230, + 3.65904288195255e+232, + 5.01288874827499e+234, + 6.91778647261949e+236, + 9.61572319694109e+238, + 1.34620124757175e+241, + 1.89814375907617e+243, + 2.69536413788816e+245, + 3.85437071718007e+247, + 5.5502938327393e+249, + 8.04792605747199e+251, + 1.17499720439091e+254, + 1.72724589045464e+256, + 2.55632391787286e+258, + 3.80892263763057e+260, + 5.71338395644585e+262, + 8.62720977423323e+264, + 1.31133588568345e+267, + 2.00634390509568e+269, + 3.08976961384735e+271, + 4.78914290146339e+273, + 7.47106292628289e+275, + 1.17295687942641e+278, + 1.85327186949373e+280, + 2.94670227249504e+282, + 4.71472363599206e+284, + 7.59070505394721e+286, + 1.22969421873945e+289, + 2.0044015765453e+291, + 3.28721858553429e+293, + 5.42391066613159e+295, + 9.00369170577843e+297, + 1.503616514865e+300, // nmaxfactorial = 167 +}; + +/* ---------------------------------------------------------------------- + factorial n vial lookup from precomputed table +------------------------------------------------------------------------- */ + +double MathSpecial::factorial(const int n) +{ + if (n < 0 || n > nmaxfactorial) + return std::numeric_limits::quiet_NaN(); + + return nfac_table[n]; +} + + /* Library libcerf: * Compute complex error functions, based on a new implementation of * Faddeeva's w_of_z. Also provide Dawson and Voigt functions. diff --git a/src/math_special.h b/src/math_special.h index 1e7b10d9fd..59517a2f76 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -20,6 +20,10 @@ namespace LAMMPS_NS { namespace MathSpecial { + // tabulated factorial function + + extern double factorial(const int); + // support function for scaled error function complement extern double erfcx_y100(const double y100); diff --git a/src/neighbor.h b/src/neighbor.h index 9ee2af9c75..b9b40bcf1a 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -235,49 +235,55 @@ class Neighbor : protected Pointers { }; namespace NeighConst { - static const int NB_INTEL = 1<<0; - static const int NB_KOKKOS_DEVICE = 1<<1; - static const int NB_KOKKOS_HOST = 1<<2; - static const int NB_SSA = 1<<3; + enum { + NB_INTEL = 1<<0, + NB_KOKKOS_DEVICE = 1<<1, + NB_KOKKOS_HOST = 1<<2, + NB_SSA = 1<<3 + }; - static const int NS_BIN = 1<<0; - static const int NS_MULTI = 1<<1; - static const int NS_HALF = 1<<2; - static const int NS_FULL = 1<<3; - static const int NS_2D = 1<<4; - static const int NS_3D = 1<<5; - static const int NS_NEWTON = 1<<6; - static const int NS_NEWTOFF = 1<<7; - static const int NS_ORTHO = 1<<8; - static const int NS_TRI = 1<<9; - static const int NS_GHOST = 1<<10; - static const int NS_SSA = 1<<11; + enum { + NS_BIN = 1<<0, + NS_MULTI = 1<<1, + NS_HALF = 1<<2, + NS_FULL = 1<<3, + NS_2D = 1<<4, + NS_3D = 1<<5, + NS_NEWTON = 1<<6, + NS_NEWTOFF = 1<<7, + NS_ORTHO = 1<<8, + NS_TRI = 1<<9, + NS_GHOST = 1<<10, + NS_SSA = 1<<11 + }; - static const int NP_NSQ = 1<<0; - static const int NP_BIN = 1<<1; - static const int NP_MULTI = 1<<2; - static const int NP_HALF = 1<<3; - static const int NP_FULL = 1<<4; - static const int NP_ORTHO = 1<<5; - static const int NP_TRI = 1<<6; - static const int NP_ATOMONLY = 1<<7; - static const int NP_MOLONLY = 1<<8; - static const int NP_NEWTON = 1<<9; - static const int NP_NEWTOFF = 1<<10; - static const int NP_GHOST = 1<<11; - static const int NP_SIZE = 1<<12; - static const int NP_ONESIDE = 1<<13; - static const int NP_RESPA = 1<<14; - static const int NP_BOND = 1<<15; - static const int NP_OMP = 1<<16; - static const int NP_INTEL = 1<<17; - static const int NP_KOKKOS_DEVICE = 1<<18; - static const int NP_KOKKOS_HOST = 1<<19; - static const int NP_SSA = 1<<20; - static const int NP_COPY = 1<<21; - static const int NP_SKIP = 1<<22; - static const int NP_HALF_FULL = 1<<23; - static const int NP_OFF2ON = 1<<24; + enum { + NP_NSQ = 1<<0, + NP_BIN = 1<<1, + NP_MULTI = 1<<2, + NP_HALF = 1<<3, + NP_FULL = 1<<4, + NP_ORTHO = 1<<5, + NP_TRI = 1<<6, + NP_ATOMONLY = 1<<7, + NP_MOLONLY = 1<<8, + NP_NEWTON = 1<<9, + NP_NEWTOFF = 1<<10, + NP_GHOST = 1<<11, + NP_SIZE = 1<<12, + NP_ONESIDE = 1<<13, + NP_RESPA = 1<<14, + NP_BOND = 1<<15, + NP_OMP = 1<<16, + NP_INTEL = 1<<17, + NP_KOKKOS_DEVICE = 1<<18, + NP_KOKKOS_HOST = 1<<19, + NP_SSA = 1<<20, + NP_COPY = 1<<21, + NP_SKIP = 1<<22, + NP_HALF_FULL = 1<<23, + NP_OFF2ON = 1<<24 + }; } } diff --git a/src/npair_half_size_multi_newtoff.cpp b/src/npair_half_size_multi_newtoff.cpp new file mode 100644 index 0000000000..813a642b96 --- /dev/null +++ b/src/npair_half_size_multi_newtoff.cpp @@ -0,0 +1,117 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "npair_half_size_multi_newtoff.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtoff::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/npair_half_size_multi_newtoff.h b/src/npair_half_size_multi_newtoff.h new file mode 100644 index 0000000000..f255f9a17d --- /dev/null +++ b/src/npair_half_size_multi_newtoff.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newtoff, + NPairHalfSizeMultiNewtoff, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtoff : public NPair { + public: + NPairHalfSizeMultiNewtoff(class LAMMPS *); + ~NPairHalfSizeMultiNewtoff() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED +*/ diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp new file mode 100644 index 0000000000..943965ab63 --- /dev/null +++ b/src/npair_half_size_multi_newton.cpp @@ -0,0 +1,143 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "npair_half_size_multi_newton.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewton::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/npair_half_size_multi_newton.h b/src/npair_half_size_multi_newton.h new file mode 100644 index 0000000000..3e3d6f4180 --- /dev/null +++ b/src/npair_half_size_multi_newton.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton, + NPairHalfSizeMultiNewton, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewton : public NPair { + public: + NPairHalfSizeMultiNewton(class LAMMPS *); + ~NPairHalfSizeMultiNewton() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED +*/ diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp new file mode 100644 index 0000000000..e59a3f088a --- /dev/null +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -0,0 +1,125 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "npair_half_size_multi_newton_tri.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonTri::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/npair_half_size_multi_newton_tri.h b/src/npair_half_size_multi_newton_tri.h new file mode 100644 index 0000000000..6afe8201a7 --- /dev/null +++ b/src/npair_half_size_multi_newton_tri.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton/tri, + NPairHalfSizeMultiNewtonTri, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonTri : public NPair { + public: + NPairHalfSizeMultiNewtonTri(class LAMMPS *); + ~NPairHalfSizeMultiNewtonTri() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED +*/ diff --git a/src/pair_coul_streitz.h b/src/pair_coul_streitz.h index d5a4a2de5c..2f62846212 100644 --- a/src/pair_coul_streitz.h +++ b/src/pair_coul_streitz.h @@ -36,7 +36,7 @@ class PairCoulStreitz : public Pair { double memory_usage(); virtual void *extract(const char *, int &); - static const int NPARAMS_PER_LINE = 6; + static constexpr int NPARAMS_PER_LINE = 6; protected: struct Param { diff --git a/src/pair_lj_cubic.cpp b/src/pair_lj_cubic.cpp index 7fe4757358..c5f1f5b257 100644 --- a/src/pair_lj_cubic.cpp +++ b/src/pair_lj_cubic.cpp @@ -26,6 +26,7 @@ #include "memory.h" #include "error.h" +#include "pair_lj_cubic_const.h" using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/pair_lj_cubic.h b/src/pair_lj_cubic.h index d0af06c0a0..4b5edf7e97 100644 --- a/src/pair_lj_cubic.h +++ b/src/pair_lj_cubic.h @@ -46,20 +46,6 @@ class PairLJCubic : public Pair { void allocate(); }; - -namespace PairLJCubicConstants { - - // LJ quantities scaled by epsilon and rmin = sigma*2^1/6 - - static const double RT6TWO = 1.1224621; // 2^1/6 - static const double SS = 1.1086834; // inflection point (13/7)^1/6 - static const double PHIS = -0.7869823; // energy at s - static const double DPHIDS = 2.6899009; // gradient at s - static const double A3 = 27.93357; // cubic coefficient - static const double SM = 1.5475375; // cubic cutoff = s*67/48 - -} - } #endif diff --git a/src/pair_lj_cubic_const.h b/src/pair_lj_cubic_const.h new file mode 100644 index 0000000000..5fc773ac3c --- /dev/null +++ b/src/pair_lj_cubic_const.h @@ -0,0 +1,48 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_PAIR_LJ_CUBIC_CONST_H +#define LMP_PAIR_LJ_CUBIC_CONST_H + +namespace LAMMPS_NS { +namespace PairLJCubicConstants { + + // LJ quantities scaled by epsilon and rmin = sigma*2^1/6 + + static constexpr double RT6TWO = 1.1224620483093730; // 2^1/6 + static constexpr double SS = 1.1086834179687215; // inflection point (13/7)^1/6 + static constexpr double PHIS = -0.7869822485207097; // energy at s + static constexpr double DPHIDS = 2.6899008972047196; // gradient at s + static constexpr double A3 = 27.9335700460986445; // cubic coefficient + static constexpr double SM = 1.5475372709146737; // cubic cutoff = s*67/48} +} +} +#endif + +// python script to compute the constants +// +// sixth = 1.0/6.0 +// rmin = pow(2.0,sixth) +// rs = pow(26.0/7.0,sixth) +// ss = rs/rmin +// pow6 = pow(1.0/rs,6.0) +// phis = 4.0*pow6*(pow6-1.0) +// dpds = -24.0*pow6*(2.0*pow6-1.0)/rs*rmin +// a3 = 8.0*pow(dpds,3)/(9.0*phis*phis) +// sm = 67.0/48.0*ss +// print("static constexpr double RT6TWO = %19.16f; // 2^1/6" % rmin) +// print("static constexpr double SS = %19.16f; // inflection point (13/7)^1/6" % ss) +// print("static constexpr double PHIS = %19.16f; // energy at s" % phis) +// print("static constexpr double DPHIDS = %19.16f; // gradient at s" % dpds) +// print("static constexpr double A3 = %19.16f; // cubic coefficient" % a3) +// print("static constexpr double SM = %19.16f; // cubic cutoff = s*67/48" % sm) diff --git a/src/pair_zbl.cpp b/src/pair_zbl.cpp index 6dff702d2c..80b623ff35 100644 --- a/src/pair_zbl.cpp +++ b/src/pair_zbl.cpp @@ -16,15 +16,18 @@ ------------------------------------------------------------------------- */ #include "pair_zbl.h" + #include + #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" +#include "memory.h" #include "neighbor.h" #include "neigh_list.h" -#include "memory.h" -#include "error.h" +#include "pair_zbl_const.h" // From J.F. Zeigler, J. P. Biersack and U. Littmark, // "The Stopping and Range of Ions in Matter" volume 1, Pergamon, 1985. diff --git a/src/pair_zbl.h b/src/pair_zbl.h index 55a5dc5fa4..6a16bc7419 100644 --- a/src/pair_zbl.h +++ b/src/pair_zbl.h @@ -54,23 +54,6 @@ class PairZBL : public Pair { double d2zbldr2(double, int, int); void set_coeff(int, int, double, double); }; - -namespace PairZBLConstants { - - // ZBL constants - - static const double pzbl = 0.23; - static const double a0 = 0.46850; - static const double c1 = 0.02817; - static const double c2 = 0.28022; - static const double c3 = 0.50986; - static const double c4 = 0.18175; - static const double d1 = 0.20162; - static const double d2 = 0.40290; - static const double d3 = 0.94229; - static const double d4 = 3.19980; -} - } #endif diff --git a/src/pair_zbl_const.h b/src/pair_zbl_const.h new file mode 100644 index 0000000000..385657693f --- /dev/null +++ b/src/pair_zbl_const.h @@ -0,0 +1,34 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_PAIR_ZBL_CONST_H +#define LMP_PAIR_ZBL_CONST_H + +namespace LAMMPS_NS { +namespace PairZBLConstants { + + // ZBL constants + + static constexpr double pzbl = 0.23; + static constexpr double a0 = 0.46850; + static constexpr double c1 = 0.02817; + static constexpr double c2 = 0.28022; + static constexpr double c3 = 0.50986; + static constexpr double c4 = 0.18175; + static constexpr double d1 = 0.20162; + static constexpr double d2 = 0.40290; + static constexpr double d3 = 0.94229; + static constexpr double d4 = 3.19980; +} +} +#endif diff --git a/src/text_file_reader.h b/src/text_file_reader.h index 0b90304911..0da21e4581 100644 --- a/src/text_file_reader.h +++ b/src/text_file_reader.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS class TextFileReader { std::string filename; std::string filetype; - static const int MAXLINE = 1024; + static constexpr int MAXLINE = 1024; char line[MAXLINE]; FILE *fp; diff --git a/unittest/c-library/test_library_mpi.cpp b/unittest/c-library/test_library_mpi.cpp index cf33ed13b5..7502da767a 100644 --- a/unittest/c-library/test_library_mpi.cpp +++ b/unittest/c-library/test_library_mpi.cpp @@ -173,7 +173,7 @@ TEST(MPI, multi_partition) MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &me); - const char *args[] = {"LAMMPS_test", "-log", "none", "-partition", "2x2", + const char *args[] = {"LAMMPS_test", "-log", "none", "-partition", "4x1", "-echo", "screen", "-nocite", "-in", "none"}; char **argv = (char **)args; int argc = sizeof(args) / sizeof(char *); @@ -183,19 +183,15 @@ TEST(MPI, multi_partition) lammps_command(lmp, "atom_style atomic"); lammps_command(lmp, "region box block 0 2 0 2 0 2"); lammps_command(lmp, "create_box 1 box"); - lammps_command(lmp, "variable partition universe 1 2"); + lammps_command(lmp, "variable partition universe 1 2 3 4"); EXPECT_EQ(lammps_extract_setting(lmp, "universe_size"), nprocs); EXPECT_EQ(lammps_extract_setting(lmp, "universe_rank"), me); - EXPECT_EQ(lammps_extract_setting(lmp, "world_size"), nprocs / 2); - EXPECT_EQ(lammps_extract_setting(lmp, "world_rank"), me % 2); + EXPECT_EQ(lammps_extract_setting(lmp, "world_size"), 1); + EXPECT_EQ(lammps_extract_setting(lmp, "world_rank"), 0); char *part_id = (char *)lammps_extract_variable(lmp, "partition", nullptr); - if (me < 2) { - ASSERT_THAT(part_id, StrEq("1")); - } else { - ASSERT_THAT(part_id, StrEq("2")); - } + ASSERT_THAT(part_id, StrEq(std::to_string(me+1))); lammps_close(lmp); }; diff --git a/unittest/force-styles/tests/fix-timestep-oneway.yaml b/unittest/force-styles/tests/fix-timestep-oneway.yaml index 4fcef41f5b..01131a3803 100644 --- a/unittest/force-styles/tests/fix-timestep-oneway.yaml +++ b/unittest/force-styles/tests/fix-timestep-oneway.yaml @@ -1,7 +1,7 @@ --- lammps_version: 24 Aug 2020 date_generated: Tue Sep 15 09:44:41 202 -epsilon: 2e-14 +epsilon: 4e-14 prerequisites: ! | atom full fix oneway diff --git a/unittest/force-styles/tests/mol-pair-lj_cubic.yaml b/unittest/force-styles/tests/mol-pair-lj_cubic.yaml index ad7f0666ca..6d23998ad5 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cubic.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cubic.yaml @@ -1,7 +1,7 @@ --- -lammps_version: 24 Aug 2020 -date_generated: Tue Sep 15 09:44:15 202 -epsilon: 2e-13 +lammps_version: 30 Nov 2020 +date_generated: Fri Dec 18 22:30:42 202 +epsilon: 1e-13 prerequisites: ! | atom full pair lj/cubic @@ -19,72 +19,72 @@ pair_coeff: ! | 5 5 0.015 3.1 extract: ! "" natoms: 29 -init_vdwl: 166.005266308562 +init_vdwl: 166.005266338502 init_coul: 0 init_stress: ! |2- - 4.5860676757908186e+02 4.8091912919212928e+02 1.0767204080701006e+03 -2.1005546139122362e+02 -2.9491286717936713e+00 1.6145675857120941e+02 + 4.5860676760613336e+02 4.8091912923656065e+02 1.0767204080957010e+03 -2.1005546139899724e+02 -2.9491286649299440e+00 1.6145675855773089e+02 init_forces: ! |2 - 1 9.1849370411551270e+00 7.6268937957720553e+01 6.1726872441625311e+01 - 2 2.2858712118514426e+01 1.8809274242266209e+01 -2.6905829837199740e+01 - 3 -3.2016987482543328e+01 -9.4135849525427091e+01 -3.4799279593035926e+01 - 4 -5.5341015869901478e-01 1.5206999898436971e-01 -3.9418368928369890e-01 - 5 -1.8042057425348118e-01 -3.0459951056385326e-01 8.7068483241007189e-01 - 6 -2.0038994438822397e+02 2.3344446299945159e+02 2.8487343926572851e+02 - 7 8.0909912172413883e+00 -7.8410849891085633e+01 -4.3214084684451740e+02 - 8 4.7943581255133857e+01 -2.1287511456246008e+01 1.4094503445180061e+02 - 9 1.1447552368270737e+01 1.2328709806786962e+01 5.0656476982000299e+01 - 10 1.3071496571967870e+02 -1.4589264560693914e+02 -4.4748155922123622e+01 - 11 -1.6551880116149281e-01 -4.1534332040572380e-01 -6.8284765241715795e-01 - 12 1.7721533626133388e+00 6.3456329073685158e-01 -8.2372301448028962e-01 - 13 5.6789360334118277e-01 -2.2634410312439054e-01 -9.7536738055328392e-03 - 14 -2.4337021468262635e-01 4.6659433642728905e-02 -6.1110664501270184e-01 - 15 -2.1936997101927893e-02 5.9238263972968364e-01 2.1493099548264527e-01 - 16 1.1121534968449923e+02 -7.8056927924992834e+01 -2.9249212971206231e+02 - 17 -1.1020604609843586e+02 7.6481296254913858e+01 2.9430701446263464e+02 - 18 -1.6570656719723909e-02 -2.7996966177077785e-02 2.6456326954440619e-02 - 19 7.4243353115058947e-04 6.3524893127716046e-04 1.8675586277048476e-04 - 20 -7.4243353115058947e-04 -6.3524893127716046e-04 -1.8675586277048476e-04 + 1 9.1849370400954342e+00 7.6268937960282486e+01 6.1726872440953251e+01 + 2 2.2858712118557477e+01 1.8809274242339992e+01 -2.6905829837209875e+01 + 3 -3.2016987482586380e+01 -9.4135849525500873e+01 -3.4799279593025787e+01 + 4 -5.5341015917141478e-01 1.5206999930354809e-01 -3.9418368927471259e-01 + 5 -1.8042057490785210e-01 -3.0459951013385317e-01 8.7068483295449139e-01 + 6 -2.0038994438983289e+02 2.3344446299845976e+02 2.8487343926562755e+02 + 7 8.0909912152126484e+00 -7.8410849888970290e+01 -4.3214084684693944e+02 + 8 4.7943581253518161e+01 -2.1287511458762772e+01 1.4094503445043824e+02 + 9 1.1447552368334570e+01 1.2328709806695397e+01 5.0656476981938759e+01 + 10 1.3071496571895125e+02 -1.4589264560750766e+02 -4.4748155921681018e+01 + 11 -1.6551880130741767e-01 -4.1534332035668403e-01 -6.8284765271893910e-01 + 12 1.7721533633976141e+00 6.3456328808910523e-01 -8.2372301296999062e-01 + 13 5.6789360517170440e-01 -2.2634410322310763e-01 -9.7536744156455080e-03 + 14 -2.4337021462140446e-01 4.6659433737428486e-02 -6.1110664517858926e-01 + 15 -2.1936996885396808e-02 5.9238264018028652e-01 2.1493099542010694e-01 + 16 1.1121534968641315e+02 -7.8056927926703892e+01 -2.9249212971001145e+02 + 17 -1.1020604609471918e+02 7.6481296251709367e+01 2.9430701446464468e+02 + 18 -1.6570656172995385e-02 -2.7996961634495443e-02 2.6456324093123054e-02 + 19 7.4243314752471426e-04 6.3524860303507949e-04 1.8675576627109685e-04 + 20 -7.4243314752471426e-04 -6.3524860303507949e-04 -1.8675576627109685e-04 21 -1.1415041486189516e+01 -1.3016363071591645e+01 3.6007276733401099e+01 - 22 -1.7227422089792942e+01 -4.1746638094950628e+00 -2.7029162034499002e+01 - 23 2.8642463575982458e+01 1.7191026881086707e+01 -8.9781146989020968e+00 - 24 5.8150644491939154e+00 -3.3774314134628064e+01 1.7867788752379695e+01 - 25 -2.3666545027773044e+01 3.8106021846559952e+00 -1.9896269873584632e+01 - 26 1.7843812244577855e+01 2.9960339884741117e+01 2.0167430316952100e+00 - 27 8.2825859209946024e+00 -3.6194570066818969e+01 1.4492694351988913e+01 - 28 -2.8773892796642542e+01 1.2366374307374247e+01 -1.9468877181285176e+01 - 29 2.0497044211022661e+01 2.3831279505404666e+01 4.9748677441078746e+00 -run_vdwl: 164.176727313193 + 22 -1.7227422090167991e+01 -4.1746638096674396e+00 -2.7029162034658786e+01 + 23 2.8642463576357507e+01 1.7191026881259084e+01 -8.9781146987423170e+00 + 24 5.8150644495067958e+00 -3.3774314132751599e+01 1.7867788754173183e+01 + 25 -2.3666545028156815e+01 3.8106021844353424e+00 -1.9896269873795571e+01 + 26 1.7843812244961626e+01 2.9960339884961773e+01 2.0167430319061479e+00 + 27 8.2825859198612442e+00 -3.6194570067428131e+01 1.4492694352248694e+01 + 28 -2.8773892797077618e+01 1.2366374307246014e+01 -1.9468877181493664e+01 + 29 2.0497044211457737e+01 2.3831279505532898e+01 4.9748677443163629e+00 +run_vdwl: 164.176727343123 run_coul: 0 run_stress: ! |2- - 4.5696669905868288e+02 4.7904134871830234e+02 1.0582153129928076e+03 -2.0794233475912671e+02 -2.0206236780739086e+00 1.5948889983617320e+02 + 4.5696669908578690e+02 4.7904134876274537e+02 1.0582153130184231e+03 -2.0794233476691440e+02 -2.0206236712071677e+00 1.5948889982266843e+02 run_forces: ! |2 - 1 9.2483453892728740e+00 7.5945844239974676e+01 6.1367289784738311e+01 - 2 2.2732852134554197e+01 1.8737493288759875e+01 -2.6669600068587577e+01 - 3 -3.1964986246503269e+01 -9.3734416340601200e+01 -3.4671255449722295e+01 - 4 -5.5001585694065913e-01 1.5070776418957152e-01 -3.9275226164990551e-01 - 5 -1.7969915375876547e-01 -3.0362292678324765e-01 8.6793365916706400e-01 - 6 -1.9796806590003533e+02 2.2990162655369829e+02 2.7501790745035373e+02 - 7 7.9109888171861886e+00 -7.6524641847864288e+01 -4.2032627555539375e+02 - 8 4.5992175211879918e+01 -1.9873148355500923e+01 1.3922798515578887e+02 - 9 1.1403303962503362e+01 1.2231165268449960e+01 5.0409090529604853e+01 - 10 1.3047784096912977e+02 -1.4556513307073578e+02 -4.4756481385998420e+01 - 11 -1.6346238349194633e-01 -4.0922088697872150e-01 -6.7303941060221573e-01 - 12 1.7689668574627495e+00 6.3166692380469824e-01 -8.3233385524693426e-01 - 13 5.6480104331071312e-01 -2.2395151872256039e-01 -9.5790993973179691e-03 - 14 -2.4030540495346994e-01 4.5179188229734012e-02 -6.0304369153488313e-01 - 15 -2.3220111317111478e-02 5.9408611078423390e-01 2.1676726911830960e-01 - 16 1.0963039832302356e+02 -7.7096357855469549e+01 -2.8842624961188653e+02 - 17 -1.0862142117090467e+02 7.5521002117836630e+01 2.9024023117219969e+02 - 18 -1.6565212275867415e-02 -2.7990268691876326e-02 2.6458602067006932e-02 - 19 7.1473709241289744e-04 6.1248437925700357e-04 1.7889258733572757e-04 - 20 -7.1473709241289744e-04 -6.1248437925700357e-04 -1.7889258733572757e-04 - 21 -1.1536904971247665e+01 -1.3021625993962397e+01 3.6108894191673429e+01 - 22 -1.7333879764643559e+01 -4.2314763344327275e+00 -2.7103019136756011e+01 - 23 2.8870784735891224e+01 1.7253102328395126e+01 -9.0058750549174214e+00 - 24 6.1437425316795213e+00 -3.4297207023632204e+01 1.8296742414004438e+01 - 25 -2.4276461284621075e+01 3.8560435260643189e+00 -2.0415720860228767e+01 - 26 1.8125049871956215e+01 3.0437790982988908e+01 2.1072387594169975e+00 - 27 8.4078124309265192e+00 -3.6323119973255714e+01 1.4505938075037919e+01 - 28 -2.8937319168272772e+01 1.2421253477627801e+01 -1.9540501416079319e+01 - 29 2.0535244350189391e+01 2.3904950625827414e+01 5.0332497948309163e+00 + 1 9.2483453882111135e+00 7.5945844242532630e+01 6.1367289784064468e+01 + 2 2.2732852134599622e+01 1.8737493288836678e+01 -2.6669600068598680e+01 + 3 -3.1964986246546172e+01 -9.3734416340673434e+01 -3.4671255449709932e+01 + 4 -5.5001585741105219e-01 1.5070776450703974e-01 -3.9275226164162935e-01 + 5 -1.7969915441332354e-01 -3.0362292635308630e-01 8.6793365971015668e-01 + 6 -1.9796806590165264e+02 2.2990162655270728e+02 2.7501790745023669e+02 + 7 7.9109888151537930e+00 -7.6524641845739424e+01 -4.2032627555780203e+02 + 8 4.5992175210254615e+01 -1.9873148358014198e+01 1.3922798515443196e+02 + 9 1.1403303962567387e+01 1.2231165268357689e+01 5.0409090529541324e+01 + 10 1.3047784096840681e+02 -1.4556513307130629e+02 -4.4756481385556341e+01 + 11 -1.6346238363922402e-01 -4.0922088693014880e-01 -6.7303941090576869e-01 + 12 1.7689668582530251e+00 6.3166692115802814e-01 -8.3233385373913271e-01 + 13 5.6480104514758445e-01 -2.2395151881989347e-01 -9.5791000096411370e-03 + 14 -2.4030540489083232e-01 4.5179188324553629e-02 -6.0304369170302041e-01 + 15 -2.3220111102749796e-02 5.9408611123108124e-01 2.1676726905655996e-01 + 16 1.0963039832494319e+02 -7.7096357857185083e+01 -2.8842624960984045e+02 + 17 -1.0862142116718823e+02 7.5521002114629724e+01 2.9024023117422536e+02 + 18 -1.6565211729771653e-02 -2.7990264151000276e-02 2.6458599205916818e-02 + 19 7.1473670261989000e-04 6.1248404522910675e-04 1.7889248977386897e-04 + 20 -7.1473670261989000e-04 -6.1248404522910675e-04 -1.7889248977386897e-04 + 21 -1.1536904971247079e+01 -1.3021625993961788e+01 3.6108894191671638e+01 + 22 -1.7333879765020438e+01 -4.2314763346057394e+00 -2.7103019136915339e+01 + 23 2.8870784736267517e+01 1.7253102328567529e+01 -9.0058750547563005e+00 + 24 6.1437425319896413e+00 -3.4297207021755547e+01 1.8296742415795482e+01 + 25 -2.4276461285000558e+01 3.8560435258438814e+00 -2.0415720860437396e+01 + 26 1.8125049872338021e+01 3.0437790983209180e+01 2.1072387596271578e+00 + 27 8.4078124297937791e+00 -3.6323119973862774e+01 1.4505938075297303e+01 + 28 -2.8937319168706676e+01 1.2421253477499173e+01 -1.9540501416287395e+01 + 29 2.0535244350622666e+01 2.3904950625953887e+01 5.0332497950390769e+00 ...