Merge pull request lammps#4719 from akohlmey/collected-small-changes

Collected small changes and fixes

Author: akohlmey

doc/src/Build_extras.rst

@@ -243,7 +243,7 @@ necessary for ``hipcc`` and the linker to work correctly.
 When compiling for HIP ROCm, GPU sorting with ``-D
 HIP_USE_DEVICE_SORT=on`` requires installing the ``hipcub`` library
 (https://github.com/ROCmSoftwarePlatform/hipCUB).  The HIP CUDA-backend
-additionally requires cub (https://nvlabs.github.io/cub).  Setting
+additionally requires cub (https://nvidia.github.io/cccl/cub/).  Setting
 ``-DDOWNLOAD_CUB=yes`` will download and compile CUB.
 
 The GPU library has some multi-thread support using OpenMP.  If LAMMPS

doc/src/Errors_details.rst

@@ -1068,3 +1068,22 @@ also increase the value for the "page" parameter to maintain the ratio
 between "one" and "page" to reduce waste of memory.  For some more
 details, please check out the documentation for the :doc:`neigh_modify
 command <neigh_modify>`.
+
+.. _err0037:
+
+Variable ...: Compute/Fix ... does not compute requested property
+-----------------------------------------------------------------
+
+Compute and fix styles can compute different kinds of properties: for
+example, global scalars, vectors, or arrays, or per-atom vectors or
+arrays.  In equal-style or similar variable, only scalar properties can
+be used, so to access a particular element in a vector one has to use
+square brackets with a suitable index to select it.  However, not all
+fixes and computes provide all types of properties.  So this error
+message will be shown if there is a mismatch, of if there are not
+enough or too many square brackets.  To differentiate between
+accessing an element of a global array or a per-atom array element of
+a specific atom, one has to use a reference with a lower case 'c'
+(e.g. 'c_name') for the former and upper case 'C' (e.g. 'C_name') for
+the latter. The same applies to fix styles.  The full details are
+in the documentation for the :doc:`variable command <variable>`.

doc/src/Howto_mdi.rst

@@ -138,7 +138,7 @@ the ``examples/QUANTUM`` sub-directories for more details:
 
 There are also at least two quantum codes which have direct MDI
 support, `Quantum ESPRESSO (QE) <https://www.quantum-espresso.org/>`_
-and `INQ <https://qsg.llnl.gov/node/101.html>`_.  There are also
+and `INQ <https://gitlab.com/npneq/inq>`_.  There are also
 several QM codes which have indirect support through QCEngine or i-PI.
 The former means they require a wrapper program (QCEngine) with MDI
 support which writes/read files to pass data to the quantum code

doc/src/fix_external.rst

@@ -71,14 +71,8 @@ which is typically a 32-bit integer unless LAMMPS is compiled with
 <size>` section of the manual.  Finally, *fexternal* are the forces
 returned by the driver program.
 
-The fix has a set_callback() method which the external driver can call
-to pass a pointer to its foo() function.  See the
-couple/lammps_quest/lmpqst.cpp file in the LAMMPS distribution for an
-example of how this is done.  This sample application performs
-classical MD using quantum forces computed by a density functional
-code `Quest <quest_>`_.
-
-.. _quest: https://www.sandia.gov/quest/
+The best way to set up the callback function is to use the C-language
+library interface function :cpp:func:`lammps_set_fix_external_callback`.
 
 ----------
 

doc/src/fix_phonon.rst

@@ -134,7 +134,7 @@ for other systems, *nasr* = 10 is typically sufficient.
 The *map_file* contains the mapping information between the lattice
 indices and the atom IDs, which tells the code which atom sits at
 which lattice point; the lattice indices start from 0. An auxiliary
-code, `latgen <https://code.google.com/p/latgen>`_, can be employed to
+code, `latgen <https://code.google.com/archive/p/latgen>`_, can be employed to
 generate the compatible map file for various crystals.
 
 In case one simulates a non-periodic system, where the whole simulation

src/EXTRA-DUMP/dump_extxyz.cpp

@@ -33,7 +33,7 @@ static constexpr int DELTA = 1048576;
 /* ---------------------------------------------------------------------- */
 
 DumpExtXYZ::DumpExtXYZ(LAMMPS *lmp, int narg, char **arg) :
-    DumpXYZ(lmp, narg, arg), properties_string(nullptr)
+    DumpXYZ(lmp, narg, arg), properties_string(nullptr), thermo_string(nullptr)
 {
   // style specific customizable settings
   with_vel = 1;
@@ -65,6 +65,7 @@ DumpExtXYZ::DumpExtXYZ(LAMMPS *lmp, int narg, char **arg) :
 DumpExtXYZ::~DumpExtXYZ()
 {
   delete[] properties_string;
+  delete[] thermo_string;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -135,37 +136,56 @@ int DumpExtXYZ::modify_param(int narg, char **arg)
 
 /* ---------------------------------------------------------------------- */
 
+int DumpExtXYZ::count()
+{
+  // Dump::count() is called for all processes before writing the header.
+  // We must update the thermodynamic data in the header here, since we are
+  // calling computes that may perform communication. We only include data
+  // from computes that were already invoked to avoid XXX not tallied errors.
+
+  delete[] thermo_string;
+  std::string thermo_buf;
+
+  if (time_flag) thermo_buf += fmt::format(" Time={:.6f}", compute_time());
+  thermo_buf += fmt::format(" pbc=\"{} {} {}\"", domain->xperiodic ? "T" : "F",
+                            domain->yperiodic ? "T" : "F", domain->zperiodic ? "T" : "F");
+  thermo_buf +=
+      fmt::format(" Lattice=\"{:g} {:g} {:g} {:g} {:g} {:g} {:g} {:g} {:g}\"", domain->xprd, 0., 0.,
+                  domain->xy, domain->yprd, 0., domain->xz, domain->yz, domain->zprd);
+
+  if (output && output->thermo) {
+    auto *pe = output->thermo->pe;
+    if (pe && (pe->invoked_scalar == update->ntimestep))
+      thermo_buf += fmt::format(" Potential_energy={}", pe->compute_scalar());
+
+    auto *temp = output->thermo->temperature;
+    if (temp && (temp->invoked_scalar == update->ntimestep))
+      thermo_buf += fmt::format(" Temperature={}", temp->compute_scalar());
+
+    auto *press = output->thermo->pressure;
+    if (press && (press->invoked_vector == update->ntimestep)) {
+      press->compute_vector();
+      thermo_buf +=
+          fmt::format(" Stress=\"{} {} {} {} {} {} {} {} {}\"", press->vector[0], press->vector[3],
+                      press->vector[4], press->vector[3], press->vector[1], press->vector[5],
+                      press->vector[4], press->vector[5], press->vector[2]);
+    }
+  }
+  thermo_string = utils::strdup(thermo_buf);
+
+  // perform real count() function in parent
+  return DumpXYZ::count();
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpExtXYZ::write_header(bigint n)
 {
   if (me == 0) {
     if (!fp)
       error->one(FLERR, Error::NOLASTLINE, "Must not use 'run pre no' after creating a new dump");
 
-    std::string header = fmt::format("{}\nTimestep={}", n, update->ntimestep);
-    if (time_flag) header += fmt::format(" Time={:.6f}", compute_time());
-    header += fmt::format(" pbc=\"{} {} {}\"", domain->xperiodic ? "T" : "F",
-                          domain->yperiodic ? "T" : "F", domain->zperiodic ? "T" : "F");
-    header +=
-        fmt::format(" Lattice=\"{:g} {:g} {:g} {:g} {:g} {:g} {:g} {:g} {:g}\"", domain->xprd, 0.,
-                    0., domain->xy, domain->yprd, 0., domain->xz, domain->yz, domain->zprd);
-
-    if (output && output->thermo) {
-      auto *pe = output->thermo->pe;
-      if (pe) header += fmt::format(" Potential_energy={}", pe->compute_scalar());
-
-      auto *temp = output->thermo->temperature;
-      if (temp) header += fmt::format(" Temperature={}", temp->compute_scalar());
-
-      auto *press = output->thermo->pressure;
-      if (press) {
-        press->compute_vector();
-        header +=
-            fmt::format(" Stress=\"{} {} {} {} {} {} {} {} {}\"", press->vector[0],
-                        press->vector[3], press->vector[4], press->vector[3], press->vector[1],
-                        press->vector[5], press->vector[4], press->vector[5], press->vector[2]);
-      }
-    }
-
+    auto header = fmt::format("{}\nTimestep={}{}", n, update->ntimestep, thermo_string);
     header += fmt::format(" Properties={}", properties_string);
     utils::print(fp, header + "\n");
   }

src/EXTRA-DUMP/dump_extxyz.h

@@ -36,9 +36,11 @@ class DumpExtXYZ : public DumpXYZ {
   int with_temp;
   int with_press;
   char *properties_string;
+  char *thermo_string;
 
   void update_properties();
   void init_style() override;
+  int count() override;
   void write_header(bigint) override;
   void pack(tagint *) override;
   int convert_string(int, double *) override;

src/KOKKOS/fix_spring_self_kokkos.cpp

@@ -80,10 +80,11 @@ void FixSpringSelfKokkos<DeviceType>::init()
   FixSpringSelf::init();
 
   if (kstyle != CONSTANT)
-    error->all(FLERR, "Fix spring/self/kk does not support variable spring constants (yet)");
+    error->all(FLERR, Error::NOLASTLINE,
+               "Fix spring/self/kk does not support variable spring constants (yet)");
 
   if (utils::strmatch(update->integrate_style,"^respa"))
-    error->all(FLERR,"Cannot (yet) use respa with Kokkos");
+    error->all(FLERR, Error::NOLASTLINE, "Cannot (yet) use respa with Kokkos");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -107,47 +108,45 @@ void FixSpringSelfKokkos<DeviceType>::post_force(int /*vflag*/)
   copymode = 1;
 
   {
-  // local variables for lambda capture
-  auto prd = Few<double,3>(domain->prd);
-  auto h = Few<double,6>(domain->h);
-  auto triclinic = domain->triclinic;
-  auto l_k = k;
-  auto l_xoriginal = d_xoriginal;
-
-  auto l_x = x;
-  auto l_f = f;
-  auto l_mask = mask;
-  auto l_image = image;
-  auto l_groupbit = groupbit;
-  auto l_xflag = xflag;
-  auto l_yflag = yflag;
-  auto l_zflag = zflag;
-
-  Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType>(0,nlocal), LAMMPS_LAMBDA(const int& i, double& espring_kk) {
-    if (l_mask[i] & l_groupbit) {
-      Few<double,3> x_i;
-      x_i[0] = l_x(i,0);
-      x_i[1] = l_x(i,1);
-      x_i[2] = l_x(i,2);
-      auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,l_image(i));
-      auto dx = unwrap[0] - l_xoriginal(i, 0);
-      auto dy = unwrap[1] - l_xoriginal(i, 1);
-      auto dz = unwrap[2] - l_xoriginal(i, 2);
-      if (!l_xflag) dx = 0.0;
-      if (!l_yflag) dy = 0.0;
-      if (!l_zflag) dz = 0.0;
-      l_f(i,0) -= l_k*dx;
-      l_f(i,1) -= l_k*dy;
-      l_f(i,2) -= l_k*dz;
-      espring_kk += l_k * (dx*dx + dy*dy + dz*dz);
-    }
-  },espring_kk);
+    // local variables for lambda capture
+    auto prd = Few<double,3>(domain->prd);
+    auto h = Few<double,6>(domain->h);
+    auto triclinic = domain->triclinic;
+    auto l_k = k;
+    auto l_xoriginal = d_xoriginal;
+
+    auto l_x = x;
+    auto l_f = f;
+    auto l_mask = mask;
+    auto l_image = image;
+    auto l_groupbit = groupbit;
+    auto l_xflag = xflag;
+    auto l_yflag = yflag;
+    auto l_zflag = zflag;
+
+    Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType>(0,nlocal), LAMMPS_LAMBDA(const int& i, double& espring_kk) {
+        if (l_mask[i] & l_groupbit) {
+          Few<double,3> x_i;
+          x_i[0] = l_x(i,0);
+          x_i[1] = l_x(i,1);
+          x_i[2] = l_x(i,2);
+          auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,l_image(i));
+          auto dx = unwrap[0] - l_xoriginal(i, 0);
+          auto dy = unwrap[1] - l_xoriginal(i, 1);
+          auto dz = unwrap[2] - l_xoriginal(i, 2);
+          if (!l_xflag) dx = 0.0;
+          if (!l_yflag) dy = 0.0;
+          if (!l_zflag) dz = 0.0;
+          l_f(i,0) -= l_k*dx;
+          l_f(i,1) -= l_k*dy;
+          l_f(i,2) -= l_k*dz;
+          espring_kk += l_k * (dx*dx + dy*dy + dz*dz);
+        }
+      },espring_kk);
   }
 
   copymode = 0;
-
   atomKK->modified(execution_space, F_MASK);
-
   espring = 0.5*espring_kk;
 }
 

src/KSPACE/fft3d.cpp

@@ -251,7 +251,7 @@ struct fft_plan_3d *fft_3d_create_plan(
        int scaled, int permute, int *nbuf, int usecollective)
 {
   struct fft_plan_3d *plan;
-  int me,nprocs,nthreads;
+  int me,nprocs;
   int flag,remapflag;
   int first_ilo,first_ihi,first_jlo,first_jhi,first_klo,first_khi;
   int second_ilo,second_ihi,second_jlo,second_jhi,second_klo,second_khi;
@@ -264,12 +264,14 @@ struct fft_plan_3d *fft_3d_create_plan(
   MPI_Comm_rank(comm,&me);
   MPI_Comm_size(comm,&nprocs);
 
+#if defined(FFT_MKL_THREADS) || defined(FFT_FFTW_THREADS)
 #if defined(_OPENMP)
   // query OpenMP info.
   // should have been initialized systemwide in Comm class constructor
-  nthreads = omp_get_max_threads();
+  int nthreads = omp_get_max_threads();
 #else
-  nthreads = 1;
+  int nthreads = 1;
+#endif
 #endif
 
   // compute division of procs in 2 dimensions not on-processor

src/KSPACE/kissfft.h

@@ -158,7 +158,7 @@ static void kiss_fft(kiss_fft_cfg, const FFT_DATA *, FFT_DATA *);
     (x)->im = KISS_FFT_SIN(phase); \
   } while (0)
 
-static void kf_bfly2(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st, int m)
+static void kf_bfly2(FFT_DATA *Fout, const size_t fstride, kiss_fft_cfg st, int m)
 {
   FFT_DATA *Fout2;
   FFT_DATA *tw1 = st->twiddles;
@@ -178,7 +178,7 @@ static void kf_bfly2(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st
   } while (--m);
 }
 
-static void kf_bfly4(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st, const size_t m)
+static void kf_bfly4(FFT_DATA *Fout, const size_t fstride, kiss_fft_cfg st, const size_t m)
 {
   FFT_DATA *tw1, *tw2, *tw3;
   FFT_DATA scratch[6];
@@ -223,7 +223,7 @@ static void kf_bfly4(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st
   } while (--k);
 }
 
-static void kf_bfly3(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st, size_t m)
+static void kf_bfly3(FFT_DATA *Fout, const size_t fstride, kiss_fft_cfg st, size_t m)
 {
   size_t k = m;
   const size_t m2 = 2 * m;
@@ -264,7 +264,7 @@ static void kf_bfly3(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st
   } while (--k);
 }
 
-static void kf_bfly5(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st, int m)
+static void kf_bfly5(FFT_DATA *Fout, const size_t fstride, kiss_fft_cfg st, int m)
 {
   FFT_DATA *Fout0, *Fout1, *Fout2, *Fout3, *Fout4;
   int u;
@@ -329,15 +329,15 @@ static void kf_bfly5(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st
 }
 
 /* perform the butterfly for one stage of a mixed radix FFT */
-static void kf_bfly_generic(FFT_DATA *Fout, const size_t fstride, const kiss_fft_cfg st, int m,
+static void kf_bfly_generic(FFT_DATA *Fout, const size_t fstride, kiss_fft_cfg st, int m,
                             int p)
 {
   int u, k, q1, q;
   FFT_DATA *twiddles = st->twiddles;
   FFT_DATA t;
   int Norig = st->nfft;
 
-  FFT_DATA *scratch = (FFT_DATA *) KISS_FFT_TMP_ALLOC(sizeof(FFT_DATA) * p);
+  auto *scratch = (FFT_DATA *) KISS_FFT_TMP_ALLOC(sizeof(FFT_DATA) * p);
   for (u = 0; u < m; ++u) {
     k = u;
     for (q1 = 0; q1 < p; ++q1) {
@@ -363,7 +363,7 @@ static void kf_bfly_generic(FFT_DATA *Fout, const size_t fstride, const kiss_fft
 }
 
 static void kf_work(FFT_DATA *Fout, const FFT_DATA *f, const size_t fstride, int in_stride,
-                    int *factors, const kiss_fft_cfg st)
+                    int *factors, kiss_fft_cfg st)
 {
   FFT_DATA *Fout_beg = Fout;
   const int p = *factors++; /* the radix  */
@@ -482,7 +482,7 @@ static void kiss_fft_stride(kiss_fft_cfg st, const FFT_DATA *fin, FFT_DATA *fout
   if (fin == fout) {
     // NOTE: this is not really an in-place FFT algorithm.
     // It just performs an out-of-place FFT into a temp buffer
-    FFT_DATA *tmpbuf = (FFT_DATA *) KISS_FFT_TMP_ALLOC(sizeof(FFT_DATA) * st->nfft);
+    auto *tmpbuf = (FFT_DATA *) KISS_FFT_TMP_ALLOC(sizeof(FFT_DATA) * st->nfft);
     kf_work(tmpbuf, fin, 1, in_stride, st->factors, st);
     memcpy(fout, tmpbuf, sizeof(FFT_DATA) * st->nfft);
     KISS_FFT_TMP_FREE(tmpbuf);