Collected small fixes for the next patch release

doc/src/bond_table.rst

@@ -119,8 +119,9 @@ the bond length r (in distance units), the 3rd value is the energy (in
 energy units), and the 4th is the force (in force units).  The bond
 lengths must range from a LO value to a HI value, and increase from
 one line to the next.  If the actual bond length is ever smaller than
-the LO value or larger than the HI value, then the bond energy and
-force is evaluated as if the bond were the LO or HI length.
+the LO value or larger than the HI value, then the calculation is
+aborted with an error, so it is advisable to cover the whole range
+of possible bond lengths.
 
 Note that one file can contain many sections, each with a tabulated
 potential.  LAMMPS reads the file section by section until it finds

src/GRANULAR/fix_pour.cpp

@@ -54,7 +54,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 6) error->all(FLERR,"Illegal fix pour command");
 
-  if (lmp->kokkos) 
+  if (lmp->kokkos)
     error->all(FLERR,"Cannot yet use fix pour with the KOKKOS package");
 
   time_depend = 1;
@@ -797,7 +797,7 @@ bool FixPour::outside(int dim, double value, double lo, double hi)
   bool outside_range = (value < lo || value > hi);
   if (!outside_range || !domain->periodicity[dim]) return outside_range;
 
-  // for periodic dimension: 
+  // for periodic dimension:
   // must perform additional tests if range wraps around the periodic box
 
   bool outside_pbc_range = true;

src/KOKKOS/pair_eam_kokkos.cpp

