diff --git a/doc/lammps.1 b/doc/lammps.1 index a374e62604..58086b1fae 100644 --- a/doc/lammps.1 +++ b/doc/lammps.1 @@ -1,4 +1,4 @@ -.TH LAMMPS "1" "14 December 2021" "2021-12-14" +.TH LAMMPS "1" "7 January 2022" "2022-1-7" .SH NAME .B LAMMPS \- Molecular Dynamics Simulator. diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 9ac4fc851c..a226d73b36 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -119,6 +119,7 @@ OPT. * :doc:`granular ` * :doc:`gw ` * :doc:`gw/zbl ` + * :doc:`harmonic/cut (o) ` * :doc:`hbond/dreiding/lj (o) ` * :doc:`hbond/dreiding/morse (o) ` * :doc:`hdnnp ` diff --git a/doc/src/Developer_org.rst b/doc/src/Developer_org.rst index b5acdf5631..ce36d9a590 100644 --- a/doc/src/Developer_org.rst +++ b/doc/src/Developer_org.rst @@ -225,7 +225,7 @@ follows: commands in an input script. - The Force class computes various forces between atoms. The Pair - parent class is for non-bonded or pair-wise forces, which in LAMMPS + parent class is for non-bonded or pairwise forces, which in LAMMPS also includes many-body forces such as the Tersoff 3-body potential if those are computed by walking pairwise neighbor lists. The Bond, Angle, Dihedral, Improper parent classes are styles for bonded diff --git a/doc/src/Errors_warnings.rst b/doc/src/Errors_warnings.rst index d7a97c1507..2d588a1b77 100644 --- a/doc/src/Errors_warnings.rst +++ b/doc/src/Errors_warnings.rst @@ -416,7 +416,7 @@ This will most likely cause errors in kinetic fluctuations. not defined for the specified atom style. *Molecule has bond topology but no special bond settings* - This means the bonded atoms will not be excluded in pair-wise + This means the bonded atoms will not be excluded in pairwise interactions. *Molecule template for create_atoms has multiple molecules* diff --git a/doc/src/Run_output.rst b/doc/src/Run_output.rst index 8adfd4b293..a988be94ad 100644 --- a/doc/src/Run_output.rst +++ b/doc/src/Run_output.rst @@ -106,7 +106,7 @@ individual ranks. Here is an example output for this section: ---------- The third section above lists the number of owned atoms (Nlocal), -ghost atoms (Nghost), and pair-wise neighbors stored per processor. +ghost atoms (Nghost), and pairwise neighbors stored per processor. The max and min values give the spread of these values across processors with a 10-bin histogram showing the distribution. The total number of histogram counts is equal to the number of processors. @@ -114,7 +114,7 @@ number of histogram counts is equal to the number of processors. ---------- The last section gives aggregate statistics (across all processors) -for pair-wise neighbors and special neighbors that LAMMPS keeps track +for pairwise neighbors and special neighbors that LAMMPS keeps track of (see the :doc:`special_bonds ` command). The number of times neighbor lists were rebuilt is tallied, as is the number of potentially *dangerous* rebuilds. If atom movement triggered neighbor diff --git a/doc/src/Speed_kokkos.rst b/doc/src/Speed_kokkos.rst index 14c2ec680e..8b9b2e99af 100644 --- a/doc/src/Speed_kokkos.rst +++ b/doc/src/Speed_kokkos.rst @@ -214,7 +214,7 @@ threads/task as Nt. The product of these two values should be N, i.e. The default for the :doc:`package kokkos ` command when running on KNL is to use "half" neighbor lists and set the Newton flag to "on" for both pairwise and bonded interactions. This will typically - be best for many-body potentials. For simpler pair-wise potentials, it + be best for many-body potentials. For simpler pairwise potentials, it may be faster to use a "full" neighbor list with Newton flag to "off". Use the "-pk kokkos" :doc:`command-line switch ` to change the default :doc:`package kokkos ` options. See its page for diff --git a/doc/src/balance.rst b/doc/src/balance.rst index 5063a502bb..1d24e467d8 100644 --- a/doc/src/balance.rst +++ b/doc/src/balance.rst @@ -383,7 +383,7 @@ multiple groups, its weight is the product of the weight factors. This weight style is useful in combination with pair style :doc:`hybrid `, e.g. when combining a more costly many-body -potential with a fast pair-wise potential. It is also useful when +potential with a fast pairwise potential. It is also useful when using :doc:`run_style respa ` where some portions of the system have many bonded interactions and others none. It assumes that the computational cost for each group remains constant over time. diff --git a/doc/src/compute_pressure_cylinder.rst b/doc/src/compute_pressure_cylinder.rst index a008254540..9913ef159b 100644 --- a/doc/src/compute_pressure_cylinder.rst +++ b/doc/src/compute_pressure_cylinder.rst @@ -61,7 +61,7 @@ Restrictions This compute currently calculates the pressure tensor contributions for pair styles only (i.e. no bond, angle, dihedral, etc. contributions and in the presence of bonded interactions, the result will be incorrect -due to exclusions for special bonds) and requires pair-wise force +due to exclusions for special bonds) and requires pairwise force calculations not available for most many-body pair styles. K-space calculations are also excluded. Note that this pressure compute outputs the configurational terms only; the kinetic contribution is not included diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index be75153f6f..352f9c61bf 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -17,13 +17,14 @@ Syntax * one or more keyword/value pairs may be appended * these keywords apply to various dump styles -* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* +* keyword = *append* or *at* or *balance* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* .. parsed-literal:: *append* arg = *yes* or *no* *at* arg = N N = index of frame written upon first dump + *balance* arg = *yes* or *no* *buffer* arg = *yes* or *no* *delay* arg = Dstep Dstep = delay output until this timestep @@ -667,6 +668,14 @@ keywords are set to non-default values (i.e. the number of dump file pieces is not equal to the number of procs), then sorting cannot be performed. +In a parallel run, the per-processor dump file pieces can have +significant imbalance in number of lines of per-atom info. The *balance* +keyword determines whether the number of lines in each processor +snapshot are balanced to be nearly the same. A balance value of *no* +means no balancing will be done, while *yes* means balancing will be +performed. This balancing preserves dump sorting order. For a serial +run, this option is ignored since the output is already balanced. + .. note:: Unless it is required by the dump style, sorting dump file @@ -832,6 +841,7 @@ Default The option defaults are * append = no +* balance = no * buffer = yes for dump styles *atom*, *custom*, *loca*, and *xyz* * element = "C" for every atom type * every = whatever it was set to via the :doc:`dump ` command diff --git a/doc/src/package.rst b/doc/src/package.rst index 2cf772ea5a..437601dc60 100644 --- a/doc/src/package.rst +++ b/doc/src/package.rst @@ -460,7 +460,7 @@ using *neigh/thread* *on*, a full neighbor list must also be used. Using is turned on by default only when there are 16K atoms or less owned by an MPI rank and when using a full neighbor list. Not all KOKKOS-enabled potentials support this keyword yet, and only thread over atoms. Many -simple pair-wise potentials such as Lennard-Jones do support threading +simple pairwise potentials such as Lennard-Jones do support threading over both atoms and neighbors. The *newton* keyword sets the Newton flags for pairwise and bonded diff --git a/doc/src/pair_charmm.rst b/doc/src/pair_charmm.rst index d45ef58060..e1469ae323 100644 --- a/doc/src/pair_charmm.rst +++ b/doc/src/pair_charmm.rst @@ -119,7 +119,7 @@ name are the older, original LAMMPS implementations. They compute the LJ and Coulombic interactions with an energy switching function (esw, shown in the formula below as S(r)), which ramps the energy smoothly to zero between the inner and outer cutoff. This can cause -irregularities in pair-wise forces (due to the discontinuous second +irregularities in pairwise forces (due to the discontinuous second derivative of energy at the boundaries of the switching region), which in some cases can result in detectable artifacts in an MD simulation. diff --git a/doc/src/pair_dpd.rst b/doc/src/pair_dpd.rst index 8a8b35e50a..8c61c789a6 100644 --- a/doc/src/pair_dpd.rst +++ b/doc/src/pair_dpd.rst @@ -50,7 +50,7 @@ Style *dpd* computes a force field for dissipative particle dynamics Style *dpd/tstat* invokes a DPD thermostat on pairwise interactions, which is equivalent to the non-conservative portion of the DPD force -field. This pair-wise thermostat can be used in conjunction with any +field. This pairwise thermostat can be used in conjunction with any :doc:`pair style `, and in leiu of per-particle thermostats like :doc:`fix langevin ` or ensemble thermostats like Nose Hoover as implemented by :doc:`fix nvt `. To use diff --git a/doc/src/pair_harmonic_cut.rst b/doc/src/pair_harmonic_cut.rst new file mode 100644 index 0000000000..877d2d97ff --- /dev/null +++ b/doc/src/pair_harmonic_cut.rst @@ -0,0 +1,90 @@ +.. index:: pair_style harmonic/cut +.. index:: pair_style harmonic/cut/omp + +pair_style harmonic/cut command +=============================== + +Accelerator Variants: *harmonic/cut/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style style + +* style = *harmonic/cut* + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style harmonic/cut 2.5 + pair_coeff * * 0.2 2.0 + pair_coeff 1 1 0.5 2.5 + +Description +""""""""""" + +Style *harmonic/cut* computes pairwise repulsive-only harmonic interactions with the formula + +.. math:: + + E = k (r_c - r)^2 \qquad r < r_c + +:math:`r_c` is the cutoff. + +The following coefficients must be defined for each pair of atoms +types via the :doc:`pair_coeff ` command as in the examples +above, or in the data file or restart files read by the +:doc:`read_data ` or :doc:`read_restart ` +commands: + +* :math:`k` (energy units) +* :math:`r_c` (distance units) + +---------- + +.. include:: accel_styles.rst + +---------- + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +For atom type pairs I,J and I != J, the :math:`k` and :math:`r_c` +coefficients can be mixed. The default mix value is *geometric*. +See the "pair_modify" command for details. + +Since the potential is zero at and beyond the cutoff parameter by +construction, there is no need to support support the :doc:`pair_modify +` shift or tail options for the energy and pressure of the +pair interaction. + +These pair styles write their information to :doc:`binary restart files `, +so pair_style and pair_coeff commands do not need to be specified in an input script +that reads a restart file. + +These pair styles can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. They do not support the +*inner*, *middle*, *outer* keywords. + +---------- + +Restrictions +"""""""""""" + +The *harmonic/cut* pair style is only enabled if LAMMPS was +built with the EXTRA-PAIR package. +See the :doc:`Build package ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_coeff ` + +Default +""""""" + +none diff --git a/doc/src/pair_python.rst b/doc/src/pair_python.rst index 35e07dbd11..65e1cd1611 100644 --- a/doc/src/pair_python.rst +++ b/doc/src/pair_python.rst @@ -164,7 +164,7 @@ Following the *LJCutMelt* example, here are the two functions: .. note:: The evaluation of scripted python code will slow down the - computation pair-wise interactions quite significantly. However, this + computation pairwise interactions quite significantly. However, this can be largely worked around through using the python pair style not for the actual simulation, but to generate tabulated potentials on the fly using the :doc:`pair_write ` command. Please see below diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 4bb3c90a8d..4f1f1e733f 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -154,10 +154,10 @@ accelerated styles exist. * :doc:`coul/wolf/cs ` - Coulomb via Wolf potential with core/shell adjustments * :doc:`dpd ` - dissipative particle dynamics (DPD) * :doc:`dpd/ext ` - generalized force field for DPD -* :doc:`dpd/ext/tstat ` - pair-wise DPD thermostatting with generalized force field +* :doc:`dpd/ext/tstat ` - pairwise DPD thermostatting with generalized force field * :doc:`dpd/fdt ` - DPD for constant temperature and pressure * :doc:`dpd/fdt/energy ` - DPD for constant energy and enthalpy -* :doc:`dpd/tstat ` - pair-wise DPD thermostatting +* :doc:`dpd/tstat ` - pairwise DPD thermostatting * :doc:`dsmc ` - Direct Simulation Monte Carlo (DSMC) * :doc:`e3b ` - Explicit-three body (E3B) water model * :doc:`drip ` - Dihedral-angle-corrected registry-dependent interlayer potential (DRIP) @@ -183,6 +183,7 @@ accelerated styles exist. * :doc:`gran/hooke/history ` - granular potential without history effects * :doc:`gw ` - Gao-Weber potential * :doc:`gw/zbl ` - Gao-Weber potential with a repulsive ZBL core +* :doc:`harmonic/cut ` - repulsive-only harmonic potential * :doc:`hbond/dreiding/lj ` - DREIDING hydrogen bonding LJ potential * :doc:`hbond/dreiding/morse ` - DREIDING hydrogen bonding Morse potential * :doc:`hdnnp ` - High-dimensional neural network potential diff --git a/doc/src/pair_sw.rst b/doc/src/pair_sw.rst index d71999b2d4..d87da43b2c 100644 --- a/doc/src/pair_sw.rst +++ b/doc/src/pair_sw.rst @@ -202,7 +202,7 @@ elements are the same. Thus the two-body parameters for Si interacting with C, comes from the SiCC entry. The three-body parameters can in principle be specific to the three elements of the configuration. In the literature, however, the three-body parameters -are usually defined by simple formulas involving two sets of pair-wise +are usually defined by simple formulas involving two sets of pairwise parameters, corresponding to the ij and ik pairs, where i is the center atom. The user must ensure that the correct combining rule is used to calculate the values of the three-body parameters for diff --git a/doc/src/pair_write.rst b/doc/src/pair_write.rst index 1d88137255..b62383cb7b 100644 --- a/doc/src/pair_write.rst +++ b/doc/src/pair_write.rst @@ -76,8 +76,10 @@ must be set before this command can be invoked. Due to how the pairwise force is computed, an inner value > 0.0 must be specified even if the potential has a finite value at r = 0.0. -For EAM potentials, the pair_write command only tabulates the -pairwise portion of the potential, not the embedding portion. +The *pair_write* command can only be used for pairwise additive +interactions for which a `Pair::single()` function can be and has +been implemented. This excludes for example manybody potentials +or TIP4P coulomb styles. Related commands """""""""""""""" diff --git a/doc/src/run_style.rst b/doc/src/run_style.rst index b4d1c22113..fd63c82b90 100644 --- a/doc/src/run_style.rst +++ b/doc/src/run_style.rst @@ -89,7 +89,7 @@ in its 3d FFTs. In this scenario, splitting your P total processors into 2 subsets of processors, P1 in the first partition and P2 in the second partition, can enable your simulation to run faster. This is because the long-range forces in PPPM can be calculated at the same -time as pair-wise and bonded forces are being calculated, and the FFTs +time as pairwise and bonded forces are being calculated, and the FFTs can actually speed up when running on fewer processors. To use this style, you must define 2 partitions where P1 is a multiple diff --git a/examples/PACKAGES/electron_stopping/Si.Si.elstop b/examples/PACKAGES/electron_stopping/Si.Si.elstop index 51618149ff..1a6952e571 100644 --- a/examples/PACKAGES/electron_stopping/Si.Si.elstop +++ b/examples/PACKAGES/electron_stopping/Si.Si.elstop @@ -1,8 +1,8 @@ -# Electronic stopping for Si in Si -# Uses metal units +# Electronic stopping for Si in Si UNITS: metal DATE: 2019-02-19 CONTRIBUTOR: Risto Toijala risto.toijala@helsinki.fi # Kinetic energy in eV, stopping power in eV/A # For other atom types, add columns. +# atom type 1 # energy Si in Si 3918.2 6.541 15672.9 13.091 diff --git a/examples/PACKAGES/gle/qt-300k.A b/examples/PACKAGES/gle/qt-300k.A index 073100ad9d..31ada37b05 100644 --- a/examples/PACKAGES/gle/qt-300k.A +++ b/examples/PACKAGES/gle/qt-300k.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 3.904138445158e-4 4.681059722010e-2 3.079778738058e-2 4.428079381336e-2 5.057825203477e-2 2.591193419597e-2 1.513907125942e-2 -4.789343294190e-2 2.037551040972e-2 6.597801861779e-2 -8.528273506602e-3 -2.230839572773e-3 6.593086307870e-3 -6.653653981891e-2 -2.905096373618e-2 -6.597801861779e-2 2.086885297843e-2 1.145471984072e-2 3.111465343867e-2 1.101562087523e-2 -3.264959166486e-2 diff --git a/examples/PACKAGES/gle/qt-300k.C b/examples/PACKAGES/gle/qt-300k.C index 50d68df251..d6239ca8dc 100644 --- a/examples/PACKAGES/gle/qt-300k.C +++ b/examples/PACKAGES/gle/qt-300k.C @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 2.999985914100e+2 5.937593337000e+0 2.144376751500e+2 5.883055908000e+1 -1.119803766000e+2 -6.793381764000e+1 1.379789732400e+1 5.937593337000e+0 3.781851303000e+2 -5.794518522000e+1 -2.178772681500e+2 -1.649310659100e+2 -6.557113452000e+1 3.833830743000e+1 2.144376751500e+2 -5.794518522000e+1 7.325759985000e+2 2.188507713000e+2 -3.704586531000e+2 -3.385193865000e+1 -4.827706758000e+0 diff --git a/examples/PACKAGES/gle/smart.A b/examples/PACKAGES/gle/smart.A index a63bd881c8..8abbdea76a 100644 --- a/examples/PACKAGES/gle/smart.A +++ b/examples/PACKAGES/gle/smart.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 4.247358737218e-3 3.231404593065e-2 -9.715629522215e-3 8.199488441225e-3 7.288427565896e-3 -3.634229949055e-3 -3.294395982200e-3 4.887738278128e-2 5.978602802893e-1 -2.079222967202e-1 1.088740438247e-1 4.666954324786e-2 -1.334147134493e-2 -3.914996645811e-2 4.015863091190e-2 2.079222967202e-1 1.645970119444e-1 8.351952991453e-2 5.114740085468e-2 1.217862237677e-2 -4.506711205974e-2 diff --git a/lib/kokkos/Makefile.kokkos b/lib/kokkos/Makefile.kokkos index 7ffea5a62c..899d6e49a2 100644 --- a/lib/kokkos/Makefile.kokkos +++ b/lib/kokkos/Makefile.kokkos @@ -1224,6 +1224,9 @@ ifeq ($(KOKKOS_INTERNAL_USE_HIP), 1) KOKKOS_SRC += $(wildcard $(KOKKOS_PATH)/core/src/HIP/*.cpp) + ifeq ($(KOKKOS_INTERNAL_ENABLE_DESUL_ATOMICS), 1) + KOKKOS_SRC += $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp + endif KOKKOS_HEADERS += $(wildcard $(KOKKOS_PATH)/core/src/HIP/*.hpp) KOKKOS_CXXFLAGS+=$(KOKKOS_INTERNAL_HIP_ARCH_FLAG) diff --git a/lib/kokkos/Makefile.targets b/lib/kokkos/Makefile.targets index 93854d0cf1..c097e80fec 100644 --- a/lib/kokkos/Makefile.targets +++ b/lib/kokkos/Makefile.targets @@ -68,6 +68,8 @@ Kokkos_HIP_Instance.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/HIP/Kokkos_ $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Instance.cpp Kokkos_HIP_Locks.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Locks.cpp $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Locks.cpp +Lock_Array_HIP.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp + $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp endif ifeq ($(KOKKOS_INTERNAL_USE_PTHREADS), 1) diff --git a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp index 6964d5b41b..418aa14c68 100644 --- a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp +++ b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp @@ -694,7 +694,7 @@ std::pair CudaInternal::resize_team_scratch_space( int current_team_scratch = 0; int zero = 0; int one = 1; - while (m_team_scratch_pool[current_team_scratch].compare_exchange_weak( + while (!m_team_scratch_pool[current_team_scratch].compare_exchange_weak( zero, one, std::memory_order_release, std::memory_order_relaxed)) { current_team_scratch = (current_team_scratch + 1) % m_n_team_scratch; } diff --git a/python/lammps/core.py b/python/lammps/core.py index fcd5c76ad5..0bc8d0cb3f 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -122,6 +122,9 @@ class lammps(object): for f in os.listdir(winpath)]): lib_ext = ".dll" modpath = winpath + elif any([f.startswith('liblammps') and f.endswith('.so') + for f in os.listdir(modpath)]): + lib_ext = ".so" else: import platform if platform.system() == "Darwin": diff --git a/src/.gitignore b/src/.gitignore index 695c9a19af..a9eef4b371 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1489,6 +1489,8 @@ /pair_dpd_fdt.h /pair_dpd_fdt_energy.cpp /pair_dpd_fdt_energy.h +/pair_harmonic_cut.cpp +/pair_harmonic_cut.h /pair_lennard_mdf.cpp /pair_lennard_mdf.h /pair_lj_cut_coul_long_cs.cpp diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index 08af6ee3c9..39cb88f81b 100644 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -18,18 +18,18 @@ #include "pair_gayberne.h" -#include -#include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" -#include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "citeme.h" -#include "memory.h" +#include "comm.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include using namespace LAMMPS_NS; diff --git a/src/EXTRA-FIX/fix_electron_stopping.cpp b/src/EXTRA-FIX/fix_electron_stopping.cpp index bd54fea97d..49fcdd41bf 100644 --- a/src/EXTRA-FIX/fix_electron_stopping.cpp +++ b/src/EXTRA-FIX/fix_electron_stopping.cpp @@ -29,6 +29,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include "region.h" #include "update.h" @@ -236,49 +237,43 @@ double FixElectronStopping::compute_scalar() return SeLoss_all; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + read electron stopping parameters. only called from MPI rank 0. + format: energy then one column per atom type + read as many lines as available. + energies must be sorted in ascending order. + ---------------------------------------------------------------------- */ void FixElectronStopping::read_table(const char *file) { - char line[MAXLINE]; - - FILE *fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open stopping range table {}: {}", file, utils::getsyserror()); - const int ncol = atom->ntypes + 1; + int nlines = 0; + PotentialFileReader reader(lmp, file, "electron stopping data table"); - int l = 0; - while (true) { - if (fgets(line, MAXLINE, fp) == nullptr) break; // end of file - if (line[0] == '#') continue; // comment + try { + char *line; + double oldvalue = 0.0; - char *pch = strtok(line, " \t\n\r"); - if (pch == nullptr) continue; // blank line + while ((line = reader.next_line())) { + if (nlines >= maxlines) grow_table(); + ValueTokenizer values(line); + elstop_ranges[0][nlines] = values.next_double(); + if (elstop_ranges[0][nlines] <= oldvalue) + throw TokenizerException("energy values must be positive and in ascending order",line); - if (l >= maxlines) grow_table(); + oldvalue = elstop_ranges[0][nlines]; + for (int i = 1; i < ncol; ++i) + elstop_ranges[i][nlines] = values.next_double(); - int i = 0; - for ( ; i < ncol && pch != nullptr; i++) { - elstop_ranges[i][l] = utils::numeric(FLERR, pch,false,lmp); - pch = strtok(nullptr, " \t\n\r"); + ++nlines; } - - if (i != ncol || pch != nullptr) // too short or too long - error->one(FLERR, "fix electron/stopping: Invalid table line"); - - if (l >= 1 && elstop_ranges[0][l] <= elstop_ranges[0][l-1]) - error->one(FLERR, - "fix electron/stopping: Energies must be in ascending order"); - - l++; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing electron stopping data: {}", e.what()); } - table_entries = l; - - if (table_entries == 0) + if (nlines == 0) error->one(FLERR, "Did not find any data in electron/stopping table file"); - fclose(fp); + table_entries = nlines; } /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-FIX/fix_gle.cpp b/src/EXTRA-FIX/fix_gle.cpp index de225c6980..acb7b84fe1 100644 --- a/src/EXTRA-FIX/fix_gle.cpp +++ b/src/EXTRA-FIX/fix_gle.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: Michele Ceriotti (EPFL), Joe Morrone (Stony Brook), - Axel Kohylmeyer (Temple U) + Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ #include "fix_gle.h" @@ -24,6 +24,7 @@ #include "error.h" #include "force.h" #include "memory.h" +#include "potential_file_reader.h" #include "random_mars.h" #include "respa.h" #include "update.h" @@ -210,6 +211,8 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : S = new double[ns1sq]; TT = new double[ns1sq]; ST = new double[ns1sq]; + memset(A,0,sizeof(double)*ns1sq); + memset(C,0,sizeof(double)*ns1sq); // start temperature (t ramp) t_start = utils::numeric(FLERR,arg[4],false,lmp); @@ -221,46 +224,16 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : int seed = utils::inumeric(FLERR,arg[6],false,lmp); // LOADING A matrix - FILE *fgle = nullptr; char *fname = arg[7]; if (comm->me == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open A-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading A-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - - char line[MAXLINE],*ptr; - int n,nwords,ndone=0,eof=0; - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle A matrix"); + try { + reader.next_dvector(A, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading A-matrix: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - - nwords = utils::count_words(line); - if (nwords == 0) continue; - - ptr = strtok(line," \t\n\r\f"); - do { - A[ndone] = atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(A,ns1sq,MPI_DOUBLE,0,world); fnoneq=0; gle_every=1; gle_step=0; for (int iarg=8; iargboltz / force->mvv2e; - memset(C,0,sizeof(double)*ns1sq); for (int i=0; ime == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open C-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading C-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - ndone = eof = 0; - const double cfac = force->boltz / force->mvv2e; - - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle C matrix"); + try { + reader.next_dvector(C, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading C-matrix: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - - nwords = utils::count_words(line); - if (nwords == 0) continue; - - ptr = strtok(line," \t\n\r\f"); - do { - C[ndone] = cfac*atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(C,ns1sq,MPI_DOUBLE,0,world); + const double cfac = force->boltz / force->mvv2e; + for (int i=0; i < ns1sq; ++i) C[i] *= cfac; } #ifdef GLE_DEBUG @@ -366,12 +313,12 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : FixGLE::~FixGLE() { delete random; - delete [] A; - delete [] C; - delete [] S; - delete [] T; - delete [] TT; - delete [] ST; + delete[] A; + delete[] C; + delete[] S; + delete[] T; + delete[] TT; + delete[] ST; memory->destroy(sqrt_m); memory->destroy(gle_s); diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp index 147a63512a..8abaf97bec 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp @@ -29,7 +29,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp) {} +BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp), nn(nullptr), mm(nullptr) {} /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.h b/src/EXTRA-MOLECULE/bond_fene_nm.h index f00394c6d8..4dbf64a331 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.h +++ b/src/EXTRA-MOLECULE/bond_fene_nm.h @@ -38,7 +38,6 @@ class BondFENENM : public BondFENE { virtual void *extract(const char *, int &); protected: - double TWO_1_3; double *nn, *mm; virtual void allocate(); diff --git a/src/EXTRA-PAIR/pair_e3b.cpp b/src/EXTRA-PAIR/pair_e3b.cpp index 7e865ac6f8..3af8e754c0 100644 --- a/src/EXTRA-PAIR/pair_e3b.cpp +++ b/src/EXTRA-PAIR/pair_e3b.cpp @@ -99,7 +99,7 @@ void PairE3B::compute(int eflag, int vflag) ev_init(eflag, vflag); //clear sumExp array - memset(sumExp, 0.0, nbytes); + memset(sumExp, 0, sizeof(double)*maxID); evdwl = 0.0; pvector[0] = pvector[1] = pvector[2] = pvector[3] = 0.0; @@ -364,7 +364,6 @@ void PairE3B::allocateE3B() maxID = find_maxID(); if (!natoms) error->all(FLERR, "No atoms found"); memory->create(sumExp, maxID, "pair:sumExp"); - nbytes = sizeof(double) * maxID; } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-PAIR/pair_e3b.h b/src/EXTRA-PAIR/pair_e3b.h index f6cf916918..1f38c18293 100644 --- a/src/EXTRA-PAIR/pair_e3b.h +++ b/src/EXTRA-PAIR/pair_e3b.h @@ -50,7 +50,6 @@ class PairE3B : public Pair { int **pairO, ***pairH; // pair lists double ***exps, ****del3, ***fpair3, *sumExp; int maxID; //size of global sumExp array - size_t nbytes; //size of sumExp array in bytes int natoms; //to make sure number of atoms is constant virtual void allocate(); diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.cpp b/src/EXTRA-PAIR/pair_harmonic_cut.cpp new file mode 100644 index 0000000000..f2e1c422fe --- /dev/null +++ b/src/EXTRA-PAIR/pair_harmonic_cut.cpp @@ -0,0 +1,314 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_harmonic_cut.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCut::PairHarmonicCut(LAMMPS *lmp) : Pair(lmp), k(nullptr), cut(nullptr) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCut::~PairHarmonicCut() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(cut); + memory->destroy(cutsq); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairHarmonicCut::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, fxtmp, fytmp, fztmp; + double delx, dely, delz, rsq, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; + + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + const double fpair = delta * philj / r; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + } + + if (evflag) ev_tally(i, j, nlocal, newton_pair, philj, 0.0, fpair, delx, dely, delz); + } + } + f[i][0] += fxtmp; + f[i][1] += fytmp; + f[i][2] += fztmp; + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairHarmonicCut::allocate() +{ + allocated = 1; + int n = atom->ntypes + 1; + + memory->create(setflag, n, n, "pair:setflag"); + for (int i = 1; i < n; i++) + for (int j = i; j < n; j++) setflag[i][j] = 0; + + memory->create(k, n, n, "pair:k"); + memory->create(cut, n, n, "pair:cut"); + memory->create(cutsq, n, n, "pair:cutsq"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairHarmonicCut::settings(int narg, char **arg) +{ + if (narg > 0) error->all(FLERR, "Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHarmonicCut::coeff(int narg, char **arg) +{ + if (narg != 4) error->all(FLERR, "Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + + double k_one = utils::numeric(FLERR, arg[2], false, lmp); + double cut_one = utils::numeric(FLERR, arg[3], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { + k[i][j] = k_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairHarmonicCut::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + cut[j][i] = cut[i][j]; + k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]); + k[j][i] = k[i][j]; + } + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i, j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j], sizeof(int), 1, fp); + if (setflag[i][j]) { + fwrite(&k[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairHarmonicCut::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i, j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_restart_settings(FILE *fp) +{ + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairHarmonicCut::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); + } + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, k[i][i], cut[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) fprintf(fp, "%d %d %g %g\n", i, j, k[i][j], cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairHarmonicCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, double &fforce) +{ + if (rsq >= cutsq[itype][jtype]) { + fforce = 0.0; + return 0.0; + } + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + fforce = delta * philj / r; + return philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairHarmonicCut::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str, "k") == 0) return (void *) k; + if (strcmp(str, "cut") == 0) return (void *) cut; + return nullptr; +} diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.h b/src/EXTRA-PAIR/pair_harmonic_cut.h new file mode 100644 index 0000000000..7dfdea995a --- /dev/null +++ b/src/EXTRA-PAIR/pair_harmonic_cut.h @@ -0,0 +1,72 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 PAIR_CLASS +// clang-format off +PairStyle(harmonic/cut,PairHarmonicCut); +// clang-format on +#else + +#ifndef LMP_PAIR_HARMONIC_CUT_H +#define LMP_PAIR_HARMONIC_CUT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairHarmonicCut : public Pair { + public: + PairHarmonicCut(class LAMMPS *); + virtual ~PairHarmonicCut(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void write_data(FILE *); + void write_data_all(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double **k, **cut; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/ diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 519181be52..a201386724 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -43,7 +43,7 @@ using namespace FixConst; #ifdef __INTEL_OFFLOAD #ifndef _LMP_INTEL_OFFLOAD -#warning "Not building Intel package with Xeon Phi offload support." +#warning "Not building INTEL package with Xeon Phi offload support." #endif #endif @@ -303,8 +303,7 @@ void FixIntel::init() if (force->pair_match("^hybrid", 0) != nullptr) { _pair_hybrid_flag = 1; if (force->newton_pair != 0 && force->pair->no_virial_fdotr_compute) - error->all(FLERR, - "Intel package requires fdotr virial with newton on."); + error->all(FLERR,"INTEL package requires fdotr virial with newton on."); } else _pair_hybrid_flag = 0; @@ -325,8 +324,7 @@ void FixIntel::init() _pair_hybrid_zero = 0; if (force->newton_pair == 0) _pair_hybrid_flag = 0; if (nstyles > 1) - error->all(FLERR, - "Currently, cannot offload more than one intel style with hybrid."); + error->all(FLERR,"Currently, cannot offload more than one intel style with hybrid."); } #endif @@ -356,15 +354,12 @@ void FixIntel::init() void FixIntel::setup(int vflag) { if (neighbor->style != Neighbor::BIN) - error->all(FLERR, - "Currently, neighbor style BIN must be used with Intel package."); + error->all(FLERR,"Currently, neighbor style BIN must be used with INTEL package."); if (vflag > 3) - error->all(FLERR, - "Cannot currently get per-atom virials with Intel package."); + error->all(FLERR,"Cannot currently get per-atom virials with INTEL package."); #ifdef _LMP_INTEL_OFFLOAD if (neighbor->exclude_setting() != 0) - error->all(FLERR, - "Currently, cannot use neigh_modify exclude with Intel package offload."); + error->all(FLERR,"Currently, cannot use neigh_modify exclude with INTEL package offload."); post_force(vflag); #endif } @@ -425,14 +420,14 @@ void FixIntel::pair_init_check(const bool cdmessage) if (_offload_balance != 0.0 && comm->me == 0) { #ifndef __INTEL_COMPILER_BUILD_DATE - error->warning(FLERR, "Unknown Intel Compiler Version\n"); + error->warning(FLERR,"Unknown Intel Compiler Version\n"); #else if (__INTEL_COMPILER_BUILD_DATE != 20131008 && __INTEL_COMPILER_BUILD_DATE < 20141023) - error->warning(FLERR, "Unsupported Intel Compiler."); + error->warning(FLERR,"Unsupported Intel Compiler."); #endif #if !defined(__INTEL_COMPILER) - error->warning(FLERR, "Unsupported Intel Compiler."); + error->warning(FLERR,"Unsupported Intel Compiler."); #endif } @@ -452,7 +447,7 @@ void FixIntel::pair_init_check(const bool cdmessage) #endif int need_tag = 0; - if (atom->molecular != Atom::ATOMIC) need_tag = 1; + if (atom->molecular != Atom::ATOMIC || three_body_neighbor()) need_tag = 1; // Clear buffers used for pair style char kmode[80]; @@ -474,27 +469,28 @@ void FixIntel::pair_init_check(const bool cdmessage) #endif if (_print_pkg_info && comm->me == 0) { - if (screen) { - fprintf(screen, - "----------------------------------------------------------\n"); - if (_offload_balance != 0.0) { - fprintf(screen,"Using Intel Coprocessor with %d threads per core, ", - _offload_tpc); - fprintf(screen,"%d threads per task\n",_offload_threads); - } else { - fprintf(screen,"Using Intel Package without Coprocessor.\n"); - } - fprintf(screen,"Precision: %s\n",kmode); - if (cdmessage) { - #ifdef LMP_USE_AVXCD - fprintf(screen,"AVX512 CD Optimizations: Enabled\n"); - #else - fprintf(screen,"AVX512 CD Optimizations: Disabled\n"); - #endif - } - fprintf(screen, - "----------------------------------------------------------\n"); + utils::logmesg(lmp, "----------------------------------------------------------\n"); + if (_offload_balance != 0.0) { + utils::logmesg(lmp,"Using Intel Coprocessor with {} threads per core, " + "{} threads per task\n",_offload_tpc, _offload_threads); + } else { + utils::logmesg(lmp,"Using INTEL Package without Coprocessor.\n"); } + utils::logmesg(lmp,"Compiler: {}\n",platform::compiler_info()); + #ifdef LMP_SIMD_COMPILER + utils::logmesg(lmp,"SIMD compiler directives: Enabled\n"); + #else + utils::logmesg(lmp,"SIMD compiler directives: Disabled\n"); + #endif + utils::logmesg(lmp,"Precision: {}\n",kmode); + if (cdmessage) { + #ifdef LMP_USE_AVXCD + utils::logmesg(lmp,"AVX512 CD Optimizations: Enabled\n"); + #else + utils::logmesg(lmp,"AVX512 CD Optimizations: Disabled\n"); + #endif + } + utils::logmesg(lmp, "----------------------------------------------------------\n"); } _print_pkg_info = 0; } @@ -505,8 +501,7 @@ void FixIntel::bond_init_check() { if ((_offload_balance != 0.0) && (atom->molecular != Atom::ATOMIC) && (force->newton_pair != force->newton_bond)) - error->all(FLERR, - "INTEL package requires same setting for newton bond and non-bond."); + error->all(FLERR,"INTEL package requires same setting for newton bond and non-bond."); int intel_pair = 0; if (force->pair_match("/intel$", 0) != nullptr) @@ -517,8 +512,7 @@ void FixIntel::bond_init_check() } if (intel_pair == 0) - error->all(FLERR, "Intel styles for bond/angle/dihedral/improper " - "require intel pair style."); + error->all(FLERR,"Intel styles for bond/angle/dihedral/improper require intel pair style."); } /* ---------------------------------------------------------------------- */ @@ -534,7 +528,7 @@ void FixIntel::kspace_init_check() } if (intel_pair == 0) - error->all(FLERR, "Intel styles for kspace require intel pair style."); + error->all(FLERR,"Intel styles for kspace require intel pair style."); } /* ---------------------------------------------------------------------- */ @@ -551,7 +545,7 @@ void FixIntel::check_neighbor_intel() _offload_noghost = 0; } if (neighbor->requests[i]->skip && _offload_balance != 0.0) - error->all(FLERR, "Cannot yet use hybrid styles with Intel offload."); + error->all(FLERR,"Cannot yet use hybrid styles with Intel offload."); // avoid flagging a neighbor list as both INTEL and OPENMP if (neighbor->requests[i]->intel) @@ -786,8 +780,7 @@ void FixIntel::add_oresults(const ft * _noalias const f_in, if (f_in[1].w == 1) error->all(FLERR,"Bad matrix inversion in mldivide3"); else - error->all(FLERR, - "Sphere particles not yet supported for gayberne/intel"); + error->all(FLERR,"Sphere particles not yet supported for gayberne/intel"); } } @@ -926,7 +919,7 @@ void FixIntel::add_off_results(const ft * _noalias const f_in, int nlocal = atom->nlocal; if (neighbor->ago == 0) { if (_off_overflow_flag[LMP_OVERFLOW]) - error->one(FLERR, "Neighbor list overflow, boost neigh_modify one"); + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); _offload_nlocal = _off_overflow_flag[LMP_LOCAL_MAX] + 1; _offload_min_ghost = _off_overflow_flag[LMP_GHOST_MIN]; _offload_nghost = _off_overflow_flag[LMP_GHOST_MAX] + 1 - @@ -938,7 +931,7 @@ void FixIntel::add_off_results(const ft * _noalias const f_in, if (atom->torque) if (f_in[1].w < 0.0) - error->all(FLERR, "Bad matrix inversion in mldivide3"); + error->all(FLERR,"Bad matrix inversion in mldivide3"); add_results(f_in, ev_global, _off_results_eatom, _off_results_vatom, 1); // Load balance? @@ -1043,8 +1036,7 @@ void FixIntel::output_timing_data() { timers[TIME_OFFLOAD_PAIR]; double tt = MAX(ht,ct); if (timers[TIME_OFFLOAD_LATENCY] / tt > 0.07 && _separate_coi == 0) - error->warning(FLERR, - "Leaving a core free can improve performance for offload"); + error->warning(FLERR,"Leaving a core free can improve performance for offload"); } fprintf(_tscreen, "------------------------------------------------\n"); } @@ -1109,7 +1101,7 @@ void FixIntel::set_offload_affinity() int ppn = get_ppn(node_rank); if (ppn % _ncops != 0) - error->all(FLERR, "MPI tasks per node must be multiple of offload_cards"); + error->all(FLERR,"MPI tasks per node must be multiple of offload_cards"); ppn = ppn / _ncops; _cop = node_rank / ppn; node_rank = node_rank % ppn; @@ -1214,8 +1206,7 @@ int FixIntel::set_host_affinity(const int nomp) int subscription = nthreads * ppn; if (subscription > ncores) { if (rank == 0) - error->warning(FLERR, - "More MPI tasks/OpenMP threads than available cores"); + error->warning(FLERR,"More MPI tasks/OpenMP threads than available cores"); return 0; } if (subscription == ncores) diff --git a/src/INTEL/intel_buffers.cpp b/src/INTEL/intel_buffers.cpp index d86570b0d3..cbbc609fc0 100644 --- a/src/INTEL/intel_buffers.cpp +++ b/src/INTEL/intel_buffers.cpp @@ -207,8 +207,6 @@ void IntelBuffers::free_nmax() template void IntelBuffers::_grow_nmax(const int offload_end) { - if (lmp->atom->molecular) _need_tag = 1; - else _need_tag = 0; #ifdef _LMP_INTEL_OFFLOAD free_nmax(); int size = lmp->atom->nmax; @@ -296,9 +294,7 @@ void IntelBuffers::free_list_ptrs() /* ---------------------------------------------------------------------- */ template -void IntelBuffers::grow_data3(NeighList *list, - int *&numneighhalf, - int *&cnumneigh) +void IntelBuffers::grow_data3(NeighList *list, int *&numneighhalf, int *&cnumneigh) { const int size = list->get_maxlocal(); int list_num; @@ -321,10 +317,8 @@ void IntelBuffers::grow_data3(NeighList *list, lmp->memory->destroy(_neigh_list_ptrs[list_num].cnumneigh); lmp->memory->destroy(_neigh_list_ptrs[list_num].numneighhalf); } - lmp->memory->create(_neigh_list_ptrs[list_num].cnumneigh, size, - "_cnumneigh"); - lmp->memory->create(_neigh_list_ptrs[list_num].numneighhalf, size, - "_cnumneigh"); + lmp->memory->create(_neigh_list_ptrs[list_num].cnumneigh, size, "_cnumneigh"); + lmp->memory->create(_neigh_list_ptrs[list_num].numneighhalf, size, "_cnumneigh"); _neigh_list_ptrs[list_num].size = size; } numneighhalf = _neigh_list_ptrs[list_num].numneighhalf; @@ -334,8 +328,7 @@ void IntelBuffers::grow_data3(NeighList *list, /* ---------------------------------------------------------------------- */ template -void IntelBuffers::_grow_list_local(NeighList *list, - const int three_body, +void IntelBuffers::_grow_list_local(NeighList *list, const int three_body, const int offload_end) { free_list_local(); diff --git a/src/INTEL/nbin_intel.cpp b/src/INTEL/nbin_intel.cpp index 94f18002a0..29bc7c9ced 100644 --- a/src/INTEL/nbin_intel.cpp +++ b/src/INTEL/nbin_intel.cpp @@ -161,10 +161,8 @@ void NBinIntel::bin_atoms(IntelBuffers * buffers) { const flt_t dx = (INTEL_BIGP - bboxhi[0]); const flt_t dy = (INTEL_BIGP - bboxhi[1]); const flt_t dz = (INTEL_BIGP - bboxhi[2]); - if (dx * dx + dy * dy + dz * dz < - static_cast(neighbor->cutneighmaxsq)) - error->one(FLERR, - "Intel package expects no atoms within cutoff of {1e15,1e15,1e15}."); + if (dx * dx + dy * dy + dz * dz < static_cast(neighbor->cutneighmaxsq)) + error->one(FLERR,"INTEL package expects no atoms within cutoff of (1e15,1e15,1e15)."); } // ---------- Grow and cast/pack buffers ------------- diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 395e50006c..2fffb2353a 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -31,8 +31,7 @@ using namespace LAMMPS_NS; NPairIntel::NPairIntel(LAMMPS *lmp) : NPair(lmp) { int ifix = modify->find_fix("package_intel"); if (ifix < 0) - error->all(FLERR, - "The 'package intel' command is required for /intel styles"); + error->all(FLERR,"The 'package intel' command is required for /intel styles"); _fix = static_cast(modify->fix[ifix]); #ifdef _LMP_INTEL_OFFLOAD _cop = _fix->coprocessor_number(); @@ -659,6 +658,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ns += n2 - pack_offset - maxnbors; #ifdef LMP_INTEL_3BODY_FAST + int alln = n; n = lane; for (int u = pack_offset; u < alln; u++) { neighptr[n] = neighptr2[u]; diff --git a/src/INTEL/pair_sw_intel.cpp b/src/INTEL/pair_sw_intel.cpp index fd043ae5a1..8743fbe54c 100644 --- a/src/INTEL/pair_sw_intel.cpp +++ b/src/INTEL/pair_sw_intel.cpp @@ -1115,8 +1115,7 @@ void PairSWIntel::init_style() int ifix = modify->find_fix("package_intel"); if (ifix < 0) - error->all(FLERR, - "The 'package intel' command is required for /intel styles"); + error->all(FLERR,"The 'package intel' command is required for /intel styles"); fix = static_cast(modify->fix[ifix]); fix->pair_init_check(true); diff --git a/src/INTEL/verlet_lrt_intel.cpp b/src/INTEL/verlet_lrt_intel.cpp index 216ba98302..f867ad73e6 100644 --- a/src/INTEL/verlet_lrt_intel.cpp +++ b/src/INTEL/verlet_lrt_intel.cpp @@ -71,7 +71,7 @@ void VerletLRTIntel::init() #ifndef LMP_INTEL_USELRT error->all(FLERR, - "LRT otion for Intel package disabled at compile time"); + "LRT otion for INTEL package disabled at compile time"); #endif } diff --git a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp index 8379dc8f46..582e9a39ce 100644 --- a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp +++ b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp @@ -633,6 +633,7 @@ void FixACKS2ReaxFFKokkos::compute_h_item(int ii, int &m_fill, const template template +KOKKOS_INLINE_FUNCTION void FixACKS2ReaxFFKokkos::compute_h_team( const typename Kokkos::TeamPolicy::member_type &team, int atoms_per_team, int vector_length) const { @@ -935,6 +936,7 @@ void FixACKS2ReaxFFKokkos::compute_x_item(int ii, int &m_fill, const template template +KOKKOS_INLINE_FUNCTION void FixACKS2ReaxFFKokkos::compute_x_team( const typename Kokkos::TeamPolicy::member_type &team, int atoms_per_team, int vector_length) const { diff --git a/src/KOKKOS/pair_snap_kokkos.h b/src/KOKKOS/pair_snap_kokkos.h index bd56d87f59..e7af536782 100644 --- a/src/KOKKOS/pair_snap_kokkos.h +++ b/src/KOKKOS/pair_snap_kokkos.h @@ -85,6 +85,19 @@ public: using complex = SNAComplex; // Static team/tile sizes for device offload + +#ifdef KOKKOS_ENABLE_HIP + static constexpr int team_size_compute_neigh = 2; + static constexpr int tile_size_compute_ck = 2; + static constexpr int tile_size_pre_ui = 2; + static constexpr int team_size_compute_ui = 2; + static constexpr int tile_size_transform_ui = 2; + static constexpr int tile_size_compute_zi = 2; + static constexpr int tile_size_compute_bi = 2; + static constexpr int tile_size_transform_bi = 2; + static constexpr int tile_size_compute_yi = 2; + static constexpr int team_size_compute_fused_deidrj = 2; +#else static constexpr int team_size_compute_neigh = 4; static constexpr int tile_size_compute_ck = 4; static constexpr int tile_size_pre_ui = 4; @@ -95,6 +108,7 @@ public: static constexpr int tile_size_transform_bi = 4; static constexpr int tile_size_compute_yi = 8; static constexpr int team_size_compute_fused_deidrj = sizeof(real_type) == 4 ? 4 : 2; +#endif // Custom MDRangePolicy, Rank3, to reduce verbosity of kernel launches // This hides the Kokkos::IndexType and Kokkos::Rank<3...> diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 2259d2f829..da4a2ad658 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -67,33 +67,7 @@ MSM::MSM(LAMMPS *lmp) factors[0] = 2; MPI_Comm_rank(world,&me); - - phi1d = dphi1d = nullptr; - nmax = 0; - part2grid = nullptr; - - g_direct = nullptr; - g_direct_top = nullptr; - - v0_direct = v1_direct = v2_direct = nullptr; - v3_direct = v4_direct = v5_direct = nullptr; - - v0_direct_top = v1_direct_top = v2_direct_top = nullptr; - v3_direct_top = v4_direct_top = v5_direct_top = nullptr; - - ngrid = nullptr; - - alpha = betax = betay = betaz = nullptr; - nx_msm = ny_msm = nz_msm = nullptr; - nxlo_in = nylo_in = nzlo_in = nullptr; - nxhi_in = nyhi_in = nzhi_in = nullptr; - nxlo_out = nylo_out = nzlo_out = nullptr; - nxhi_out = nyhi_out = nzhi_out = nullptr; - delxinv = delyinv = delzinv = nullptr; - qgrid = nullptr; - egrid = nullptr; - v0grid = v1grid = v2grid = v3grid = v4grid = v5grid = nullptr; peratom_allocate_flag = 0; scalar_pressure_flag = 1; @@ -116,7 +90,7 @@ void MSM::settings(int narg, char **arg) MSM::~MSM() { - delete [] factors; + delete[] factors; deallocate(); if (peratom_allocate_flag) deallocate_peratom(); deallocate_levels(); @@ -311,6 +285,11 @@ double MSM::estimate_total_error() void MSM::setup() { + // change_box may trigger MSM::setup() before MSM::init() was called + // error out and request full initialization. + + if (!delxinv) error->all(FLERR, "MSM must be fully initialized for this operation"); + double *prd; double a = cutoff; @@ -626,15 +605,19 @@ void MSM::allocate() gcall->setup(ngcall_buf1,ngcall_buf2); npergrid = 1; + memory->destroy(gcall_buf1); + memory->destroy(gcall_buf2); memory->create(gcall_buf1,npergrid*ngcall_buf1,"msm:gcall_buf1"); memory->create(gcall_buf2,npergrid*ngcall_buf2,"msm:gcall_buf2"); // allocate memory for each grid level for (int n=0; ndestroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); memory->create3d_offset(qgrid[n],nzlo_out[n],nzhi_out[n], nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:qgrid"); + memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n], nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid"); @@ -654,6 +637,8 @@ void MSM::allocate() gc[n]->setup(ngc_buf1[n],ngc_buf2[n]); npergrid = 1; + memory->destroy(gc_buf1[n]); + memory->destroy(gc_buf2[n]); memory->create(gc_buf1[n],npergrid*ngc_buf1[n],"msm:gc_buf1"); memory->create(gc_buf2[n],npergrid*ngc_buf2[n],"msm:gc_buf2"); } else { @@ -835,11 +820,15 @@ void MSM::deallocate_levels() { if (world_levels) { for (int n=0; n < levels; ++n) { - if (qgrid[n]) + if (qgrid[n]) { memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + qgrid[n] = nullptr; + } - if (egrid[n]) + if (egrid[n]) { memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + egrid[n] = nullptr; + } if (gc) { if (gc[n]) { @@ -857,57 +846,57 @@ void MSM::deallocate_levels() } } - delete [] ngrid; + delete[] ngrid; ngrid = nullptr; memory->destroy(procneigh_levels); - delete [] world_levels; - delete [] active_flag; + delete[] world_levels; + delete[] active_flag; - delete [] gc; - delete [] gc_buf1; - delete [] gc_buf2; - delete [] ngc_buf1; - delete [] ngc_buf2; + delete[] gc; + delete[] gc_buf1; + delete[] gc_buf2; + delete[] ngc_buf1; + delete[] ngc_buf2; - delete [] alpha; - delete [] betax; - delete [] betay; - delete [] betaz; + delete[] alpha; + delete[] betax; + delete[] betay; + delete[] betaz; - delete [] nx_msm; - delete [] ny_msm; - delete [] nz_msm; + delete[] nx_msm; + delete[] ny_msm; + delete[] nz_msm; - delete [] nxlo_in; - delete [] nylo_in; - delete [] nzlo_in; + delete[] nxlo_in; + delete[] nylo_in; + delete[] nzlo_in; - delete [] nxhi_in; - delete [] nyhi_in; - delete [] nzhi_in; + delete[] nxhi_in; + delete[] nyhi_in; + delete[] nzhi_in; - delete [] nxlo_out; - delete [] nylo_out; - delete [] nzlo_out; + delete[] nxlo_out; + delete[] nylo_out; + delete[] nzlo_out; - delete [] nxhi_out; - delete [] nyhi_out; - delete [] nzhi_out; + delete[] nxhi_out; + delete[] nyhi_out; + delete[] nzhi_out; - delete [] delxinv; - delete [] delyinv; - delete [] delzinv; + delete[] delxinv; + delete[] delyinv; + delete[] delzinv; - delete [] qgrid; - delete [] egrid; + delete[] qgrid; + delete[] egrid; - delete [] v0grid; - delete [] v1grid; - delete [] v2grid; - delete [] v3grid; - delete [] v4grid; - delete [] v5grid; + delete[] v0grid; + delete[] v1grid; + delete[] v2grid; + delete[] v3grid; + delete[] v4grid; + delete[] v5grid; world_levels = nullptr; active_flag = nullptr; @@ -1060,8 +1049,7 @@ void MSM::set_grid_global() } if (flag && gridflag && me == 0) - error->warning(FLERR, - "Number of MSM mesh points changed to be a multiple of 2"); + error->warning(FLERR, "Number of MSM mesh points changed to be a multiple of 2"); // adjust Coulombic cutoff to give desired error (if requested) @@ -1087,8 +1075,7 @@ void MSM::set_grid_global() *p_cutoff = cutoff; if (me == 0) - error->warning(FLERR,"Adjusting Coulombic cutoff for " - "MSM, new cutoff = {:.8}", cutoff); + error->warning(FLERR,"Adjusting Coulombic cutoff for MSM, new cutoff = {:.8}", cutoff); } if (triclinic == 0) { diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 06d90e5d0e..d531233a3a 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -197,11 +197,9 @@ void PPPM::init() error->all(FLERR,"Must redefine kspace_style after changing to triclinic box"); if (domain->triclinic && differentiation_flag == 1) - error->all(FLERR,"Cannot (yet) use PPPM with triclinic box " - "and kspace_modify diff ad"); + error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and kspace_modify diff ad"); if (domain->triclinic && slabflag) - error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and " - "slab correction"); + error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and slab correction"); if (domain->dimension == 2) error->all(FLERR,"Cannot use PPPM with 2d simulation"); diff --git a/src/MAKE/MACHINES/Makefile.white b/src/MAKE/MACHINES/Makefile.crusher_kokkos similarity index 85% rename from src/MAKE/MACHINES/Makefile.white rename to src/MAKE/MACHINES/Makefile.crusher_kokkos index cb101998b3..7dc1447d4e 100644 --- a/src/MAKE/MACHINES/Makefile.white +++ b/src/MAKE/MACHINES/Makefile.crusher_kokkos @@ -1,4 +1,4 @@ -# white = KOKKOS/CUDA package, OpenMPI with nvcc compiler, Pascal GPU, Power8 CPU +# crusher_kokkos = KOKKOS/HIP, AMD MI250X GPU and AMD EPYC 7A53 "Optimized 3rd Gen EPYC" CPU, Cray MPICH, hipcc compiler SHELL = /bin/sh @@ -7,13 +7,13 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) -export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper -CC = mpicxx -CCFLAGS = -g -O3 + +CC = hipcc +CCFLAGS = -g -O3 -munsafe-fp-atomics -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M -LINK = mpicxx +LINK = hipcc LINKFLAGS = -g -O3 LIB = SIZE = size @@ -21,8 +21,8 @@ SIZE = size ARCHIVE = ar ARFLAGS = -rc SHLIBFLAGS = -shared -KOKKOS_DEVICES = Cuda -KOKKOS_ARCH = Pascal60,Power8 +KOKKOS_DEVICES = HIP +KOKKOS_ARCH = Zen3,Vega90A # --------------------------------------------------------------------- # LAMMPS-specific settings, all OPTIONAL @@ -43,9 +43,9 @@ LMP_INC = -DLAMMPS_GZIP # PATH = path for MPI library # LIB = name of MPI library -MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 +MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 -I${MPICH_DIR}/include MPI_PATH = -MPI_LIB = +MPI_LIB = -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa # FFT library # see discussion in Section 3.5.2 of manual @@ -69,11 +69,6 @@ JPG_INC = JPG_PATH = JPG_LIB = -# library for loading shared objects (defaults to -ldl, should be empty on Windows) -# uncomment to change the default - -# override DYN_LIB = - # --------------------------------------------------------------------- # build rules and dependencies # do not edit this section @@ -83,7 +78,7 @@ include Makefile.package EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC) EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH) -EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) $(DYN_LIB) +EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS) EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS) @@ -123,6 +118,6 @@ depend : fastdep.exe $(SRC) @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1 fastdep.exe: ../DEPEND/fastdep.c - gcc -O -o $@ $< + cc -O -o $@ $< sinclude .depend diff --git a/src/MAKE/MACHINES/Makefile.spock_kokkos b/src/MAKE/MACHINES/Makefile.spock_kokkos new file mode 100644 index 0000000000..a85ebb3039 --- /dev/null +++ b/src/MAKE/MACHINES/Makefile.spock_kokkos @@ -0,0 +1,123 @@ +# spock_kokkos = KOKKOS/HIP, AMD MI100 GPU and AMD EPYC 7662 "Rome" CPU, Cray MPICH, hipcc compiler + +SHELL = /bin/sh + +# --------------------------------------------------------------------- +# compiler/linker settings +# specify flags and libraries needed for your compiler + +KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) + +CC = hipcc +CCFLAGS = -g -O3 -DNDEBUG +SHFLAGS = -fPIC +DEPFLAGS = -M + +LINK = hipcc +LINKFLAGS = -g -O3 +LIB = +SIZE = size + +ARCHIVE = ar +ARFLAGS = -rc +SHLIBFLAGS = -shared +KOKKOS_DEVICES = HIP +KOKKOS_ARCH = Zen2,Vega908 + +# --------------------------------------------------------------------- +# LAMMPS-specific settings, all OPTIONAL +# specify settings for LAMMPS features you will use +# if you change any -D setting, do full re-compile after "make clean" + +# LAMMPS ifdef settings +# see possible settings in Section 3.5 of the manual + +LMP_INC = -DLAMMPS_GZIP + +# MPI library +# see discussion in Section 3.4 of the manual +# MPI wrapper compiler/linker can provide this info +# can point to dummy MPI library in src/STUBS as in Makefile.serial +# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts +# INC = path for mpi.h, MPI compiler settings +# PATH = path for MPI library +# LIB = name of MPI library + +MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 -I${MPICH_DIR}/include +MPI_PATH = +MPI_LIB = -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa + +# FFT library +# see discussion in Section 3.5.2 of manual +# can be left blank to use provided KISS FFT library +# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings +# PATH = path for FFT library +# LIB = name of FFT library + +FFT_INC = +FFT_PATH = +FFT_LIB = + +# JPEG and/or PNG library +# see discussion in Section 3.5.4 of manual +# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC +# INC = path(s) for jpeglib.h and/or png.h +# PATH = path(s) for JPEG library and/or PNG library +# LIB = name(s) of JPEG library and/or PNG library + +JPG_INC = +JPG_PATH = +JPG_LIB = + +# --------------------------------------------------------------------- +# build rules and dependencies +# do not edit this section + +include Makefile.package.settings +include Makefile.package + +EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC) +EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH) +EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) +EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS) +EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS) + +# Path to src files + +vpath %.cpp .. +vpath %.h .. + +# Link target + +$(EXE): main.o $(LMPLIB) $(EXTRA_LINK_DEPENDS) + $(LINK) $(LINKFLAGS) main.o $(EXTRA_PATH) $(LMPLINK) $(EXTRA_LIB) $(LIB) -o $@ + $(SIZE) $@ + +# Library targets + +$(ARLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + @rm -f ../$(ARLIB) + $(ARCHIVE) $(ARFLAGS) ../$(ARLIB) $(OBJ) + @rm -f $(ARLIB) + @ln -s ../$(ARLIB) $(ARLIB) + +$(SHLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + $(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o ../$(SHLIB) \ + $(OBJ) $(EXTRA_LIB) $(LIB) + @rm -f $(SHLIB) + @ln -s ../$(SHLIB) $(SHLIB) + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $< + +# Individual dependencies + +depend : fastdep.exe $(SRC) + @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1 + +fastdep.exe: ../DEPEND/fastdep.c + cc -O -o $@ $< + +sinclude .depend diff --git a/src/MAKE/MACHINES/Makefile.summit_kokkos b/src/MAKE/MACHINES/Makefile.summit_kokkos index 87f8c75da2..f22b27cc74 100644 --- a/src/MAKE/MACHINES/Makefile.summit_kokkos +++ b/src/MAKE/MACHINES/Makefile.summit_kokkos @@ -9,7 +9,7 @@ SHELL = /bin/sh KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) CC = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi b/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi index cb3ef0e442..42a8236c7c 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi +++ b/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi @@ -10,7 +10,7 @@ KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) export MPICH_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only b/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only index 6d5e8d779e..0adb53eef0 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only +++ b/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_omp b/src/MAKE/OPTIONS/Makefile.kokkos_omp index e505da8ae2..82144652dd 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_omp +++ b/src/MAKE/OPTIONS/Makefile.kokkos_omp @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_phi b/src/MAKE/OPTIONS/Makefile.kokkos_phi index b825ad691a..9d5691251c 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_phi +++ b/src/MAKE/OPTIONS/Makefile.kokkos_phi @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp index 46d439b879..dc1ec360c6 100644 --- a/src/MANYBODY/pair_eam_cd.cpp +++ b/src/MANYBODY/pair_eam_cd.cpp @@ -499,10 +499,13 @@ void PairEAMCD::read_h_coeff(char *filename) // Seek to end of file, read last part into a buffer and // then skip over lines in buffer until reaching the end. - platform::fseek(fptr, platform::END_OF_FILE); - platform::fseek(fptr, platform::ftell(fptr) - MAXLINE); + if ( (platform::fseek(fptr, platform::END_OF_FILE) < 0) + || (platform::fseek(fptr, platform::ftell(fptr) - MAXLINE) < 0)) + error->one(FLERR,"Failure to seek to end-of-file for reading h(x) coeffs: {}", + utils::getsyserror()); + char *buf = new char[MAXLINE+1]; - fread(buf, 1, MAXLINE, fptr); + utils::sfread(FLERR, buf, 1, MAXLINE, fptr, filename, error); buf[MAXLINE] = '\0'; // must 0-terminate buffer for string processing Tokenizer lines(buf, "\n"); delete[] buf; diff --git a/src/MANYBODY/pair_edip.cpp b/src/MANYBODY/pair_edip.cpp index af2314b8b9..611b1af455 100644 --- a/src/MANYBODY/pair_edip.cpp +++ b/src/MANYBODY/pair_edip.cpp @@ -31,6 +31,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -368,7 +369,7 @@ void PairEDIP::compute(int eflag, int vflag) directorCos_ik_z = invR_ik * dr_ik[2]; cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + - directorCos_ij_z * directorCos_ik_z; + directorCos_ij_z * directorCos_ik_z; cosTetaDiff = cosTeta + tauFunction; cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; @@ -376,33 +377,33 @@ void PairEDIP::compute(int eflag, int vflag) expMinusQFunctionCosTetaDiffCosTetaDiff = exp(-qFunctionCosTetaDiffCosTetaDiff); potentia3B_factor = lambda * - ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + - eta * qFunctionCosTetaDiffCosTetaDiff); + ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + + eta * qFunctionCosTetaDiffCosTetaDiff); exp3B_ik = preExp3B_ij[neighbor_k]; exp3BDerived_ik = preExp3BDerived_ij[neighbor_k]; forceMod3B_factor1_ij = -exp3BDerived_ij * exp3B_ik * potentia3B_factor; forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff * - (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); + (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); forceMod3B_factor1_ik = -exp3BDerived_ik * exp3B_ij * potentia3B_factor; forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); forceModCoord += (forceMod3B_factor2 * (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); @@ -765,136 +766,92 @@ double PairEDIP::init_one(int i, int j) void PairEDIP::read_file(char *file) { - int params_per_line = 20; - char **words = new char *[params_per_line + 1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file, lmp, nullptr); - if (fp == nullptr) - error->one(FLERR, "Cannot open EDIP potential file {}: {}", file, utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n, nwords, ielement, jelement, kelement; - char line[MAXLINE], *ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - ptr = fgets(line, MAXLINE, fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else - n = strlen(line) + 1; - } - MPI_Bcast(&eof, 1, MPI_INT, 0, world); - if (eof) break; - MPI_Bcast(&n, 1, MPI_INT, 0, world); - MPI_Bcast(line, n, MPI_CHAR, 0, world); + for (ielement = 0; ielement < nelements; ielement++) + if (iname == elements[ielement]) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (jname == elements[jelement]) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (kname == elements[kelement]) break; + if (kelement == nelements) continue; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n], MAXLINE - n, fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else - n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - MPI_Bcast(&eof, 1, MPI_INT, 0, world); - if (eof) break; - MPI_Bcast(&n, 1, MPI_INT, 0, world); - MPI_Bcast(line, n, MPI_CHAR, 0, world); - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) error->all(FLERR, "Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line, " \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr, " \t\n\r\f"))) continue; - - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next entry in file - - for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0], elements[ielement]) == 0) break; - if (ielement == nelements) continue; - for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1], elements[jelement]) == 0) break; - if (jelement == nelements) continue; - for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2], elements[kelement]) == 0) break; - if (kelement == nelements) continue; - - // load up parameter settings and error check their values - - if (nparams == maxparam) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params"); - - // make certain all addional allocated storage is initialized - // to avoid false positives when checking with valgrind - - memset(params + nparams, 0, DELTA * sizeof(Param)); - } - - params[nparams].ielement = ielement; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamm = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || - params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || - params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamm < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || - params[nparams].sigma < 0.0) - error->all(FLERR, "Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete[] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -946,7 +903,7 @@ void PairEDIP::setup_params() cutoffC = params[0].cutoffC; sigma = params[0].sigma; lambda = params[0].lambda; - gamm = params[0].gamm; + gamm = params[0].gamma; eta = params[0].eta; Q0 = params[0].Q0; mu = params[0].mu; diff --git a/src/MANYBODY/pair_edip.h b/src/MANYBODY/pair_edip.h index 5812768d55..a0d1e45613 100644 --- a/src/MANYBODY/pair_edip.h +++ b/src/MANYBODY/pair_edip.h @@ -34,14 +34,23 @@ class PairEDIP : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { - double A, B; - double cutoffA, cutoffC, cutsq; - double alpha, beta; - double eta, gamm, lambda, mu, rho, sigma, Q0; - double u1, u2, u3, u4; - int ielement, jelement, kelement; + 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, dummy; // dummy added for better alignment }; double *preInvR_ij; diff --git a/src/MANYBODY/pair_edip_multi.cpp b/src/MANYBODY/pair_edip_multi.cpp index 8017fa4f8e..ffbc135a34 100644 --- a/src/MANYBODY/pair_edip_multi.cpp +++ b/src/MANYBODY/pair_edip_multi.cpp @@ -30,6 +30,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -118,8 +119,8 @@ void PairEDIPMulti::compute(int eflag, int vflag) double costheta; double dpairZ,dtripleZ; - // eflag != 0 means compute energy contributions in this step - // vflag != 0 means compute virial contributions in this step + // eflag != 0 means compute energy contributions in this step + // vflag != 0 means compute virial contributions in this step evdwl = 0.0; ev_init(eflag,vflag); @@ -155,39 +156,40 @@ void PairEDIPMulti::compute(int eflag, int vflag) // pre-loop to compute environment coordination f(Z) for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + j = jlist[jj]; + j &= NEIGHMASK; - double delx, dely, delz, r_ij; + double delx, dely, delz, r_ij; - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp; - r_ij = delx * delx + dely * dely + delz * delz; + delx = x[j][0] - xtmp; + dely = x[j][1] - ytmp; + delz = x[j][2] - ztmp; + r_ij = delx * delx + dely * dely + delz * delz; - jtype = map[type[j]]; - ijparam = elem3param[itype][jtype][jtype]; - if (r_ij > params[ijparam].cutsq) continue; + jtype = map[type[j]]; + const Param ¶m = params[elem3param[itype][jtype][jtype]]; - r_ij = sqrt(r_ij); + if (r_ij > param.cutsq) continue; - // zeta and its derivative dZ/dr + r_ij = sqrt(r_ij); - if (r_ij < params[ijparam].cutoffC) zeta_i += 1.0; - else { - double f, fdr; - edip_fc(r_ij, ¶ms[ijparam], f, fdr); - zeta_i += f; - dzetair = -fdr / r_ij; + // zeta and its derivative dZ/dr - preForceCoord_counter=numForceCoordPairs*5; - preForceCoord[preForceCoord_counter+0]=dzetair; - preForceCoord[preForceCoord_counter+1]=delx; - preForceCoord[preForceCoord_counter+2]=dely; - preForceCoord[preForceCoord_counter+3]=delz; - preForceCoord[preForceCoord_counter+4]=j; - numForceCoordPairs++; - } + if (r_ij < param.cutoffC) zeta_i += 1.0; + else { + double f, fdr; + edip_fc(r_ij, param, f, fdr); + zeta_i += f; + dzetair = -fdr / r_ij; + + preForceCoord_counter=numForceCoordPairs*5; + preForceCoord[preForceCoord_counter+0]=dzetair; + preForceCoord[preForceCoord_counter+1]=delx; + preForceCoord[preForceCoord_counter+2]=dely; + preForceCoord[preForceCoord_counter+3]=delz; + preForceCoord[preForceCoord_counter+4]=j; + numForceCoordPairs++; + } } // two-body interactions @@ -217,7 +219,7 @@ void PairEDIPMulti::compute(int eflag, int vflag) // already considered in constructing the potential double fdr, fdZ; - edip_pair(r_ij, zeta_i, ¶ms[ijparam], evdwl, fdr, fdZ); + edip_pair(r_ij, zeta_i, params[ijparam], evdwl, fdr, fdZ); fpair = -fdr / r_ij; dpairZ += fdZ; @@ -234,102 +236,102 @@ void PairEDIPMulti::compute(int eflag, int vflag) // three-body Forces for (kk = jj + 1; kk < jnum; kk++) { - double dr_ik[3], r_ik, f_ik[3]; + double dr_ik[3], r_ik, f_ik[3]; - k = jlist[kk]; - k &= NEIGHMASK; - ktype = map[type[k]]; - ikparam = elem3param[itype][ktype][ktype]; - ijkparam = elem3param[itype][jtype][ktype]; + k = jlist[kk]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem3param[itype][ktype][ktype]; + ijkparam = elem3param[itype][jtype][ktype]; - dr_ik[0] = x[k][0] - xtmp; - dr_ik[1] = x[k][1] - ytmp; - dr_ik[2] = x[k][2] - ztmp; - r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + dr_ik[0] = x[k][0] - xtmp; + dr_ik[1] = x[k][1] - ytmp; + dr_ik[2] = x[k][2] - ztmp; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; - if (r_ik > params[ikparam].cutsq) continue; + if (r_ik > params[ikparam].cutsq) continue; - r_ik = sqrt(r_ik); + r_ik = sqrt(r_ik); - costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; + costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; - double v1, v2, v3, v4, v5, v6, v7; + double v1, v2, v3, v4, v5, v6, v7; - edip_fcut3(r_ij, ¶ms[ijparam], v1, v2); - edip_fcut3(r_ik, ¶ms[ikparam], v3, v4); - edip_h(costheta, zeta_i, ¶ms[ijkparam], v5, v6, v7); + edip_fcut3(r_ij, params[ijparam], v1, v2); + edip_fcut3(r_ik, params[ikparam], v3, v4); + edip_h(costheta, zeta_i, params[ijkparam], v5, v6, v7); - // potential energy and forces - evdwl = v1 * v3 * v5; - dtripleZ += v1 * v3 * v7; + // potential energy and forces + evdwl = v1 * v3 * v5; + dtripleZ += v1 * v3 * v7; - double dri[3], drj[3], drk[3]; - double dhl, dfr; + double dri[3], drj[3], drk[3]; + double dhl, dfr; - dhl = v1 * v3 * v6; + dhl = v1 * v3 * v6; - costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); + costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); - f_ij[0] = -dhl * drj[0]; - f_ij[1] = -dhl * drj[1]; - f_ij[2] = -dhl * drj[2]; - f_ik[0] = -dhl * drk[0]; - f_ik[1] = -dhl * drk[1]; - f_ik[2] = -dhl * drk[2]; + f_ij[0] = -dhl * drj[0]; + f_ij[1] = -dhl * drj[1]; + f_ij[2] = -dhl * drj[2]; + f_ik[0] = -dhl * drk[0]; + f_ik[1] = -dhl * drk[1]; + f_ik[2] = -dhl * drk[2]; - dfr = v2 * v3 * v5; - fpair = -dfr / r_ij; + dfr = v2 * v3 * v5; + fpair = -dfr / r_ij; - f_ij[0] += fpair * dr_ij[0]; - f_ij[1] += fpair * dr_ij[1]; - f_ij[2] += fpair * dr_ij[2]; + f_ij[0] += fpair * dr_ij[0]; + f_ij[1] += fpair * dr_ij[1]; + f_ij[2] += fpair * dr_ij[2]; - dfr = v1 * v4 * v5; - fpair = -dfr / r_ik; + dfr = v1 * v4 * v5; + fpair = -dfr / r_ik; - f_ik[0] += fpair * dr_ik[0]; - f_ik[1] += fpair * dr_ik[1]; - f_ik[2] += fpair * dr_ik[2]; + f_ik[0] += fpair * dr_ik[0]; + f_ik[1] += fpair * dr_ik[1]; + f_ik[2] += fpair * dr_ik[2]; - f[j][0] += f_ij[0]; - f[j][1] += f_ij[1]; - f[j][2] += f_ij[2]; + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; - f[k][0] += f_ik[0]; - f[k][1] += f_ik[1]; - f[k][2] += f_ik[2]; + f[k][0] += f_ik[0]; + f[k][1] += f_ik[1]; + f[k][2] += f_ik[2]; - f[i][0] -= f_ij[0] + f_ik[0]; - f[i][1] -= f_ij[1] + f_ik[1]; - f[i][2] -= f_ij[2] + f_ik[2]; + f[i][0] -= f_ij[0] + f_ik[0]; + f[i][1] -= f_ij[1] + f_ik[1]; + f[i][2] -= f_ij[2] + f_ik[2]; - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); } } // forces due to environment coordination f(Z) for (int idx = 0; idx < numForceCoordPairs; idx++) { - double delx, dely, delz; + double delx, dely, delz; - preForceCoord_counter = idx * 5; - dzetair = preForceCoord[preForceCoord_counter+0]; - delx = preForceCoord[preForceCoord_counter+1]; - dely = preForceCoord[preForceCoord_counter+2]; - delz = preForceCoord[preForceCoord_counter+3]; - j = static_cast (preForceCoord[preForceCoord_counter+4]); + preForceCoord_counter = idx * 5; + dzetair = preForceCoord[preForceCoord_counter+0]; + delx = preForceCoord[preForceCoord_counter+1]; + dely = preForceCoord[preForceCoord_counter+2]; + delz = preForceCoord[preForceCoord_counter+3]; + j = static_cast (preForceCoord[preForceCoord_counter+4]); - dzetair *= (dpairZ + dtripleZ); + dzetair *= (dpairZ + dtripleZ); - f[j][0] += dzetair * delx; - f[j][1] += dzetair * dely; - f[j][2] += dzetair * delz; + f[j][0] += dzetair * delx; + f[j][1] += dzetair * dely; + f[j][2] += dzetair * delz; - f[i][0] -= dzetair * delx; - f[i][1] -= dzetair * dely; - f[i][2] -= dzetair * delz; + f[i][0] -= dzetair * delx; + f[i][1] -= dzetair * dely; + f[i][2] -= dzetair * delz; - evdwl = 0.0; - if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); + evdwl = 0.0; + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); } } @@ -342,13 +344,13 @@ double sqr(double x) } //pair Vij, partial derivatives dVij(r,Z)/dr and dVij(r,Z)/dZ -void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, +void PairEDIPMulti::edip_pair(double r, double z, const Param ¶m, double &eng, double &fdr, double &fZ) { - double A = param->A; - double B = param->B; - double rho = param->rho; - double beta = param->beta; + double A = param.A; + double B = param.B; + double rho = param.rho; + double beta = param.beta; double v1,v2,v3,v4; v1 = pow(B / r, rho); @@ -361,11 +363,11 @@ void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, } //function fc(r) in calculating coordination Z and derivative fc'(r) -void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fc(double r, const Param ¶m, double &f, double &fdr) { - double a = param->cutoffA; - double c = param->cutoffC; - double alpha = param->alpha; + double a = param.cutoffA; + double c = param.cutoffC; + double alpha = param.alpha; double x; double v1, v2; @@ -392,10 +394,10 @@ void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) } //cut-off function for Vij and its derivative fcut2'(r) -void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut2(double r, const Param ¶m, double &f, double &fdr) { - double sigma = param->sigma; - double a = param->cutoffA; + double sigma = param.sigma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -411,12 +413,12 @@ void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) } //function tau(Z) and its derivative tau'(Z) -void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) +void PairEDIPMulti::edip_tau(double z, const Param ¶m, double &f, double &fdZ) { - double u1 = param->u1; - double u2 = param->u2; - double u3 = param->u3; - double u4 = param->u4; + double u1 = param.u1; + double u2 = param.u2; + double u3 = param.u3; + double u4 = param.u4; double v1, v2; v1 = exp(-u4 * z); @@ -427,13 +429,13 @@ void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) } //function h(l,Z) and its partial derivatives dh(l,Z)/dl and dh(l,Z)/dZ -void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, +void PairEDIPMulti::edip_h(double l, double z, const Param ¶m, double &f, double &fdl, double &fdZ) { - double lambda = param->lambda; - double eta = param->eta; - double Q0 = param->Q0; - double mu = param->mu; + double lambda = param.lambda; + double eta = param.eta; + double Q0 = param.Q0; + double mu = param.mu; double Q, QdZ, Tau, TaudZ; double u2, du2l, du2Z; double v1, v2, v3; @@ -464,10 +466,10 @@ void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, } //cut-off function for Vijk and its derivative fcut3'(r) -void PairEDIPMulti::edip_fcut3(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut3(double r, const Param ¶m, double &f, double &fdr) { - double gamma = param->gamma; - double a = param->cutoffA; + double gamma = param.gamma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -578,137 +580,92 @@ double PairEDIPMulti::init_one(int i, int j) void PairEDIPMulti::read_file(char *file) { - int params_per_line = 20; - char **words = new char*[params_per_line+1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open EDIP potential file {}: {}",file,utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + for (ielement = 0; ielement < nelements; ielement++) + if (iname == elements[ielement]) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (jname == elements[jelement]) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (kname == elements[kelement]) break; + if (kelement == nelements) continue; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line," \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue; - - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next entry in file - - for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0],elements[ielement]) == 0) break; - if (ielement == nelements) continue; - for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1],elements[jelement]) == 0) break; - if (jelement == nelements) continue; - for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2],elements[kelement]) == 0) break; - if (kelement == nelements) continue; - - // load up parameter settings and error check their values - - if (nparams == maxparam) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); - - // make certain all addional allocated storage is initialized - // to avoid false positives when checking with valgrind - - memset(params + nparams, 0, DELTA*sizeof(Param)); - } - - params[nparams].ielement = ielement; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamma = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || - params[nparams].cutoffA < 0.0 || params[nparams].cutoffC < 0.0 || - params[nparams].alpha < 0.0 || params[nparams].beta < 0.0 || - params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || - params[nparams].rho < 0.0 || params[nparams].sigma < 0.0) - error->all(FLERR,"Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete [] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -753,5 +710,4 @@ void PairEDIPMulti::setup() rtmp = sqrt(params[m].cutsq); if (rtmp > cutmax) cutmax = rtmp; } - } diff --git a/src/MANYBODY/pair_edip_multi.h b/src/MANYBODY/pair_edip_multi.h index 3ee7347a56..1070956bc2 100644 --- a/src/MANYBODY/pair_edip_multi.h +++ b/src/MANYBODY/pair_edip_multi.h @@ -34,6 +34,8 @@ class PairEDIPMulti : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { double A, B; // coefficients for pair interaction I-J @@ -48,7 +50,7 @@ class PairEDIPMulti : public Pair { double mu, Q0; // coefficients for function Q(Z) double u1, u2, u3, u4; // coefficients for function tau(Z) double cutsq; - int ielement, jelement, kelement; + int ielement, jelement, kelement, dummy; // dummy added for better alignment }; double *preForceCoord; @@ -63,12 +65,12 @@ class PairEDIPMulti : public Pair { void read_file(char *); void setup(); - void edip_pair(double, double, Param *, double &, double &, double &); - void edip_fc(double, Param *, double &, double &); - void edip_fcut2(double, Param *, double &, double &); - void edip_tau(double, Param *, double &, double &); - void edip_h(double, double, Param *, double &, double &, double &); - void edip_fcut3(double, Param *, double &, double &); + void edip_pair(double, double, const Param &, double &, double &, double &); + void edip_fc(double, const Param &, double &, double &); + void edip_fcut2(double, const Param &, double &, double &); + void edip_tau(double, const Param &, double &, double &); + void edip_h(double, double, const Param &, double &, double &, double &); + void edip_fcut3(double, const Param &, double &, double &); }; } // namespace LAMMPS_NS diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index f21c3bffca..484adc7436 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -431,8 +431,7 @@ void PairTersoff::read_file(char *file) // transparently convert units for supported conversions int unit_convert = reader.get_unit_convert(); - double conversion_factor = utils::get_conversion_factor(utils::ENERGY, - unit_convert); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); while ((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 5f7ffcbbe5..c0978f51cb 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; using namespace FixConst; static const char cite_fix_bond_swap[] = - "neighbor multi command:\n\n" + "fix bond/swap command:\n\n" "@Article{Auhl03,\n" " author = {R. Auhl, R. Everaers, G. S. Grest, K. Kremer, S. J. Plimpton},\n" " title = {Equilibration of long chain polymer melts in computer simulations},\n" diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index d1051ecb11..f5c82ed5e6 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -18,19 +18,28 @@ #include "pair_list.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" -#include "error.h" +#include "text_file_reader.h" +#include "tokenizer.h" +#include +#include +#include +#include using namespace LAMMPS_NS; -static const char * const stylename[] = { - "none", "harmonic", "morse", "lj126", nullptr +enum { NONE = 0, HARM, MORSE, LJ126 }; + +static std::map stylename = { + { "none", NONE }, + { "harmonic", HARM }, + { "morse", MORSE }, + { "lj126", LJ126 } }; // fast power function for integer exponent > 0 @@ -55,7 +64,6 @@ PairList::PairList(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; respa_enable = 0; cut_global = 0.0; - style = nullptr; params = nullptr; check_flag = 1; } @@ -66,7 +74,6 @@ PairList::~PairList() { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(style); memory->destroy(params); } @@ -90,7 +97,7 @@ void PairList::compute(int eflag, int vflag) int pc = 0; for (int n=0; n < npairs; ++n) { - const list_parm_t &par = params[n]; + const list_param &par = params[n]; i = atom->map(par.id1); j = atom->map(par.id2); @@ -122,34 +129,34 @@ void PairList::compute(int eflag, int vflag) if (rsq < par.cutsq) { const double r2inv = 1.0/rsq; - if (style[n] == HARM) { + if (par.style == HARM) { const double r = sqrt(rsq); - const double dr = par.parm.harm.r0 - r; - fpair = 2.0*par.parm.harm.k*dr/r; + const double dr = par.param.harm.r0 - r; + fpair = 2.0*par.param.harm.k*dr/r; if (eflag_either) - epair = par.parm.harm.k*dr*dr - par.offset; + epair = par.param.harm.k*dr*dr - par.offset; - } else if (style[n] == MORSE) { + } else if (par.style == MORSE) { const double r = sqrt(rsq); - const double dr = par.parm.morse.r0 - r; - const double dexp = exp(par.parm.morse.alpha * dr); - fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha + const double dr = par.param.morse.r0 - r; + const double dexp = exp(par.param.morse.alpha * dr); + fpair = 2.0*par.param.morse.d0*par.param.morse.alpha * (dexp*dexp - dexp) / r; if (eflag_either) - epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; + epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = r2inv*r2inv*r2inv; - const double sig6 = mypow(par.parm.lj126.sigma,6); - fpair = 24.0*par.parm.lj126.epsilon*r6inv + const double sig6 = mypow(par.param.lj126.sigma,6); + fpair = 24.0*par.param.lj126.epsilon*r6inv * (2.0*sig6*sig6*r6inv - sig6) * r2inv; if (eflag_either) - epair = 4.0*par.parm.lj126.epsilon*r6inv + epair = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6) - par.offset; } @@ -210,130 +217,74 @@ void PairList::settings(int narg, char **arg) if (strcmp(arg[2],"check") == 0) check_flag = 1; } - FILE *fp = utils::open_potential(arg[0],lmp,nullptr); - char line[1024]; - if (fp == nullptr) - error->all(FLERR,"Cannot open pair list file"); + std::vector mystyles; + std::vector myparams; - // count lines in file for upper limit of storage needed - int num = 1; - while (fgets(line,1024,fp)) ++num; - rewind(fp); - memory->create(style,num,"pair_list:style"); - memory->create(params,num,"pair_list:params"); - - // now read and parse pair list file for real - npairs = 0; - char *ptr; - tagint id1, id2; - int nharm=0, nmorse=0, nlj126=0; - - while (fgets(line,1024,fp)) { - ptr = strtok(line," \t\n\r\f"); - - // skip empty lines - if (!ptr) continue; - - // skip comment lines starting with # - if (*ptr == '#') continue; - - // get atom ids of pair - id1 = ATOTAGINT(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - id2 = ATOTAGINT(ptr); - - // get potential type - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - - style[npairs] = NONE; - list_parm_t &par = params[npairs]; - par.id1 = id1; - par.id2 = id2; - - // harmonic potential - if (strcmp(ptr,stylename[HARM]) == 0) { - style[npairs] = HARM; - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.k = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nharm; - - // morse potential - } else if (strcmp(ptr,stylename[MORSE]) == 0) { - style[npairs] = MORSE; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.d0 = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.alpha = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nmorse; - - } else if (strcmp(ptr,stylename[LJ126]) == 0) { - // 12-6 lj potential - style[npairs] = LJ126; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.epsilon = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.sigma = utils::numeric(FLERR,ptr,false,lmp); - - ++nlj126; - - } else { - error->all(FLERR,"Unknown pair list potential style"); - } - - // optional cutoff parameter. if not specified use global value - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr != nullptr) && (*ptr != '#')) { - double cut = utils::numeric(FLERR,ptr,false,lmp); - par.cutsq = cut*cut; - } else { - par.cutsq = cut_global*cut_global; - } - - // count complete entry - ++npairs; - } - fclose(fp); - - // informative output + // read and parse potential file only on MPI rank 0. if (comm->me == 0) { - if (screen) - fprintf(screen,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); - if (logfile) - fprintf(logfile,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); + int nharm, nmorse, nlj126, nskipped; + FILE *fp = utils::open_potential(arg[0],lmp,nullptr); + TextFileReader reader(fp,"pair list coeffs"); + npairs = nharm = nmorse = nlj126 = nskipped = 0; + + try { + char *line; + while ((line = reader.next_line())) { + ValueTokenizer values(line); + list_param oneparam; + oneparam.offset = 0.0; + oneparam.id1 = values.next_tagint(); + oneparam.id2 = values.next_tagint(); + oneparam.style = stylename[values.next_string()]; + ++npairs; + + switch (oneparam.style) { + + case HARM: + oneparam.param.harm.k = values.next_double(); + oneparam.param.harm.r0 = values.next_double(); + ++nharm; + break; + + case MORSE: + oneparam.param.morse.d0 = values.next_double(); + oneparam.param.morse.alpha = values.next_double(); + oneparam.param.morse.r0 = values.next_double(); + ++nmorse; + break; + + case LJ126: + oneparam.param.lj126.epsilon = values.next_double(); + oneparam.param.lj126.sigma = values.next_double(); + ++nlj126; + break; + + case NONE: // fallthrough + error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", + utils::trim(line)); + ++nskipped; + break; + } + if (values.has_next()) + oneparam.cutsq = values.next_double(); + else + oneparam.cutsq = cut_global*cut_global; + + myparams.push_back(oneparam); + } + } catch (std::exception &e) { + error->one(FLERR,"Error reading pair list coeffs file: {}", e.what()); + } + utils::logmesg(lmp, "Read {} ({}/{}/{}) interacting pair lines from {}. " + "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, arg[0], nskipped); + + memory->create(params,npairs,"pair_list:params"); + memcpy(params, myparams.data(),npairs*sizeof(list_param)); + fclose(fp); } + MPI_Bcast(&npairs, 1, MPI_INT, 0, world); + if (comm->me != 0) memory->create(params,npairs,"pair_list:params"); + MPI_Bcast(params, npairs*sizeof(list_param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- @@ -374,21 +325,21 @@ void PairList::init_style() if (offset_flag) { for (int n=0; n < npairs; ++n) { - list_parm_t &par = params[n]; + list_param &par = params[n]; - if (style[n] == HARM) { - const double dr = sqrt(par.cutsq) - par.parm.harm.r0; - par.offset = par.parm.harm.k*dr*dr; + if (par.style == HARM) { + const double dr = sqrt(par.cutsq) - par.param.harm.r0; + par.offset = par.param.harm.k*dr*dr; - } else if (style[n] == MORSE) { - const double dr = par.parm.morse.r0 - sqrt(par.cutsq); - const double dexp = exp(par.parm.morse.alpha * dr); - par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp); + } else if (par.style == MORSE) { + const double dr = par.param.morse.r0 - sqrt(par.cutsq); + const double dexp = exp(par.param.morse.alpha * dr); + par.offset = par.param.morse.d0 * (dexp*dexp - 2.0*dexp); - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = par.cutsq*par.cutsq*par.cutsq; - const double sig6 = mypow(par.parm.lj126.sigma,6); - par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + const double sig6 = mypow(par.param.lj126.sigma,6); + par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); } } } @@ -411,7 +362,7 @@ double PairList::init_one(int, int) double PairList::memory_usage() { double bytes = (double)npairs * sizeof(int); - bytes += (double)npairs * sizeof(list_parm_t); + bytes += (double)npairs * sizeof(list_param); const int n = atom->ntypes+1; bytes += (double)n*(n*sizeof(int) + sizeof(int *)); bytes += (double)n*(n*sizeof(double) + sizeof(double *)); diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 60fd3a7a24..b7161c4f2c 100644 --- a/src/MISC/pair_list.h +++ b/src/MISC/pair_list.h @@ -39,8 +39,6 @@ class PairList : public Pair { protected: void allocate(); - enum { NONE = 0, HARM, MORSE, LJ126 }; - // potential specific parameters struct harm_p { double k, r0; @@ -52,23 +50,23 @@ class PairList : public Pair { double epsilon, sigma; }; - union parm_u { - struct harm_p harm; - struct morse_p morse; - struct lj126_p lj126; + union param_u { + harm_p harm; + morse_p morse; + lj126_p lj126; }; - typedef struct { + struct list_param { + int style; // potential style indicator tagint id1, id2; // global atom ids double cutsq; // cutoff**2 for this pair double offset; // energy offset - union parm_u parm; // parameters for style - } list_parm_t; + union param_u param; // parameters for style + }; protected: double cut_global; // global cutoff distance - int *style; // list of styles for pair interactions - list_parm_t *params; // lisf of pair interaction parameters + list_param *params; // lisf of pair interaction parameters int npairs; // # of atom pairs in global list int check_flag; // 1 if checking for missing pairs }; diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index b3eea79064..30fd460dcc 100644 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -18,24 +18,25 @@ #include "bond_quartic.h" -#include #include "atom.h" -#include "neighbor.h" #include "comm.h" -#include "force.h" -#include "pair.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "neighbor.h" +#include "pair.h" +#include + +#include "math_const.h" using namespace LAMMPS_NS; +using namespace MathConst; /* ---------------------------------------------------------------------- */ -BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp) -{ - TWO_1_3 = pow(2.0,(1.0/3.0)); -} +BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp), k(nullptr), + b1(nullptr), b2(nullptr), rc(nullptr), u0(nullptr) {} /* ---------------------------------------------------------------------- */ @@ -120,7 +121,7 @@ void BondQuartic::compute(int eflag, int vflag) rb = dr - b2[type]; fbond = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; fbond += 48.0*sr6*(sr6-0.5)/rsq; @@ -128,7 +129,7 @@ void BondQuartic::compute(int eflag, int vflag) if (eflag) { ebond = k[type]*r2*ra*rb + u0[type]; - if (rsq < TWO_1_3) ebond += 4.0*sr6*(sr6-1.0) + 1.0; + if (rsq < MY_CUBEROOT2) ebond += 4.0*sr6*(sr6-1.0) + 1.0; } // apply force to each of 2 atoms @@ -336,7 +337,7 @@ double BondQuartic::single(int type, double rsq, int i, int j, eng += k[type]*r2*ra*rb + u0[type]; fforce = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; eng += 4.0*sr6*(sr6-1.0) + 1.0; diff --git a/src/MOLECULE/bond_quartic.h b/src/MOLECULE/bond_quartic.h index 35ec705849..3973477e68 100644 --- a/src/MOLECULE/bond_quartic.h +++ b/src/MOLECULE/bond_quartic.h @@ -38,7 +38,6 @@ class BondQuartic : public Bond { double single(int, double, int, int, double &); protected: - double TWO_1_3; double *k, *b1, *b2, *rc, *u0; void allocate(); diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index a763c5d14c..34f54ef26c 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -37,6 +37,7 @@ #include "force.h" #include "math_const.h" #include "memory.h" +#include "potential_file_reader.h" #include "respa.h" #include "tokenizer.h" #include "update.h" @@ -634,150 +635,23 @@ double FixCMAP::compute_scalar() void FixCMAP::read_grid_map(char *cmapfile) { - char linebuf[MAXLINE]; - char *chunk,*line; - int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter; - - FILE *fp = nullptr; if (comm->me == 0) { - fp = utils::open_potential(cmapfile,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open fix cmap file {}: {}", - cmapfile, utils::getsyserror()); + try { + memset(&cmapgrid[0][0][0], 0, 6*CMAPDIM*CMAPDIM*sizeof(double)); + PotentialFileReader reader(lmp, cmapfile, "cmap grid"); - } + // there are six maps in this order. + // alanine, alanine-proline, proline, proline-proline, glycine, glycine-proline. + // read as one big blob of numbers while ignoring comments - for (int ix1 = 0; ix1 < 6; ix1++) - for (int ix2 = 0; ix2 < CMAPDIM; ix2++) - for (int ix3 = 0; ix3 < CMAPDIM; ix3++) - cmapgrid[ix1][ix2][ix3] = 0.0; + reader.next_dvector(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM); - counter = 0; - i1 = i2 = i3 = i4 = i5 = i6 = 0; - j1 = j2 = j3 = j4 = j5 = j6 = 0; - - int done = 0; - - while (!done) { - // only read on rank 0 and broadcast to all other ranks - if (comm->me == 0) - done = (fgets(linebuf,MAXLINE,fp) == nullptr); - - MPI_Bcast(&done,1,MPI_INT,0,world); - if (done) continue; - - MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world); - - // remove leading whitespace - line = linebuf; - while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line; - - // skip if empty line or comment - if (!line || *line =='\n' || *line == '\0' || *line == '#') continue; - - // read in the cmap grid point values - // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors - // will occur if content of the file "cmap.data" is altered - // - // Reading order of the maps: - // 1. Alanine map - // 2. Alanine before proline map - // 3. Proline map - // 4. Two adjacent prolines map - // 5. Glycine map - // 6. Glycine before proline map - - chunk = strtok(line, " \r\n"); - while (chunk != nullptr) { - - // alanine map - - if (counter < CMAPDIM*CMAPDIM) { - cmapgrid[0][i1][j1] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j1++; - if (j1 == CMAPDIM) { - j1 = 0; - i1++; - } - counter++; - } - - // alanine-proline map - - else if (counter >= CMAPDIM*CMAPDIM && - counter < 2*CMAPDIM*CMAPDIM) { - cmapgrid[1][i2][j2]= atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j2++; - if (j2 == CMAPDIM) { - j2 = 0; - i2++; - } - counter++; - } - - // proline map - - else if (counter >= 2*CMAPDIM*CMAPDIM && - counter < 3*CMAPDIM*CMAPDIM) { - cmapgrid[2][i3][j3] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j3++; - if (j3 == CMAPDIM) { - j3 = 0; - i3++; - } - counter++; - } - - // 2 adjacent prolines map - - else if (counter >= 3*CMAPDIM*CMAPDIM && - counter < 4*CMAPDIM*CMAPDIM) { - cmapgrid[3][i4][j4] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j4++; - if (j4 == CMAPDIM) { - j4 = 0; - i4++; - } - counter++; - } - - // glycine map - - else if (counter >= 4*CMAPDIM*CMAPDIM && - counter < 5*CMAPDIM*CMAPDIM) { - cmapgrid[4][i5][j5] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j5++; - if (j5 == CMAPDIM) { - j5 = 0; - i5++; - } - counter++; - } - - // glycine-proline map - - else if (counter >= 5*CMAPDIM*CMAPDIM && - counter < 6*CMAPDIM*CMAPDIM) { - cmapgrid[5][i6][j6] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j6++; - if (j6 == CMAPDIM) { - j6 = 0; - i6++; - } - counter++; - } - - else break; + } catch (std::exception &e) { + error->one(FLERR,"Error reading CMAP potential file: {}", e.what()); } } - if (comm->me == 0) fclose(fp); + MPI_Bcast(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ diff --git a/src/Makefile b/src/Makefile index 2d651c4986..d23058447a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -439,7 +439,7 @@ clean: clean-all: rm -rf Obj_* - rm style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h + rm -f style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h clean-%: @if [ $@ = "clean-serial" ]; \ diff --git a/src/OPENMP/bond_quartic_omp.cpp b/src/OPENMP/bond_quartic_omp.cpp index 00e8357644..427293dd3e 100644 --- a/src/OPENMP/bond_quartic_omp.cpp +++ b/src/OPENMP/bond_quartic_omp.cpp @@ -22,13 +22,15 @@ #include "comm.h" #include "force.h" #include "neighbor.h" - #include "pair.h" #include #include "suffix.h" +#include "math_const.h" + using namespace LAMMPS_NS; +using namespace MathConst; /* ---------------------------------------------------------------------- */ @@ -143,7 +145,7 @@ void BondQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) rb = dr - b2[type]; fbond = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; fbond += 48.0*sr6*(sr6-0.5)/rsq; @@ -151,7 +153,7 @@ void BondQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) if (EFLAG) { ebond = k[type]*r2*ra*rb + u0[type]; - if (rsq < TWO_1_3) ebond += 4.0*sr6*(sr6-1.0) + 1.0; + if (rsq < MY_CUBEROOT2) ebond += 4.0*sr6*(sr6-1.0) + 1.0; } // apply force to each of 2 atoms diff --git a/src/OPENMP/pair_harmonic_cut_omp.cpp b/src/OPENMP/pair_harmonic_cut_omp.cpp new file mode 100644 index 0000000000..6fcdc50a63 --- /dev/null +++ b/src/OPENMP/pair_harmonic_cut_omp.cpp @@ -0,0 +1,152 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_harmonic_cut_omp.h" + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neigh_list.h" +#include "suffix.h" + +#include "omp_compat.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCutOMP::PairHarmonicCutOMP(LAMMPS *lmp) : PairHarmonicCut(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void PairHarmonicCutOMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) + eval<1, 1, 1>(ifrom, ito, thr); + else + eval<1, 1, 0>(ifrom, ito, thr); + } else { + if (force->newton_pair) + eval<1, 0, 1>(ifrom, ito, thr); + else + eval<1, 0, 0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) + eval<0, 0, 1>(ifrom, ito, thr); + else + eval<0, 0, 0>(ifrom, ito, thr); + } + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void PairHarmonicCutOMP::eval(int iifrom, int iito, ThrData *const thr) +{ + const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0]; + const int *_noalias const type = atom->type; + const double *_noalias const special_lj = force->special_lj; + const int *_noalias const ilist = list->ilist; + const int *_noalias const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; + + double xtmp, ytmp, ztmp, delx, dely, delz, fxtmp, fytmp, fztmp; + double rsq, factor_lj; + + const int nlocal = atom->nlocal; + int j, jj, jnum, jtype; + + // loop over neighbors of my atoms + + for (int ii = iifrom; ii < iito; ++ii) { + const int i = ilist[ii]; + const int itype = type[i]; + const int *_noalias const jlist = firstneigh[i]; + const double *_noalias const cutsqi = cutsq[itype]; + + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsqi[jtype]) { + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + const double fpair = delta * philj / r; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx * fpair; + f[j].y -= dely * fpair; + f[j].z -= delz * fpair; + } + + if (EVFLAG) + ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, philj, 0.0, fpair, delx, dely, delz, thr); + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairHarmonicCutOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairHarmonicCut::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_harmonic_cut_omp.h b/src/OPENMP/pair_harmonic_cut_omp.h new file mode 100644 index 0000000000..4f40d2cf2a --- /dev/null +++ b/src/OPENMP/pair_harmonic_cut_omp.h @@ -0,0 +1,48 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(harmonic/cut/omp,PairHarmonicCutOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HARMONIC_CUT_OMP_H +#define LMP_PAIR_HARMONIC_CUT_OMP_H + +#include "pair_harmonic_cut.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairHarmonicCutOMP : public PairHarmonicCut, public ThrOMP { + + public: + PairHarmonicCutOMP(class LAMMPS *); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/REAXFF/fix_acks2_reaxff.cpp b/src/REAXFF/fix_acks2_reaxff.cpp index b6789b1b2e..21cce5b55a 100644 --- a/src/REAXFF/fix_acks2_reaxff.cpp +++ b/src/REAXFF/fix_acks2_reaxff.cpp @@ -429,7 +429,7 @@ void FixACKS2ReaxFF::compute_X() double **x = atom->x; int *mask = atom->mask; - memset(X_diag,0.0,atom->nmax*sizeof(double)); + memset(X_diag,0,atom->nmax*sizeof(double)); // fill in the X matrix m_fill = 0; diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index fa8fa79e00..08a4e21fe5 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -1084,7 +1084,7 @@ void FixQEqReaxFF::vector_add(double* dest, double c, double* v, int k) void FixQEqReaxFF::get_chi_field() { - memset(&chi_field[0],0.0,atom->nmax*sizeof(double)); + memset(&chi_field[0],0,atom->nmax*sizeof(double)); if (!efield) return; const double * const *x = (const double * const *)atom->x; diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 4f44cc7c64..bed7fabb20 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -246,8 +246,10 @@ FixReaxFFSpecies::~FixReaxFFSpecies() if (posflag && multipos_opened) fclose(pos); } - modify->delete_compute(fmt::format("SPECATOM_{}",id)); - modify->delete_fix(fmt::format("SPECBOND_{}",id)); + try { + modify->delete_compute(fmt::format("SPECATOM_{}",id)); + modify->delete_fix(fmt::format("SPECBOND_{}",id)); + } catch (std::exception &) {} } /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 14742155db..cbaef073da 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -220,7 +220,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"infile") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); - delete [] inpfile; + delete[] inpfile; inpfile = utils::strdup(arg[iarg+1]); restart_file = 1; reinitflag = 0; @@ -327,7 +327,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; else { allremap = 0; - delete [] id_dilate; + delete[] id_dilate; id_dilate = utils::strdup(arg[iarg+1]); int idilate = group->find(id_dilate); if (idilate == -1) @@ -354,7 +354,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"gravity") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); - delete [] id_gravity; + delete[] id_gravity; id_gravity = utils::strdup(arg[iarg+1]); iarg += 2; @@ -395,7 +395,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : double time1 = platform::walltime(); create_bodies(bodyID); - if (customflag) delete [] bodyID; + if (customflag) delete[] bodyID; if (comm->me == 0) utils::logmesg(lmp," create bodies CPU = {:.3f} seconds\n", @@ -501,9 +501,9 @@ FixRigidSmall::~FixRigidSmall() memory->destroy(dorient); delete random; - delete [] inpfile; - delete [] id_dilate; - delete [] id_gravity; + delete[] inpfile; + delete[] id_dilate; + delete[] id_gravity; memory->destroy(langextra); memory->destroy(mass_body); @@ -2578,8 +2578,8 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody) if (me == 0) fclose(fp); - delete [] buffer; - delete [] values; + delete[] buffer; + delete[] values; } /* ---------------------------------------------------------------------- diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index e289c179d9..27076ccac3 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -84,8 +84,9 @@ class FixRigidSmall : public Fix { double maxextent; // furthest distance from body owner to body atom struct Body { - double mass; // total mass of body int natoms; // total number of atoms in body + int ilocal; // index of owning atom + double mass; // total mass of body double xcm[3]; // COM position double xgc[3]; // geometric center position double vcm[3]; // COM velocity @@ -97,12 +98,12 @@ class FixRigidSmall : public Fix { double ey_space[3]; double ez_space[3]; double xgc_body[3]; // geometric center relative to xcm in body coords - double angmom[3]; // space-frame angular momentum of body - double omega[3]; // space-frame omega of body - double conjqm[4]; // conjugate quaternion momentum - imageint image; // image flags of xcm - int remapflag[4]; // PBC remap flags - int ilocal; // index of owning atom + double angmom[3]; // space-frame angular momentum of body + double omega[3]; // space-frame omega of body + double conjqm[4]; // conjugate quaternion momentum + int remapflag[4]; // PBC remap flags + imageint image; // image flags of xcm + imageint dummy; // dummy entry for better alignment }; Body *body; // list of rigid bodies, owned and ghost diff --git a/src/angle_hybrid.cpp b/src/angle_hybrid.cpp index 12bc9aa97a..41f46ae878 100644 --- a/src/angle_hybrid.cpp +++ b/src/angle_hybrid.cpp @@ -240,7 +240,7 @@ void AngleHybrid::settings(int narg, char **arg) error->all(FLERR, "Angle style hybrid cannot have none as an argument"); styles[nstyles] = force->new_angle(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/bond_hybrid.cpp b/src/bond_hybrid.cpp index dad1010a1e..3139d9ecd8 100644 --- a/src/bond_hybrid.cpp +++ b/src/bond_hybrid.cpp @@ -242,7 +242,7 @@ void BondHybrid::settings(int narg, char **arg) if (strncmp(arg[i], "quartic", 7) == 0) has_quartic = m; styles[nstyles] = force->new_bond(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index ff9acdc4ef..b98f6eec33 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -277,10 +277,13 @@ int ComputePairLocal::compute_pairs(int flag) break; case DX: ptr[n] = delx*directionCorrection; + break; case DY: ptr[n] = dely*directionCorrection; + break; case DZ: ptr[n] = delz*directionCorrection; + break; case ENG: ptr[n] = eng; break; diff --git a/src/dihedral_hybrid.cpp b/src/dihedral_hybrid.cpp index 24a18e6783..0f72f999c9 100644 --- a/src/dihedral_hybrid.cpp +++ b/src/dihedral_hybrid.cpp @@ -241,7 +241,7 @@ void DihedralHybrid::settings(int narg, char **arg) error->all(FLERR, "Dihedral style hybrid cannot have none as an argument"); styles[nstyles] = force->new_dihedral(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/dump.cpp b/src/dump.cpp index e5652c210c..181a56ff46 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -75,6 +75,7 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp) clearstep = 0; sort_flag = 0; + balance_flag = 0; append_flag = 0; buffer_allow = 0; buffer_flag = 0; @@ -417,6 +418,7 @@ void Dump::write() if (sort_flag && sortcol == 0) pack(ids); else pack(nullptr); if (sort_flag) sort(); + if (balance_flag) balance(); // write timestep header // for multiproc, @@ -888,6 +890,160 @@ int Dump::bufcompare_reverse(const int i, const int j, void *ptr) #endif +/* ---------------------------------------------------------------------- + parallel load balance of buf across all procs + must come after sort +------------------------------------------------------------------------- */ + +void Dump::balance() +{ + bigint *proc_offsets,*proc_new_offsets; + memory->create(proc_offsets,nprocs+1,"dump:proc_offsets"); + memory->create(proc_new_offsets,nprocs+1,"dump:proc_new_offsets"); + + // compute atom offset for this proc + + bigint offset; + bigint bnme = nme; + MPI_Scan(&bnme,&offset,1,MPI_LMP_BIGINT,MPI_SUM,world); + + // gather atom offsets for all procs + + MPI_Allgather(&offset,1,MPI_LMP_BIGINT,&proc_offsets[1],1,MPI_LMP_BIGINT,world); + + proc_offsets[0] = 0; + + // how many atoms should I own after balance + + int nme_balance = static_cast(ntotal/nprocs); + + // include remainder atoms on first procs + + int remainder = ntotal % nprocs; + if (me < remainder) nme_balance += 1; + + // compute new atom offset for this proc + + bigint offset_balance; + bigint bnme_balance = nme_balance; + MPI_Scan(&bnme_balance,&offset_balance,1,MPI_LMP_BIGINT,MPI_SUM,world); + + // gather new atom offsets for all procs + + MPI_Allgather(&offset_balance,1,MPI_LMP_BIGINT,&proc_new_offsets[1],1,MPI_LMP_BIGINT,world); + + proc_new_offsets[0] = 0; + + // reset buf size to largest of any post-balance nme values + // this insures proc 0 can receive everyone's info + // cannot shrink buf to nme_balance, must use previous maxbuf value + + int nmax; + MPI_Allreduce(&nme_balance,&nmax,1,MPI_INT,MPI_MAX,world); + if (nmax > maxbuf) maxbuf = nmax; + + // allocate a second buffer for balanced data + + double* buf_balance; + memory->create(buf_balance,maxbuf*size_one,"dump:buf_balance"); + + // compute from which procs I am receiving atoms + // post recvs first + + int nswap = 0; + MPI_Request *request = new MPI_Request[nprocs]; + int procstart = 0; + int iproc = me; + int iproc_prev; + + for (int i = 0; i < nme_balance; i++) { + + // find which proc this atom belongs to + + while (proc_new_offsets[me] + i < proc_offsets[iproc]) iproc--; + while (proc_new_offsets[me] + i > proc_offsets[iproc+1]-1) iproc++; + + if (i != 0 && (iproc != iproc_prev || i == nme_balance - 1)) { + + // finished with proc + + int procrecv = iproc; + if (iproc != iproc_prev) procrecv = iproc_prev; + + int procnrecv = i - procstart + 1; + if (iproc != iproc_prev) procnrecv--; + + // post receive for this proc + + if (iproc_prev != me) + MPI_Irecv(&buf_balance[procstart*size_one],procnrecv*size_one,MPI_DOUBLE, + procrecv,0,world,&request[nswap++]); + + procstart = i; + } + + iproc_prev = iproc; + } + + // compute which atoms I am sending and to which procs + + procstart = 0; + iproc = me; + for (int i = 0; i < nme; i++) { + + // find which proc this atom should belong to + + while (proc_offsets[me] + i < proc_new_offsets[iproc]) iproc--; + while (proc_offsets[me] + i > proc_new_offsets[iproc+1] - 1) iproc++; + + if (i != 0 && (iproc != iproc_prev || i == nme - 1)) { + + // finished with proc + + int procsend = iproc; + if (iproc != iproc_prev) procsend = iproc_prev; + + int procnsend = i - procstart + 1; + if (iproc != iproc_prev) procnsend--; + + // send for this proc + + if (iproc_prev != me) { + MPI_Send(&buf[procstart*size_one],procnsend*size_one,MPI_DOUBLE,procsend,0,world); + } else { + + // sending to self, copy buffers + + int offset_me = proc_offsets[me] - proc_new_offsets[me]; + memcpy(&buf_balance[(offset_me + procstart)*size_one],&buf[procstart*size_one],procnsend*size_one*sizeof(double)); + } + + procstart = i; + } + + iproc_prev = iproc; + } + + // wait for all recvs + + for (int n = 0; n < nswap; n++) + MPI_Wait(&request[n],MPI_STATUS_IGNORE); + + nme = nme_balance; + + // swap buffers + + double *tmp = buf; + buf = buf_balance; + + // cleanup + + memory->destroy(tmp); + memory->destroy(proc_offsets); + memory->destroy(proc_new_offsets); + delete [] request; +} + /* ---------------------------------------------------------------------- process params common to all dumps here if unknown param, call modify_param specific to the dump @@ -1108,6 +1264,12 @@ void Dump::modify_params(int narg, char **arg) } iarg += 2; + } else if (strcmp(arg[iarg],"balance") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + if (nprocs > 1) + balance_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + iarg += 2; + } else if (strcmp(arg[iarg],"time") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); time_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); diff --git a/src/dump.h b/src/dump.h index 35da154d7c..83edda3810 100644 --- a/src/dump.h +++ b/src/dump.h @@ -19,6 +19,7 @@ namespace LAMMPS_NS { class Dump : protected Pointers { + friend class Output; public: char *id; // user-defined name of Dump char *style; // style of Dump @@ -66,6 +67,7 @@ class Dump : protected Pointers { int header_flag; // 0 = item, 2 = xyz int flush_flag; // 0 if no flush, 1 if flush every dump int sort_flag; // 1 if sorted output + int balance_flag; // 1 if load balanced output int append_flag; // 1 if open file in append mode, 0 if not int buffer_allow; // 1 if style allows for buffer_flag, 0 if not int buffer_flag; // 1 if buffer output as one big string, 0 if not @@ -160,6 +162,7 @@ class Dump : protected Pointers { static int bufcompare(const int, const int, void *); static int bufcompare_reverse(const int, const int, void *); #endif + void balance(); }; } // namespace LAMMPS_NS diff --git a/src/dump_image.cpp b/src/dump_image.cpp index c073d152f8..9aecb58120 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -1200,11 +1200,11 @@ int DumpImage::modify_param(int narg, char **arg) utils::bounds(FLERR,arg[1],1,atom->ntypes,nlo,nhi,error); // get list of colors + // assign colors in round-robin fashion to types + auto colors = Tokenizer(arg[2],"/").as_vector(); const int ncolors = colors.size(); - // assign colors in round-robin fashion to types - int m = 0; for (int i = nlo; i <= nhi; i++) { colortype[i] = image->color2rgb(colors[m%ncolors].c_str()); @@ -1249,32 +1249,19 @@ int DumpImage::modify_param(int narg, char **arg) int nlo,nhi; utils::bounds(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,error); - // ptrs = list of ncount colornames separated by '/' + // process list of ncount colornames separated by '/' + // assign colors in round-robin fashion to bond types - int ncount = 1; - char *nextptr; - char *ptr = arg[2]; - while ((nextptr = strchr(ptr,'/'))) { - ptr = nextptr + 1; - ncount++; - } - char **ptrs = new char*[ncount+1]; - ncount = 0; - ptrs[ncount++] = strtok(arg[2],"/"); - while ((ptrs[ncount++] = strtok(nullptr,"/"))); - ncount--; - - // assign each of ncount colors in round-robin fashion to types + auto colors = Tokenizer(arg[2],"/").as_vector(); + const int ncolors = colors.size(); int m = 0; for (int i = nlo; i <= nhi; i++) { - bcolortype[i] = image->color2rgb(ptrs[m%ncount]); + bcolortype[i] = image->color2rgb(colors[m%ncolors].c_str()); if (bcolortype[i] == nullptr) error->all(FLERR,"Invalid color in dump_modify command"); m++; } - - delete [] ptrs; return 3; } diff --git a/src/fix_balance.cpp b/src/fix_balance.cpp index 2c76f13bf0..6b5c4b6ab1 100644 --- a/src/fix_balance.cpp +++ b/src/fix_balance.cpp @@ -13,20 +13,22 @@ ------------------------------------------------------------------------- */ #include "fix_balance.h" -#include -#include "balance.h" -#include "update.h" + #include "atom.h" +#include "balance.h" #include "comm.h" #include "domain.h" -#include "neighbor.h" -#include "irregular.h" +#include "error.h" +#include "fix_store.h" #include "force.h" +#include "irregular.h" #include "kspace.h" #include "modify.h" -#include "fix_store.h" +#include "neighbor.h" #include "rcb.h" -#include "error.h" +#include "update.h" + +#include using namespace LAMMPS_NS; using namespace FixConst; diff --git a/src/fix_move.cpp b/src/fix_move.cpp index 72fd9b75d2..f7bc4d3640 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -522,7 +522,7 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for wiggle: X = X0 + A sin(w*dt) + // for wiggle: X = X0 + A sin(w*dt) } else if (mstyle == WIGGLE) { double arg = omega_rotate * delta; @@ -578,19 +578,19 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for rotate by right-hand rule around omega: - // P = point = vector = point of rotation - // R = vector = axis of rotation - // w = omega of rotation (from period) - // X0 = xoriginal = initial coord of atom - // R0 = runit = unit vector for R - // D = X0 - P = vector from P to X0 - // C = (D dot R0) R0 = projection of atom coord onto R line - // A = D - C = vector from R line to X0 - // B = R0 cross A = vector perp to A in plane of rotation - // A,B define plane of circular rotation around R line - // X = P + C + A cos(w*dt) + B sin(w*dt) - // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) + // for rotate by right-hand rule around omega: + // P = point = vector = point of rotation + // R = vector = axis of rotation + // w = omega of rotation (from period) + // X0 = xoriginal = initial coord of atom + // R0 = runit = unit vector for R + // D = X0 - P = vector from P to X0 + // C = (D dot R0) R0 = projection of atom coord onto R line + // A = D - C = vector from R line to X0 + // B = R0 cross A = vector perp to A in plane of rotation + // A,B define plane of circular rotation around R line + // X = P + C + A cos(w*dt) + B sin(w*dt) + // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) } else if (mstyle == ROTATE) { double arg = omega_rotate * delta; @@ -707,10 +707,10 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for variable: compute x,v from variables - // NOTE: also allow for changes to extra attributes? - // omega, angmom, theta, quat - // only necessary if prescribed motion involves rotation + // for variable: compute x,v from variables + // NOTE: also allow for changes to extra attributes? + // omega, angmom, theta, quat + // only necessary if prescribed motion involves rotation } else if (mstyle == VARIABLE) { @@ -778,21 +778,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vxvarstr) { if (vxvarstyle == EQUAL) v[i][0] = vx; else v[i][0] = velocity[i][0]; - if (rmass) { - x[i][0] += dtv * v[i][0]; - } else { - x[i][0] += dtv * v[i][0]; - } + x[i][0] += dtv * v[i][0]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; - x[i][0] += dtv * v[i][0]; } else { dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; - x[i][0] += dtv * v[i][0]; } + x[i][0] += dtv * v[i][0]; } if (yvarstr && vyvarstr) { @@ -806,21 +801,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vyvarstr) { if (vyvarstyle == EQUAL) v[i][1] = vy; else v[i][1] = velocity[i][1]; - if (rmass) { - x[i][1] += dtv * v[i][1]; - } else { - x[i][1] += dtv * v[i][1]; - } + x[i][1] += dtv * v[i][1]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][1] += dtfm * f[i][1]; - x[i][1] += dtv * v[i][1]; } else { dtfm = dtf / mass[type[i]]; v[i][1] += dtfm * f[i][1]; - x[i][1] += dtv * v[i][1]; } + x[i][1] += dtv * v[i][1]; } if (zvarstr && vzvarstr) { @@ -834,21 +824,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vzvarstr) { if (vzvarstyle == EQUAL) v[i][2] = vz; else v[i][2] = velocity[i][2]; - if (rmass) { - x[i][2] += dtv * v[i][2]; - } else { - x[i][2] += dtv * v[i][2]; - } + x[i][2] += dtv * v[i][2]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][2] += dtfm * f[i][2]; - x[i][2] += dtv * v[i][2]; } else { dtfm = dtf / mass[type[i]]; v[i][2] += dtfm * f[i][2]; - x[i][2] += dtv * v[i][2]; } + x[i][2] += dtv * v[i][2]; } domain->remap_near(x[i],xold); diff --git a/src/force.cpp b/src/force.cpp index 6b7e9033ca..333a6e4715 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -133,14 +133,14 @@ void _noopt Force::create_factories() Force::~Force() { - delete [] pair_style; - delete [] bond_style; - delete [] angle_style; - delete [] dihedral_style; - delete [] improper_style; - delete [] kspace_style; + delete[] pair_style; + delete[] bond_style; + delete[] angle_style; + delete[] dihedral_style; + delete[] improper_style; + delete[] kspace_style; - delete [] pair_restart; + delete[] pair_restart; if (pair) delete pair; if (bond) delete bond; @@ -220,16 +220,16 @@ void Force::setup() void Force::create_pair(const std::string &style, int trysuffix) { - delete [] pair_style; + delete[] pair_style; if (pair) delete pair; - if (pair_restart) delete [] pair_restart; + if (pair_restart) delete[] pair_restart; pair_style = nullptr; pair = nullptr; pair_restart = nullptr; int sflag; pair = new_pair(style,trysuffix,sflag); - store_style(pair_style,style,sflag); + pair_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -345,12 +345,12 @@ char *Force::pair_match_ptr(Pair *ptr) void Force::create_bond(const std::string &style, int trysuffix) { - delete [] bond_style; + delete[] bond_style; if (bond) delete bond; int sflag; bond = new_bond(style,trysuffix,sflag); - store_style(bond_style,style,sflag); + bond_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -422,12 +422,12 @@ Bond *Force::bond_match(const std::string &style) void Force::create_angle(const std::string &style, int trysuffix) { - delete [] angle_style; + delete[] angle_style; if (angle) delete angle; int sflag; angle = new_angle(style,trysuffix,sflag); - store_style(angle_style,style,sflag); + angle_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -499,12 +499,12 @@ Angle *Force::angle_match(const std::string &style) void Force::create_dihedral(const std::string &style, int trysuffix) { - delete [] dihedral_style; + delete[] dihedral_style; if (dihedral) delete dihedral; int sflag; dihedral = new_dihedral(style,trysuffix,sflag); - store_style(dihedral_style,style,sflag); + dihedral_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -576,12 +576,12 @@ Dihedral *Force::dihedral_match(const std::string &style) void Force::create_improper(const std::string &style, int trysuffix) { - delete [] improper_style; + delete[] improper_style; if (improper) delete improper; int sflag; improper = new_improper(style,trysuffix,sflag); - store_style(improper_style,style,sflag); + improper_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -653,12 +653,12 @@ Improper *Force::improper_match(const std::string &style) void Force::create_kspace(const std::string &style, int trysuffix) { - delete [] kspace_style; + delete[] kspace_style; if (kspace) delete kspace; int sflag; kspace = new_kspace(style,trysuffix,sflag); - store_style(kspace_style,style,sflag); + kspace_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -729,14 +729,14 @@ KSpace *Force::kspace_match(const std::string &word, int exact) if sflag = 1/2/3, append suffix or suffix2 or suffixp to style ------------------------------------------------------------------------- */ -void Force::store_style(char *&str, const std::string &style, int sflag) +char *Force::store_style(const std::string &style, int sflag) { std::string estyle = style; if (sflag == 1) estyle += std::string("/") + lmp->suffix; else if (sflag == 2) estyle += std::string("/") + lmp->suffix2; else if (sflag == 3) estyle += std::string("/") + lmp->suffixp; - str = utils::strdup(estyle); + return utils::strdup(estyle); } /* ---------------------------------------------------------------------- diff --git a/src/force.h b/src/force.h index 09e13c7bd1..cc9ad06282 100644 --- a/src/force.h +++ b/src/force.h @@ -146,7 +146,7 @@ class Force : protected Pointers { KSpace *new_kspace(const std::string &, int, int &); KSpace *kspace_match(const std::string &, int); - void store_style(char *&, const std::string &, int); + char *store_style(const std::string &, int); void set_special(int, char **); double memory_usage(); diff --git a/src/improper_hybrid.cpp b/src/improper_hybrid.cpp index d3a9403a6b..0354f8e92e 100644 --- a/src/improper_hybrid.cpp +++ b/src/improper_hybrid.cpp @@ -241,7 +241,7 @@ void ImproperHybrid::settings(int narg, char **arg) error->all(FLERR, "Improper style hybrid cannot have none as an argument"); styles[nstyles] = force->new_improper(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/molecule.cpp b/src/molecule.cpp index 40ca218ecc..d6c839dfc4 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -57,8 +57,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : id = utils::strdup(arg[0]); if (!utils::is_id(id)) - error->all(FLERR,"Molecule template ID must have only " - "alphanumeric or underscore characters"); + error->all(FLERR,"Molecule template ID must have only alphanumeric or underscore characters"); // parse args until reach unknown arg (next file) @@ -130,8 +129,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : if (me == 0) { fp = fopen(arg[ifile],"r"); if (fp == nullptr) - error->one(FLERR,"Cannot open molecule file {}: {}", - arg[ifile], utils::getsyserror()); + error->one(FLERR,"Cannot open molecule file {}: {}",arg[ifile], utils::getsyserror()); } read(0); if (me == 0) fclose(fp); @@ -496,8 +494,7 @@ void Molecule::read(int flag) if (nmatch != nwant) error->one(FLERR,"Invalid header line format in molecule file"); } catch (TokenizerException &e) { - error->one(FLERR, "Invalid header in molecule file\n" - "{}", e.what()); + error->one(FLERR,"Invalid header in molecule file: {}", e.what()); } } @@ -618,11 +615,9 @@ void Molecule::read(int flag) // Error: Either a too long/short section or a typo in the keyword if (utils::strmatch(keyword,"^[A-Za-z ]+$")) - error->one(FLERR,"Unknown section '{}' in molecule " - "file\n",keyword); + error->one(FLERR,"Unknown section '{}' in molecule file\n",keyword); else error->one(FLERR,"Unexpected line in molecule file " - "while looking for the next " - "section:\n{}",line); + "while looking for the next section:\n{}",line); } keyword = parse_keyword(1,line); } @@ -648,8 +643,7 @@ void Molecule::read(int flag) if (bondflag && specialflag == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot auto-generate special bonds before " - "simulation box is defined"); + error->all(FLERR,"Cannot auto-generate special bonds before simulation box is defined"); if (flag) { special_generate(); @@ -691,8 +685,7 @@ void Molecule::coords(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Coords section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Coords section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -707,19 +700,16 @@ void Molecule::coords(char *line) x[iatom][2] *= sizescale; } } catch (TokenizerException &e) { - error->all(FLERR,"Invalid line in Coords section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Coords section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Coords " - "section of molecule file",i+1); + if (count[i] == 0) error->all(FLERR,"Atom {} missing in Coords section of molecule file",i+1); if (domain->dimension == 2) { for (int i = 0; i < natoms; i++) if (x[i][2] != 0.0) - error->all(FLERR,"Z coord in molecule file for atom {} " - "must be 0.0 for 2d-simulation.",i+1); + error->all(FLERR,"Z coord in molecule file for atom {} must be 0.0 for 2d-simulation",i+1); } } @@ -737,8 +727,7 @@ void Molecule::types(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Types section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Types section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -748,17 +737,14 @@ void Molecule::types(char *line) type[iatom] += toffset; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Types section of " - "molecule file: {}\n{}", e.what(),line); + error->all(FLERR,"Invalid line in Types section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Types " - "section of molecule file",i+1); + if (count[i] == 0) error->all(FLERR,"Atom {} missing in Types section of molecule file",i+1); if ((type[i] <= 0) || (domain->box_exist && (type[i] > atom->ntypes))) - error->all(FLERR,"Invalid atom type {} for atom {} " - "in molecule file",type[i],i+1); + error->all(FLERR,"Invalid atom type {} for atom {} in molecule file",type[i],i+1); ntypes = MAX(ntypes,type[i]); } @@ -777,8 +763,7 @@ void Molecule::molecules(char *line) readline(line); ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Molecules section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Molecules section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -788,19 +773,17 @@ void Molecule::molecules(char *line) // molecule[iatom] += moffset; // placeholder for possible molecule offset } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Molecules section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Molecules section of molecule file: {}\n{}",e.what(),line); } - for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Molecules " - "section of molecule file",i+1); - - for (int i = 0; i < natoms; i++) + for (int i = 0; i < natoms; i++) { + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Molecules section of molecule file",i+1); + } + for (int i = 0; i < natoms; i++) { if (molecule[i] < 0) - error->all(FLERR,"Invalid molecule ID {} for atom {} " - "in molecule file",molecule[i],i+1); - + error->all(FLERR,"Invalid molecule ID {} for atom {} in molecule file",molecule[i],i+1); + } for (int i = 0; i < natoms; i++) nmolecules = MAX(nmolecules,molecule[i]); } @@ -818,23 +801,21 @@ void Molecule::fragments(char *line) ValueTokenizer values(utils::trim_comment(line)); if ((int)values.count() > natoms+1) - error->all(FLERR,"Too many atoms per fragment in Fragments " - "section of molecule file"); + error->all(FLERR,"Too many atoms per fragment in Fragments section of molecule file"); fragmentnames[i] = values.next_string(); while (values.has_next()) { int iatom = values.next_int()-1; if (iatom < 0 || iatom >= natoms) - error->all(FLERR,"Invalid atom ID {} for fragment {} in " - "Fragments section of molecule file", - iatom+1, fragmentnames[i]); + error->all(FLERR,"Invalid atom ID {} for fragment {} in Fragments section of " + "molecule file", iatom+1, fragmentnames[i]); fragmentmask[i][iatom] = 1; } } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid atom ID in Fragments section of " - "molecule file: {}\n{}", e.what(),line); + error->all(FLERR,"Invalid atom ID in Fragments section of " + "molecule file: {}\n{}", e.what(),line); } } @@ -851,8 +832,7 @@ void Molecule::charges(char *line) ValueTokenizer values(utils::trim_comment(line)); if ((int)values.count() != 2) - error->all(FLERR,"Invalid line in Charges section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Charges section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -862,13 +842,13 @@ void Molecule::charges(char *line) q[iatom] = values.next_double(); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Charges section of " - "molecule file: {}.\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Charges section of molecule file: {}\n{}",e.what(),line); } - for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Charges " - "section of molecule file",i+1); + for (int i = 0; i < natoms; i++) { + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Charges section of molecule file",i+1); + } } /* ---------------------------------------------------------------------- @@ -885,8 +865,7 @@ void Molecule::diameters(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Diameters section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Diameters section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) error->all(FLERR,"Invalid atom index in Diameters section of molecule file"); @@ -897,16 +876,14 @@ void Molecule::diameters(char *line) maxradius = MAX(maxradius,radius[iatom]); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Diameters section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Diameters section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Diameters " - "section of molecule file",i+1); + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Diameters section of molecule file",i+1); if (radius[i] < 0.0) - error->all(FLERR,"Invalid atom diameter {} for atom {} " - "in molecule file", radius[i], i+1); + error->all(FLERR,"Invalid atom diameter {} for atom {} in molecule file", radius[i], i+1); } } @@ -923,8 +900,7 @@ void Molecule::masses(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Masses section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Masses section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -934,8 +910,7 @@ void Molecule::masses(char *line) rmass[iatom] *= sizescale*sizescale*sizescale; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Masses section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Masses section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { @@ -972,15 +947,13 @@ void Molecule::bonds(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Bonds section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Bonds section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); atom2 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Bonds section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Bonds section of molecule file: {}\n{}",e.what(),line); } itype += boffset; @@ -1042,16 +1015,14 @@ void Molecule::angles(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 5) - error->all(FLERR,"Invalid line in Angles section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Angles section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); atom2 = values.next_tagint(); atom3 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Angles section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Angles section of molecule file: {}\n{}",e.what(),line); } itype += aoffset; @@ -1128,8 +1099,7 @@ void Molecule::dihedrals(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 6) - error->all(FLERR,"Invalid line in Dihedrals section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Dihedrals section of molecule file: {}",line); values.next_int(); itype = values.next_int(); @@ -1138,8 +1108,7 @@ void Molecule::dihedrals(int flag, char *line) atom3 = values.next_tagint(); atom4 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Dihedrals section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Dihedrals section of molecule file: {}\n{}",e.what(),line); } itype += doffset; @@ -1150,8 +1119,7 @@ void Molecule::dihedrals(int flag, char *line) (atom4 <= 0) || (atom4 > natoms) || (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->all(FLERR, - "Invalid atom ID in dihedrals section of molecule file"); + error->all(FLERR,"Invalid atom ID in dihedrals section of molecule file"); if ((itype <= 0) || (domain->box_exist && (itype > atom->ndihedraltypes))) error->all(FLERR,"Invalid dihedral type in Dihedrals section of molecule file"); @@ -1230,8 +1198,7 @@ void Molecule::impropers(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 6) - error->all(FLERR,"Invalid line in Impropers section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Impropers section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); @@ -1239,8 +1206,7 @@ void Molecule::impropers(int flag, char *line) atom3 = values.next_tagint(); atom4 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Impropers section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Impropers section of molecule file: {}\n{}",e.what(),line); } itype += ioffset; @@ -1251,8 +1217,7 @@ void Molecule::impropers(int flag, char *line) (atom4 <= 0) || (atom4 > natoms) || (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->all(FLERR, - "Invalid atom ID in impropers section of molecule file"); + error->all(FLERR,"Invalid atom ID in impropers section of molecule file"); if ((itype <= 0) || (domain->box_exist && (itype > atom->nimpropertypes))) error->all(FLERR,"Invalid improper type in Impropers section of molecule file"); @@ -1325,15 +1290,14 @@ void Molecule::nspecial_read(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Special Bond Counts section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Special Bond Counts section of molecule file: {}",line); values.next_int(); c1 = values.next_tagint(); c2 = values.next_tagint(); c3 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Special Bond Counts section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Special Bond Counts section of " + "molecule file: {}\n{}",e.what(),line); } if (flag) { @@ -1358,8 +1322,7 @@ void Molecule::special_read(char *line) int nwords = values.count(); if (nwords != nspecial[i][2]+1) - error->all(FLERR,"Molecule file special list " - "does not match special count"); + error->all(FLERR,"Molecule file special list does not match special count"); values.next_int(); // ignore @@ -1367,13 +1330,12 @@ void Molecule::special_read(char *line) special[i][m-1] = values.next_tagint(); if (special[i][m-1] <= 0 || special[i][m-1] > natoms || special[i][m-1] == i+1) - error->all(FLERR,"Invalid atom index in Special Bonds " - "section of molecule file"); + error->all(FLERR,"Invalid atom index in Special Bonds section of molecule file"); } } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Special Bonds section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Special Bonds section of " + "molecule file: {}\n{}",e.what(),line); } } @@ -1502,8 +1464,7 @@ void Molecule::shakeflag_read(char *line) shake_flag[i] = values.next_int(); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid Shake Flags section in molecule file\n" - "{}", e.what()); + error->all(FLERR,"Invalid Shake Flags section in molecule file: {}", e.what()); } for (int i = 0; i < natoms; i++) @@ -1572,8 +1533,7 @@ void Molecule::shakeatom_read(char *line) } } catch (TokenizerException &e) { - error->all(FLERR,"Invalid shake atom in molecule file\n" - "{}", e.what()); + error->all(FLERR,"Invalid shake atom in molecule file: {}", e.what()); } for (int i = 0; i < natoms; i++) { @@ -1642,8 +1602,7 @@ void Molecule::shaketype_read(char *line) error->all(FLERR,"Invalid shake type data in molecule file"); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid shake type data in molecule file\n", - "{}", e.what()); + error->all(FLERR,"Invalid shake type data in molecule file: {}",e.what()); } for (int i = 0; i < natoms; i++) { @@ -1695,8 +1654,7 @@ void Molecule::body(int flag, int pflag, char *line) } else nword += ncount; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid body params in molecule file\n", - "{}", e.what()); + error->all(FLERR,"Invalid body params in molecule file: {}", e.what()); } } @@ -1735,8 +1693,7 @@ void Molecule::check_attributes(int flag) if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1; if (mismatch && me == 0) - error->warning(FLERR, - "Molecule attributes do not match system attributes"); + error->warning(FLERR,"Molecule attributes do not match system attributes"); // for all atom styles, check nbondtype,etc @@ -1770,8 +1727,7 @@ void Molecule::check_attributes(int flag) // warn if molecule topology defined but no special settings if (onemol->bondflag && !onemol->specialflag) - if (me == 0) error->warning(FLERR,"Molecule has bond topology " - "but no special bond settings"); + if (me == 0) error->warning(FLERR,"Molecule has bond topology but no special bond settings"); } } @@ -1878,47 +1834,31 @@ void Molecule::allocate() memory->create(special,natoms,maxspecial,"molecule:special"); if (bondflag) { - memory->create(bond_type,natoms,bond_per_atom, - "molecule:bond_type"); - memory->create(bond_atom,natoms,bond_per_atom, - "molecule:bond_atom"); + memory->create(bond_type,natoms,bond_per_atom,"molecule:bond_type"); + memory->create(bond_atom,natoms,bond_per_atom,"molecule:bond_atom"); } if (angleflag) { - memory->create(angle_type,natoms,angle_per_atom, - "molecule:angle_type"); - memory->create(angle_atom1,natoms,angle_per_atom, - "molecule:angle_atom1"); - memory->create(angle_atom2,natoms,angle_per_atom, - "molecule:angle_atom2"); - memory->create(angle_atom3,natoms,angle_per_atom, - "molecule:angle_atom3"); + memory->create(angle_type,natoms,angle_per_atom,"molecule:angle_type"); + memory->create(angle_atom1,natoms,angle_per_atom,"molecule:angle_atom1"); + memory->create(angle_atom2,natoms,angle_per_atom,"molecule:angle_atom2"); + memory->create(angle_atom3,natoms,angle_per_atom,"molecule:angle_atom3"); } if (dihedralflag) { - memory->create(dihedral_type,natoms,dihedral_per_atom, - "molecule:dihedral_type"); - memory->create(dihedral_atom1,natoms,dihedral_per_atom, - "molecule:dihedral_atom1"); - memory->create(dihedral_atom2,natoms,dihedral_per_atom, - "molecule:dihedral_atom2"); - memory->create(dihedral_atom3,natoms,dihedral_per_atom, - "molecule:dihedral_atom3"); - memory->create(dihedral_atom4,natoms,dihedral_per_atom, - "molecule:dihedral_atom4"); + memory->create(dihedral_type,natoms,dihedral_per_atom,"molecule:dihedral_type"); + memory->create(dihedral_atom1,natoms,dihedral_per_atom,"molecule:dihedral_atom1"); + memory->create(dihedral_atom2,natoms,dihedral_per_atom,"molecule:dihedral_atom2"); + memory->create(dihedral_atom3,natoms,dihedral_per_atom,"molecule:dihedral_atom3"); + memory->create(dihedral_atom4,natoms,dihedral_per_atom,"molecule:dihedral_atom4"); } if (improperflag) { - memory->create(improper_type,natoms,improper_per_atom, - "molecule:improper_type"); - memory->create(improper_atom1,natoms,improper_per_atom, - "molecule:improper_atom1"); - memory->create(improper_atom2,natoms,improper_per_atom, - "molecule:improper_atom2"); - memory->create(improper_atom3,natoms,improper_per_atom, - "molecule:improper_atom3"); - memory->create(improper_atom4,natoms,improper_per_atom, - "molecule:improper_atom4"); + memory->create(improper_type,natoms,improper_per_atom,"molecule:improper_type"); + memory->create(improper_atom1,natoms,improper_per_atom,"molecule:improper_atom1"); + memory->create(improper_atom2,natoms,improper_per_atom,"molecule:improper_atom2"); + memory->create(improper_atom3,natoms,improper_per_atom,"molecule:improper_atom3"); + memory->create(improper_atom4,natoms,improper_per_atom,"molecule:improper_atom4"); } if (shakeflag) { @@ -2058,7 +1998,7 @@ void Molecule::skip_lines(int n, char *line, const std::string §ion) readline(line); if (utils::strmatch(utils::trim(utils::trim_comment(line)),"^[A-Za-z ]+$")) error->one(FLERR,"Unexpected line in molecule file while " - "skipping {} section:\n{}",section,line); + "skipping {} section:\n{}",section,line); } } diff --git a/src/output.cpp b/src/output.cpp index a8fda4173d..1376d365d8 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -41,6 +41,8 @@ using namespace LAMMPS_NS; #define DELTA 1 #define EPSDT 1.0e-6 +enum {SETUP, WRITE, RESET_DT}; + /* ---------------------------------------------------------------------- initialize all output ------------------------------------------------------------------------- */ @@ -232,7 +234,7 @@ void Output::setup(int memflag) // only do this if dump written or dump has not been written yet if (writeflag || last_dump[idump] < 0) - calculate_next_dump(0,idump,ntimestep); + calculate_next_dump(SETUP,idump,ntimestep); // if dump not written now, use addstep_compute_all() // since don't know what computes the dump will invoke @@ -361,7 +363,7 @@ void Output::setup(int memflag) dump[idump]->write(); last_dump[idump] = ntimestep; - calculate_next_dump(1,idump,ntimestep); + calculate_next_dump(WRITE,idump,ntimestep); if (mode_dump[idump] == 0 && (dump[idump]->clearstep || var_dump[idump])) @@ -468,12 +470,12 @@ void Output::setup(int memflag) for timestep mode, set next_dump for simulation time mode, set next_time_dump and next_dump which flag depends on caller - 0 = from setup() at start of run - 1 = from write() during run each time a dump file is written - 2 = from reset_dt() called from fix dt/reset when it changes timestep size + SETUP = from setup() at start of run + WRITE = from write() during run each time a dump file is written + RESET_DT = from reset_dt() called from fix dt/reset when it changes timestep size ------------------------------------------------------------------------- */ - void Output::calculate_next_dump(int which, int idump, bigint ntimestep) +void Output::calculate_next_dump(int which, int idump, bigint ntimestep) { // dump mode is by timestep // just set next_dump @@ -482,13 +484,13 @@ void Output::setup(int memflag) if (every_dump[idump]) { - // which = 0: nextdump = next multiple of every_dump - // which = 1: increment nextdump by every_dump + // which = SETUP: nextdump = next multiple of every_dump + // which = WRITE: increment nextdump by every_dump - if (which == 0) + if (which == SETUP) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; - else if (which == 1) + else if (which == WRITE) next_dump[idump] += every_dump[idump]; } else { @@ -510,17 +512,28 @@ void Output::setup(int memflag) if (every_time_dump[idump] > 0.0) { - // which = 0: nexttime = next multiple of every_time_dump - // which = 1: increment nexttime by every_time_dump - // which = 2: no change to previous nexttime (only timestep has changed) + // which = SETUP: nexttime = next multiple of every_time_dump + // which = WRITE: increment nexttime by every_time_dump + // which = RESET_DT: no change to previous nexttime (only timestep has changed) - if (which == 0) + switch (which) { + case SETUP: nexttime = static_cast (tcurrent/every_time_dump[idump]) * every_time_dump[idump] + every_time_dump[idump]; - else if (which == 1) + break; + + case WRITE: nexttime = next_time_dump[idump] + every_time_dump[idump]; - else if (which == 2) + break; + + case RESET_DT: nexttime = next_time_dump[idump]; + break; + + default: + nexttime = 0; + error->all(FLERR,"Unexpected argument to calculate_next_dump"); + } nextdump = ntimestep + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / @@ -541,10 +554,10 @@ void Output::setup(int memflag) } else { - // do not re-evaulate variable for which = 2, leave nexttime as-is + // do not re-evaulate variable for which = RESET_DT, leave nexttime as-is // unless next_time_dump < 0.0, which means variable never yet evaluated - if (which < 2 || next_time_dump[idump] < 0.0) { + if (which < RESET_DT || next_time_dump[idump] < 0.0) { nexttime = input->variable->compute_equal(ivar_dump[idump]); } else nexttime = next_time_dump[idump]; @@ -619,9 +632,8 @@ int Output::check_time_dumps(bigint ntimestep) { next_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) - if (last_dump[idump] >= 0) - error->all(FLERR, - "Cannot reset timestep with active dump - must undump first"); + if ((last_dump[idump] >= 0) && !update->whichflag && !dump[idump]->multifile) + error->all(FLERR, "Cannot reset timestep with active dump - must undump first"); if (restart_flag_single) { if (restart_every_single) { @@ -704,7 +716,7 @@ void Output::reset_dt() // since timestep change affects next step if (next_dump[idump] != ntimestep) - calculate_next_dump(2,idump,update->ntimestep); + calculate_next_dump(RESET_DT,idump,update->ntimestep); if (dump[idump]->clearstep || var_dump[idump]) next_time_dump_any = MIN(next_time_dump_any,next_dump[idump]); diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index e962e02c9e..6287fb97db 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -50,20 +50,20 @@ PairHybrid::~PairHybrid() if (nstyles > 0) { for (int m = 0; m < nstyles; m++) { delete styles[m]; - delete [] keywords[m]; - if (special_lj[m]) delete [] special_lj[m]; - if (special_coul[m]) delete [] special_coul[m]; + delete[] keywords[m]; + if (special_lj[m]) delete[] special_lj[m]; + if (special_coul[m]) delete[] special_coul[m]; } } - delete [] styles; - delete [] keywords; - delete [] multiple; + delete[] styles; + delete[] keywords; + delete[] multiple; - delete [] special_lj; - delete [] special_coul; - delete [] compute_tally; + delete[] special_lj; + delete[] special_coul; + delete[] compute_tally; - delete [] svector; + delete[] svector; if (allocated) { memory->destroy(setflag); @@ -187,7 +187,7 @@ void PairHybrid::compute(int eflag, int vflag) } - delete [] saved_special; + delete[] saved_special; if (vflag_fdotr) virial_fdotr_compute(); } @@ -274,9 +274,9 @@ void PairHybrid::settings(int narg, char **arg) if (nstyles > 0) { for (int m = 0; m < nstyles; m++) { delete styles[m]; - delete [] keywords[m]; - if (special_lj[m]) delete [] special_lj[m]; - if (special_coul[m]) delete [] special_coul[m]; + delete[] keywords[m]; + if (special_lj[m]) delete[] special_lj[m]; + if (special_coul[m]) delete[] special_coul[m]; } delete[] styles; delete[] keywords; @@ -297,8 +297,8 @@ void PairHybrid::settings(int narg, char **arg) // allocate list of sub-styles as big as possibly needed if no extra args - styles = new Pair*[narg]; - keywords = new char*[narg]; + styles = new Pair *[narg]; + keywords = new char *[narg]; multiple = new int[narg]; special_lj = new double*[narg]; @@ -322,7 +322,7 @@ void PairHybrid::settings(int narg, char **arg) error->all(FLERR,"Pair style hybrid cannot have none as an argument"); styles[nstyles] = force->new_pair(arg[iarg],1,dummy); - force->store_style(keywords[nstyles],arg[iarg],0); + keywords[nstyles] = force->store_style(arg[iarg],0); special_lj[nstyles] = special_coul[nstyles] = nullptr; compute_tally[nstyles] = 1; @@ -454,7 +454,7 @@ void PairHybrid::init_svector() single_extra = MAX(single_extra,styles[m]->single_extra); if (single_extra) { - delete [] svector; + delete[] svector; svector = new double[single_extra]; } } @@ -667,7 +667,7 @@ void PairHybrid::init_style() neighbor->requests[i]->iskip = iskip; neighbor->requests[i]->ijskip = ijskip; } else { - delete [] iskip; + delete[] iskip; memory->destroy(ijskip); } } diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 5bf593d147..68a6199e19 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -328,7 +328,7 @@ void PairHybridScaled::settings(int narg, char **arg) error->all(FLERR, "Pair style hybrid/scaled cannot have none as an argument"); styles[nstyles] = force->new_pair(arg[iarg], 1, dummy); - force->store_style(keywords[nstyles], arg[iarg], 0); + keywords[nstyles] = force->store_style(arg[iarg], 0); special_lj[nstyles] = special_coul[nstyles] = nullptr; compute_tally[nstyles] = 1; @@ -524,7 +524,7 @@ void PairHybridScaled::write_restart(FILE *fp) int n = scalevars.size(); fwrite(&n, sizeof(int), 1, fp); - for (auto var : scalevars) { + for (auto &var : scalevars) { n = var.size() + 1; fwrite(&n, sizeof(int), 1, fp); fwrite(var.c_str(), sizeof(char), n, fp); diff --git a/src/read_data.cpp b/src/read_data.cpp index 89b3ecadec..dda418b8a2 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -903,23 +903,28 @@ void ReadData::command(int narg, char **arg) // restore old styles, when reading with nocoeff flag given if (coeffflag == 0) { - if (force->pair) delete force->pair; + delete force->pair; + delete[] force->pair_style; force->pair = saved_pair; force->pair_style = saved_pair_style; - if (force->bond) delete force->bond; + delete force->bond; + delete[] force->bond_style; force->bond = saved_bond; force->bond_style = saved_bond_style; - if (force->angle) delete force->angle; + delete force->angle; + delete[] force->angle_style; force->angle = saved_angle; force->angle_style = saved_angle_style; - if (force->dihedral) delete force->dihedral; + delete force->dihedral; + delete[] force->dihedral_style; force->dihedral = saved_dihedral; force->dihedral_style = saved_dihedral_style; - if (force->improper) delete force->improper; + delete force->improper; + delete[] force->improper_style; force->improper = saved_improper; force->improper_style = saved_improper_style; diff --git a/src/reader_native.cpp b/src/reader_native.cpp index ba7a576a50..32b2279a60 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -130,6 +130,7 @@ void ReaderNative::skip() // read chunk and skip them read_buf(&nchunk, sizeof(int), 1); + if (nchunk < 0) error->one(FLERR,"Dump file is invalid or corrupted"); int n; for (int i = 0; i < nchunk; i++) { @@ -141,8 +142,7 @@ void ReaderNative::skip() read_lines(2); bigint natoms; int rv = sscanf(line,BIGINT_FORMAT,&natoms); - if (rv != 1) - error->one(FLERR,"Dump file is incorrectly formatted"); + if (rv != 1) error->one(FLERR,"Dump file is incorrectly formatted"); read_lines(5); @@ -163,20 +163,17 @@ void ReaderNative::skip_reading_magic_str() if (is_known_magic_str() && revision > 0x0001) { int len; read_buf(&len, sizeof(int), 1); + if (len < 0) error->one(FLERR,"Dump file is invalid or corrupted"); - if (len > 0) { - // has units - skip_buf(sizeof(char)*len); - } + // has units + if (len > 0) skip_buf(sizeof(char)*len); char flag = 0; read_buf(&flag, sizeof(char), 1); - - if (flag) { - skip_buf(sizeof(double)); - } + if (flag) skip_buf(sizeof(double)); read_buf(&len, sizeof(int), 1); + if (len < 0) error->one(FLERR,"Dump file is invalid or corrupted"); skip_buf(sizeof(char)*len); } } diff --git a/src/tabular_function.cpp b/src/tabular_function.cpp index a3a904d644..07df49729a 100644 --- a/src/tabular_function.cpp +++ b/src/tabular_function.cpp @@ -1,6 +1,6 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/ Sandia National Laboratories + https://www.lammps.org/, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract diff --git a/src/version.h b/src/version.h index 6bb32291be..c33ca6ee15 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "14 Dec 2021" +#define LAMMPS_VERSION "7 Jan 2022" diff --git a/tools/singularity/fedora34_mingw.def b/tools/singularity/fedora35_mingw.def similarity index 99% rename from tools/singularity/fedora34_mingw.def rename to tools/singularity/fedora35_mingw.def index 8a1f108aec..c2108a378c 100644 --- a/tools/singularity/fedora34_mingw.def +++ b/tools/singularity/fedora35_mingw.def @@ -1,5 +1,5 @@ BootStrap: docker -From: fedora:34 +From: fedora:35 %post dnf -y update @@ -94,7 +94,7 @@ EOF CUSTOM_PROMPT_ENV=/.singularity.d/env/99-zz_custom_prompt.sh cat >$CUSTOM_PROMPT_ENV <