Merge pull request #1769 from akohlmey/collected-small-fixes

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
+the LO value or larger than the HI value, then the calculation is
-force is evaluated as if the bond were the LO or HI length.
+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/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/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->destroy(rho);
-    memory->sfree(fp);
+    memory->destroy(fp);
+    memory->destroy(numforce);
     nmax = atom->nmax;
-    rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho");
+    memory->create(rho,nmax,"pair:rho");
-    fp = (double *) memory->smalloc(nmax*sizeof(double),"pair:fp");
+    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-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;