Merge pull request lammps#4037 from akohlmey/bigint_dof_computation

Make computation of PPPM grid counts and DOFs removed by fixes compatible with large systems

cmake/CMakeLists.txt

@@ -662,9 +662,14 @@ foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT
   endif()
 endforeach()
 
-if(PKG_PLUGIN)
-  target_compile_definitions(lammps PRIVATE -DLMP_PLUGIN)
-endif()
+######################################################################
+# packages with defines to disable package specific code
+######################################################################
+foreach(PKG_WITH_DEF  BPM PLUGIN)
+  if(PKG_${PKG_WITH_DEF})
+    target_compile_definitions(lammps PRIVATE -DLMP_${PKG_WITH_DEF})
+  endif()
+endforeach()
 
 # link with -ldl or equivalent for plugin loading; except on Windows
 if(NOT ${CMAKE_SYSTEM_NAME} STREQUAL "Windows")

src/BPM/Install.sh

@@ -0,0 +1,49 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+mode=$1
+
+# enforce using portable C locale
+LC_ALL=C
+export LC_ALL
+
+# arg1 = file, arg2 = file it depends on
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+# all package files with no dependencies
+
+for file in *.cpp *.h; do
+  test -f ${file} && action $file
+done
+
+# edit 2 Makefile.package files to include/exclude package info
+
+if (test $1 = 1) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's|^PKG_INC =[ \t]*|&-DLMP_BPM |' ../Makefile.package
+  fi
+
+elif (test $1 = 0) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*LMP_BPM[^ \t]* //' ../Makefile.package
+  fi
+
+fi

src/EFF/fix_langevin_eff.cpp

@@ -137,7 +137,7 @@ void FixLangevinEff::post_force_no_tally()
   dof = domain->dimension * particles;
   fix_dof = 0;
   for (int i = 0; i < modify->nfix; i++)
-    fix_dof += modify->fix[i]->dof(igroup);
+    fix_dof += (int)modify->fix[i]->dof(igroup);
 
   // extra_dof = domain->dimension
   dof -= domain->dimension + fix_dof;
@@ -306,7 +306,7 @@ void FixLangevinEff::post_force_tally()
   dof = domain->dimension * particles;
   fix_dof = 0;
   for (int i = 0; i < modify->nfix; i++)
-    fix_dof += modify->fix[i]->dof(igroup);
+    fix_dof += (int)modify->fix[i]->dof(igroup);
 
   // extra_dof = domain->dimension
   dof -= domain->dimension + fix_dof;

src/ELECTRODE/pppm_electrode.cpp

@@ -633,7 +633,9 @@ void PPPMElectrode::project_psi(double *vec, int sensor_grpbit)
   // project u_brick with weight matrix
   double **x = atom->x;
   int *mask = atom->mask;