@@ -37,6 +37,7 @@ template<class DeviceType>
 PairEAMKokkos<DeviceType>::PairEAMKokkos(LAMMPS *lmp) : PairEAM(lmp)
 {
   respa_enable = 0;
+  single_enable = 0;
 
   atomKK = (AtomKokkos *) atom;
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;

src/MANYBODY/pair_eam.cpp

@@ -43,6 +43,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
   nmax = 0;
   rho = NULL;
   fp = NULL;
+  numforce = NULL;
   map = NULL;
   type2frho = NULL;
 
@@ -77,6 +78,7 @@ PairEAM::~PairEAM()
 
   memory->destroy(rho);
   memory->destroy(fp);
+  memory->destroy(numforce);
 
   if (allocated) {
     memory->destroy(setflag);
@@ -151,9 +153,11 @@ void PairEAM::compute(int eflag, int vflag)
   if (atom->nmax > nmax) {
     memory->destroy(rho);
     memory->destroy(fp);
+    memory->destroy(numforce);
     nmax = atom->nmax;
     memory->create(rho,nmax,"pair:rho");
     memory->create(fp,nmax,"pair:fp");
+    memory->create(numforce,nmax,"pair:numforce");
   }
 
   double **x = atom->x;
@@ -255,6 +259,7 @@ void PairEAM::compute(int eflag, int vflag)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    numforce[i] = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -266,6 +271,7 @@ void PairEAM::compute(int eflag, int vflag)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq < cutforcesq) {
+        ++numforce[i];
         jtype = type[j];
         r = sqrt(rsq);
         p = r*rdr + 1.0;
@@ -802,6 +808,18 @@ double PairEAM::single(int i, int j, int itype, int jtype,
   double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip;
   double *coeff;
 
+  if (numforce[i] > 0) {
+    p = rho[i]*rdrho + 1.0;
+    m = static_cast<int> (p);
+    m = MAX(1,MIN(m,nrho-1));
+    p -= m;
+    p = MIN(p,1.0);
+    coeff = frho_spline[type2frho[itype]][m];
+    phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+    if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
+    phi *= 1.0/static_cast<double>(numforce[i]);
+  } else phi = 0.0;
+
   r = sqrt(rsq);
   p = r*rdr + 1.0;
   m = static_cast<int> (p);
@@ -818,7 +836,7 @@ double PairEAM::single(int i, int j, int itype, int jtype,
   z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
 
   recip = 1.0/r;
-  phi = z2*recip;
+  phi += z2*recip;
   phip = z2p*recip - phi*recip;
   psip = fp[i]*rhojp + fp[j]*rhoip + phip;
   fforce = -psip*recip;

src/MANYBODY/pair_eam.h

@@ -70,6 +70,7 @@ class PairEAM : public Pair {
   // per-atom arrays
 
   double *rho,*fp;
+  int *numforce;
 
   // potentials as file data
 

src/MC/fix_gcmc.cpp

@@ -945,7 +945,7 @@ void FixGCMC::attempt_atomic_insertion()
   ninsertion_attempts += 1.0;
 
   if (ngas >= max_ngas) return;
-  
+
   // pick coordinates for insertion point
 
   double coord[3];
@@ -1300,7 +1300,7 @@ void FixGCMC::attempt_molecule_insertion()
   ninsertion_attempts += 1.0;
 
   if (ngas >= max_ngas) return;
-  
+
   double com_coord[3];
   if (regionflag) {
     int region_attempt = 0;
@@ -1634,7 +1634,7 @@ void FixGCMC::attempt_atomic_insertion_full()
   ninsertion_attempts += 1.0;
 
   if (ngas >= max_ngas) return;
-  
+
   double energy_before = energy_stored;
 
   double coord[3];

src/OPT/pair_eam_opt.cpp

@@ -80,11 +80,13 @@ void PairEAMOpt::eval()
   // grow energy array if necessary
 
   if (atom->nmax > nmax) {
-    memory->sfree(rho);
-    memory->sfree(fp);
+    memory->destroy(rho);
+    memory->destroy(fp);
+    memory->destroy(numforce);
     nmax = atom->nmax;
-    rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho");
-    fp = (double *) memory->smalloc(nmax*sizeof(double),"pair:fp");
+    memory->create(rho,nmax,"pair:rho");
+    memory->create(fp,nmax,"pair:fp");
+    memory->create(numforce,nmax,"pair:numforce");
   }
 
   double** _noalias x = atom->x;
@@ -269,6 +271,7 @@ void PairEAMOpt::eval()
 
     fast_gamma_t* _noalias tabssi = &tabss[itype1*ntypes*nr];
     double* _noalias scale_i = scale[itype1+1]+1;
+    numforce[i] = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -280,6 +283,7 @@ void PairEAMOpt::eval()
       double rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq < tmp_cutforcesq) {
+        ++numforce[i];
         jtype = type[j] - 1;
         double r = sqrt(rsq);
         double rhoip,rhojp,z2,z2p;

src/USER-MISC/compute_gyration_shape_chunk.cpp

@@ -39,7 +39,7 @@ ComputeGyrationShapeChunk::ComputeGyrationShapeChunk(LAMMPS *lmp, int narg, char
   int n = strlen(arg[3]) + 1;
   id_gyration_chunk = new char[n];
   strcpy(id_gyration_chunk,arg[3]);
- 
+
   init();
 
   array_flag = 1;

src/USER-OMP/pair_eam_omp.cpp

@@ -51,9 +51,11 @@ void PairEAMOMP::compute(int eflag, int vflag)
   if (atom->nmax > nmax) {
     memory->destroy(rho);
     memory->destroy(fp);
+    memory->destroy(numforce);
     nmax = atom->nmax;
     memory->create(rho,nthreads*nmax,"pair:rho");
     memory->create(fp,nmax,"pair:fp");
+    memory->create(numforce,nmax,"pair:numforce");
   }
 
 #if defined(_OPENMP)
@@ -232,6 +234,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    numforce[i] = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -243,6 +246,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq < cutforcesq) {
+        ++numforce[i];
         jtype = type[j];
         r = sqrt(rsq);
         p = r*rdr + 1.0;

src/fix_print.cpp

@@ -136,7 +136,7 @@ void FixPrint::init()
   } else {
     if (update->ntimestep % nevery)
       next_print = (update->ntimestep/nevery)*nevery + nevery;
-    else 
+    else
       next_print = update->ntimestep;
   }