-  double const scaleinv = 1.0 / (nx_pppm * ny_pppm * nz_pppm);
+  const bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  const double scaleinv = 1.0 / ngridtotal;
+
   for (int i = 0; i < atom->nlocal; i++) {
     if (!(mask[i] & sensor_grpbit)) continue;
     double v = 0.;
@@ -1362,7 +1364,7 @@ double PPPMElectrode::compute_qopt()
   // each proc calculates contributions from every Pth grid point
 
   bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
-  int nxy_pppm = nx_pppm * ny_pppm;
+  bigint nxy_pppm = (bigint) nx_pppm * ny_pppm;
 
   double qopt = 0.0;
 

src/GPU/pppm_gpu.cpp

@@ -405,7 +405,8 @@ void PPPMGPU::poisson_ik()
 
   // if requested, compute energy and virial contribution
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0 / ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {

src/INTEL/pppm_electrode_intel.cpp

@@ -420,7 +420,9 @@ void PPPMElectrodeIntel::project_psi(IntelBuffers<flt_t, acc_t> *buffers, double
 #endif
   {
     int *mask = atom->mask;
-    const flt_t scaleinv = 1.0 / (nx_pppm * ny_pppm * nz_pppm);
+
+    const bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+    const flt_t scaleinv = 1.0 / ngridtotal;
 
     const flt_t lo0 = boxlo[0];
     const flt_t lo1 = boxlo[1];

src/KOKKOS/fix_shake_kokkos.cpp

@@ -525,7 +525,7 @@ void FixShakeKokkos<DeviceType>::operator()(TagFixShakePostForce<NEIGHFLAG,EVFLA
 ------------------------------------------------------------------------- */
 
 template<class DeviceType>
-int FixShakeKokkos<DeviceType>::dof(int igroup)
+bigint FixShakeKokkos<DeviceType>::dof(int igroup)
 {
   d_mask = atomKK->k_mask.view<DeviceType>();
   d_tag = atomKK->k_tag.view<DeviceType>();
@@ -538,7 +538,7 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
   // count dof in a cluster if and only if
   // the central atom is in group and atom i is the central atom
 
-  int n = 0;
+  bigint n = 0;
   {
     // local variables for lambda capture
 
@@ -549,7 +549,7 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
     auto groupbit = group->bitmask[igroup];
 
     Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType>(0,nlocal),
-     LAMMPS_LAMBDA(const int& i, int& n) {
+     LAMMPS_LAMBDA(const int& i, bigint& n) {
       if (!(mask[i] & groupbit)) return;
       if (d_shake_flag[i] == 0) return;
       if (d_shake_atom(i,0) != tag[i]) return;
@@ -560,8 +560,8 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
     },n);
   }
 
-  int nall;
-  MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
+  bigint nall;
+  MPI_Allreduce(&n,&nall,1,MPI_LMP_BIGINT,MPI_SUM,world);
   return nall;
 }
 

src/KOKKOS/fix_shake_kokkos.h

@@ -44,8 +44,6 @@ struct TagFixShakeUnpackExchange{};
 template<class DeviceType>
 class FixShakeKokkos : public FixShake, public KokkosBase {
 
- //friend class FixEHEX;
-
  public:
   typedef DeviceType device_type;
   typedef EV_FLOAT value_type;
@@ -77,7 +75,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
   void shake_end_of_step(int vflag) override;
   void correct_coordinates(int vflag) override;
 
-  int dof(int) override;
+  bigint dof(int) override;
 
   void unconstrained_update() override;
 

src/KOKKOS/pppm_kokkos.cpp

@@ -1381,7 +1381,8 @@ void PPPMKokkos<DeviceType>::poisson_ik()
 
   // global energy and virial contribution
 
-  scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  scaleinv = 1.0/ngridtotal;
   s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {

src/KSPACE/pppm.cpp

@@ -1188,7 +1188,7 @@ double PPPM::compute_qopt()
   // each proc calculates contributions from every Pth grid point
 
   bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
-  int nxy_pppm = nx_pppm * ny_pppm;
+  bigint nxy_pppm = (bigint) nx_pppm * ny_pppm;
 
   double qopt = 0.0;
 
@@ -1944,7 +1944,8 @@ void PPPM::poisson_ik()
 
   // global energy and virial contribution
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {
@@ -2145,7 +2146,8 @@ void PPPM::poisson_ad()
 
   // global energy and virial contribution
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {
@@ -3259,7 +3261,8 @@ void PPPM::poisson_groups(int AA_flag)
   //  keep everything in reciprocal space so
   //  no inverse FFTs needed
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   // energy

src/KSPACE/pppm_dipole.cpp

@@ -1338,7 +1338,8 @@ void PPPMDipole::poisson_ik_dipole()
 
   // global energy and virial contribution
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {

src/KSPACE/pppm_disp.cpp

@@ -4556,7 +4556,8 @@ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
 
   // if requested, compute energy and virial contribution
 
-  double scaleinv = 1.0/(nx_p*ny_p*nz_p);
+  bigint ngridtotal = (bigint) nx_p * ny_p * nz_p;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {
@@ -4696,7 +4697,8 @@ void PPPMDisp::poisson_ad(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
 
   // if requested, compute energy and virial contribution
 
-  double scaleinv = 1.0/(nx_p*ny_p*nz_p);
+  bigint ngridtotal = (bigint) nx_p * ny_p * nz_p;
+  double scaleinv = 1.0/ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {
@@ -4844,7 +4846,8 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
   int i,j,k,n;
   double eng;
 
-  double scaleinv = 1.0/(nx_pppm_6*ny_pppm_6*nz_pppm_6);
+  bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6;
+  double scaleinv = 1.0/ngridtotal;
 
   // transform charge/dispersion density (r -> k)
   // only one transform when energies and pressures not calculated
@@ -5017,7 +5020,8 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
   int i,j,k,n;
   double eng;
 
-  double scaleinv = 1.0/(nx_pppm_6*ny_pppm_6*nz_pppm_6);
+  bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6;
+  double scaleinv = 1.0/ngridtotal;
 
   // transform charge/dispersion density (r -> k)
   // only one transform required when energies and pressures not needed
@@ -5191,7 +5195,8 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
   int i,j,k,n;
   double eng;
 
-  double scaleinv = 1.0/(nx_pppm_6*ny_pppm_6*nz_pppm_6);
+  bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6;
+  double scaleinv = 1.0/ngridtotal;
 
   // transform charge/dispersion density (r -> k)
   // only one tansform required when energies and pressures not needed
@@ -5289,7 +5294,8 @@ poisson_none_ad(int n1, int n2, FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
   int i,j,k,n;
   double eng;
 
-  double scaleinv = 1.0/(nx_pppm_6*ny_pppm_6*nz_pppm_6);
+  bigint ngridtotal = (bigint) nx_pppm_6 * ny_pppm_6 * nz_pppm_6;
+  double scaleinv = 1.0/ngridtotal;
 
   // transform charge/dispersion density (r -> k)
   // only one tansform required when energies and pressures not needed

src/KSPACE/pppm_stagger.cpp

@@ -302,7 +302,7 @@ double PPPMStagger::compute_qopt()
   // each proc calculates contributions from every Pth grid point
 
   bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
-  int nxy_pppm = nx_pppm * ny_pppm;
+  bigint nxy_pppm = (bigint) nx_pppm * ny_pppm;
 
   double qopt = 0.0;
 
@@ -398,7 +398,7 @@ double PPPMStagger::compute_qopt_ad()
   // each proc calculates contributions from every Pth grid point
 
   bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
-  int nxy_pppm = nx_pppm * ny_pppm;
+  bigint nxy_pppm = (bigint) nx_pppm * ny_pppm;
 
   double qopt = 0.0;
 

src/LATBOLTZ/fix_lb_fluid.cpp

@@ -4430,9 +4430,9 @@ void FixLbFluid::calc_MPT(double &totalmass, double totalmomentum[3], double &Ta
 ------------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 
-int FixLbFluid::adjust_dof_fix() /* Based on same private method in compute class */
-{                                /* altered to return fix_dof */
-  int fix_dof = 0;
+bigint FixLbFluid::adjust_dof_fix() /* Based on same private method in compute class */
+{                                   /* altered to return fix_dof */
+  bigint fix_dof = 0;
   for (auto &ifix : modify->get_fix_list())
     if (ifix->dof_flag) fix_dof += ifix->dof(igroup);
   return fix_dof;

src/LATBOLTZ/fix_lb_fluid.h

@@ -182,7 +182,7 @@ class FixLbFluid : public Fix {
   void calc_fluidforceII(void);
   void calc_fluidforceweight(void);
 
-  int adjust_dof_fix();
+  bigint adjust_dof_fix();
   double dof_compute();
 
   /* nanopit parameters */

src/MANIFOLD/fix_nve_manifold_rattle.cpp

@@ -287,21 +287,21 @@ void FixNVEManifoldRattle::update_var_params()
 
 /* -----------------------------------------------------------------------------
    ---------------------------------------------------------------------------*/
-int FixNVEManifoldRattle::dof(int /*igroup*/)
+bigint FixNVEManifoldRattle::dof(int /*igroup*/)
 {
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
-  int natoms = 0;
+  bigint natoms = 0;
   for (int i = 0; i < nlocal; ++i) {
     if (mask[i] & groupbit) ++natoms;
   }
 
-  int dofs;
-  MPI_Allreduce( &natoms, &dofs, 1, MPI_INT, MPI_SUM, world );
+  bigint dofs;
+  MPI_Allreduce( &natoms, &dofs, 1, MPI_LMP_BIGINT, MPI_SUM, world );
 
   // Make sure that, if there is just no or one atom, no dofs are subtracted,
   // since for the first atom already 3 dofs are subtracted because of the
-  // centre of mass corrections:
+  // center of mass corrections:
   if (dofs <= 1) dofs = 0;
   stats.dofs_removed = dofs;
 

src/MANIFOLD/fix_nve_manifold_rattle.h

@@ -75,7 +75,7 @@ class FixNVEManifoldRattle : public Fix {
   void init() override;
   void reset_dt() override;
   void end_of_step() override;
-  int dof(int) override;
+  bigint dof(int) override;
   void setup(int) override {}    // Not needed for fixNVE but is for fixNVT
   double memory_usage() override;