Merge remote-tracking branch 'github/master' into Ncreate_atoms

# Conflicts:
#	src/create_atoms.cpp

doc/src/Commands_all.rst

@@ -132,10 +132,8 @@ An alphabetic list of all general LAMMPS commands.
    * :doc:`units <units>`
    * :doc:`variable <variable>`
    * :doc:`velocity <velocity>`
-   * :doc:`write\_coeff <write_coeff>`
+   * :doc:`write_coeff <write_coeff>`
-   * :doc:`write\_data <write_data>`
+   * :doc:`write_data <write_data>`
-   * :doc:`write\_dump <write_dump>`
+   * :doc:`write_dump <write_dump>`
-   * :doc:`write\_restart <write_restart>`
+   * :doc:`write_restart <write_restart>`
-   *
+
-   *
-   *

doc/src/Commands_bond.rst

@@ -47,7 +47,6 @@ OPT.
    * :doc:`oxrna2/fene <bond_oxdna>`
    * :doc:`quartic (o) <bond_quartic>`
    * :doc:`table (o) <bond_table>`
-   *
 
 .. _angle:
 
@@ -89,7 +88,6 @@ OPT.
    * :doc:`quartic (o) <angle_quartic>`
    * :doc:`sdk (o) <angle_sdk>`
    * :doc:`table (o) <angle_table>`
-   *
 
 .. _dihedral:
 
@@ -127,8 +125,6 @@ OPT.
    * :doc:`spherical <dihedral_spherical>`
    * :doc:`table (o) <dihedral_table>`
    * :doc:`table/cut <dihedral_table_cut>`
-   *
-   *
 
 .. _improper:
 
@@ -162,4 +158,3 @@ OPT.
    * :doc:`ring (o) <improper_ring>`
    * :doc:`sqdistharm <improper_sqdistharm>`
    * :doc:`umbrella (o) <improper_umbrella>`
-   *

doc/src/Commands_compute.rst

@@ -163,6 +163,4 @@ KOKKOS, o = USER-OMP, t = OPT.
    * :doc:`vcm/chunk <compute_vcm_chunk>`
    * :doc:`voronoi/atom <compute_voronoi_atom>`
    * :doc:`xrd <compute_xrd>`
-   *
+
-   *
-   *

doc/src/Commands_fix.rst

@@ -240,7 +240,3 @@ OPT.
    * :doc:`wall/region <fix_wall_region>`
    * :doc:`wall/region/ees <fix_wall_ees>`
    * :doc:`wall/srd <fix_wall_srd>`
-   *
-   *
-   *
-   *

doc/src/Commands_kspace.rst

@@ -37,4 +37,3 @@ OPT.
    * :doc:`pppm/stagger <kspace_style>`
    * :doc:`pppm/tip4p (o) <kspace_style>`
    * :doc:`scafacos <kspace_style>`
-   *

doc/src/Commands_pair.rst

@@ -257,4 +257,3 @@ OPT.
    * :doc:`yukawa (gko) <pair_yukawa>`
    * :doc:`yukawa/colloid (go) <pair_yukawa_colloid>`
    * :doc:`zbl (gko) <pair_zbl>`
-   *

doc/src/_ext/table_from_list.py

@@ -1,6 +1,6 @@
 from docutils import nodes
 from sphinx.util.docutils import SphinxDirective
-from docutils.nodes import Element, Node
+from docutils.nodes import Element, Node, list_item
 from typing import Any, Dict, List
 from sphinx import addnodes
 from sphinx.util import logging
@@ -26,8 +26,12 @@ class TableFromList(SphinxDirective):
             raise SphinxError('table_from_list content is not a list')
         fulllist = node.children[0]
 
+        # fill list with empty items to have a number of entries
+        # that is divisible by ncolumns
         if (len(fulllist) % ncolumns) != 0:
-            raise SphinxError('number of list elements not a multiple of column number')
+            missing = int(ncolumns - (len(fulllist) % ncolumns))
+            for i in range(0,missing):
+                fulllist += list_item()
 
         table = nodes.table()
         tgroup = nodes.tgroup(cols=ncolumns)

doc/src/pair_lj.rst

@@ -97,7 +97,7 @@ pair\_style lj/cut/tip4p/long command
 =====================================
 
 pair\_style lj/cut/tip4p/long/gpu command
-========================================
+=========================================
 
 pair\_style lj/cut/tip4p/long/omp command
 =========================================

doc/src/read_dump.rst

@@ -144,7 +144,7 @@ i.e. the current directory.
 The *adios* format supports reading data that was written by the 
 :doc:`dump adios <dump_adios>` command. The 
 entire dump is read in parallel across all the processes, dividing
-the atoms evenly amongs the processes. The number of writers that
+the atoms evenly among the processes. The number of writers that
 has written the dump file does not matter. Using the adios style for
 dump and read_dump is a convenient way to dump all atoms from *N* 
 writers and read it back by *M* readers. If one is running two 

python/lammps.py

@@ -1029,10 +1029,19 @@ def get_thermo_data(output):
             r = {'thermo' : thermo_data }
             runs.append(namedtuple('Run', list(r.keys()))(*list(r.values())))
         elif in_run and len(columns) > 0:
-            values = [float(x) for x in line.split()]
+            items = line.split()
+            # Convert thermo output and store it.
+            # It must have the same number of columns and
+            # all of them must be convertable to floats.
+            # Otherwise we ignore the line
+            if len(items) == len(columns):
+                try:
+                    values = [float(x) for x in items]
+                    for i, col in enumerate(columns):
+                        current_run[col].append(values[i])
+                except ValueError:
+                  pass
 
-            for i, col in enumerate(columns):
-                current_run[col].append(values[i])
     return runs
 
 class PyLammps(object):

src/GPU/pair_eam_alloy_gpu.cpp

@@ -189,6 +189,8 @@ void PairEAMAlloyGPU::init_style()
     fp_single = false;
   else
     fp_single = true;
+
+  embedstep = -1;
 }
 
 /* ---------------------------------------------------------------------- */

src/GPU/pair_eam_fs_gpu.cpp

@@ -188,6 +188,8 @@ void PairEAMFSGPU::init_style()
     fp_single = false;
   else
     fp_single = true;
+
+  embedstep = -1;
 }
 
 /* ---------------------------------------------------------------------- */

src/GPU/pair_eam_gpu.cpp

@@ -191,6 +191,8 @@ void PairEAMGPU::init_style()
     fp_single = false;
   else
     fp_single = true;
+
+  embedstep = -1;
 }
 
 /* ---------------------------------------------------------------------- */

src/MANYBODY/pair_eam.cpp

@@ -27,6 +27,7 @@
 #include "neigh_list.h"
 #include "memory.h"
 #include "error.h"
+#include "update.h"
 #include "utils.h"
 
 using namespace LAMMPS_NS;
@@ -39,6 +40,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
 {
   restartinfo = 0;
   manybody_flag = 1;
+  embedstep = -1;
 
   nmax = 0;
   rho = NULL;
@@ -246,6 +248,7 @@ void PairEAM::compute(int eflag, int vflag)
   // communicate derivative of embedding function
 
   comm->forward_comm_pair(this);
+  embedstep = update->ntimestep;
 
   // compute forces on each atom
   // loop over neighbors of my atoms
@@ -423,6 +426,7 @@ void PairEAM::init_style()
   array2spline();
 
   neighbor->request(this,instance_me);
+  embedstep = -1;
 }
 
 /* ----------------------------------------------------------------------
@@ -808,6 +812,14 @@ 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)
+    error->all(FLERR,"EAM embedding data required for this calculation is missing");
+
+  if ((comm->me == 0) && (embedstep != update->ntimestep)) {
+    error->warning(FLERR,"EAM embedding data not computed for this time step ");
+    embedstep = update->ntimestep;
+  }
+
   if (numforce[i] > 0) {
     p = rho[i]*rdrho + 1.0;
     m = static_cast<int> (p);
@@ -916,6 +928,9 @@ void PairEAM::swap_eam(double *fp_caller, double **fp_caller_hold)
   double *tmp = fp;
   fp = fp_caller;
   *fp_caller_hold = tmp;
+
+  // skip warning about out-of-sync timestep, since we already warn in the caller
+  embedstep = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */

src/MANYBODY/pair_eam.h

@@ -66,6 +66,7 @@ class PairEAM : public Pair {
   int nmax;                   // allocated size of per-atom arrays
   double cutforcesq;
   double **scale;
+  bigint embedstep;           // timestep, the embedding term was computed
 
   // per-atom arrays
 

src/RIGID/fix_rigid.cpp

@@ -541,7 +541,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
       p_chain = force->inumeric(FLERR,arg[iarg+1]);
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"inpfile") == 0) {
+    } else if (strcmp(arg[iarg],"infile") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
       delete [] inpfile;
       int n = strlen(arg[iarg+1]) + 1;

src/RIGID/fix_rigid_small.cpp

@@ -216,7 +216,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
       if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
       iarg += 5;
 
-    } else if (strcmp(arg[iarg],"inpfile") == 0) {
+    } else if (strcmp(arg[iarg],"infile") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
       delete [] inpfile;
       int n = strlen(arg[iarg+1]) + 1;

src/USER-MISC/fix_wall_reflect_stochastic.cpp

@@ -198,9 +198,9 @@ FixWallReflectStochastic::~FixWallReflectStochastic()
 void FixWallReflectStochastic::wall_particle(int m, int which, double coord)
 {
   int i, dir, dim, side, sign;
-  double vsave,factor,timecol,difftest,theta;
+  double factor,timecol,difftest,theta;
 
-  double *rmass;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
 
   // coord = current position of wall

src/USER-MISC/pair_e3b.cpp

@@ -539,7 +539,7 @@ void PairE3B::presetParam(const int flag,bool &repeatFlag,double &bondL) {
       bondL!=0.0 || rs!=0.0 || rc3!=0.0 || rc2!=0.0 )
     error->all(FLERR,"Preset keyword will overwrite another keyword setting");
 
-  double econv,lconv;
+  double econv=1.0,lconv=1.0;
   if (strcmp(update->unit_style,"real") == 0) {
     econv=1.0/4.184;
     lconv=1.0;

src/USER-MISC/pair_mesocnt.cpp

@@ -16,12 +16,14 @@
                         pak37@cam.ac.uk
 ------------------------------------------------------------------------- */
 
+#include "pair_mesocnt.h"
+
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
 #include <string>
-#include "pair_mesocnt.h"
+
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
@@ -125,11 +127,7 @@ void PairMesoCNT::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   int **bondlist = neighbor->bondlist;
-  tagint *tag = atom->tag;
-  tagint *mol = atom->molecule;
-  int nlocal = atom->nlocal;
   int nbondlist = neighbor->nbondlist;
-  int newton_pair = force->newton_pair;
 
   // update bond neighbor list when necessary
 
@@ -435,7 +433,7 @@ void PairMesoCNT::allocate()
    global settings
 ------------------------------------------------------------------------- */
 
-void PairMesoCNT::settings(int narg, char **arg)
+void PairMesoCNT::settings(int narg, char ** /* arg */)
 {
   if (narg != 0) error->all(FLERR,"Illegal pair_style command");
 }
@@ -519,7 +517,7 @@ void PairMesoCNT::init_style()
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
 
-double PairMesoCNT::init_one(int i, int j)
+double PairMesoCNT::init_one(int /* i */, int /* j */)
 {
   return cutoff;
 }
@@ -730,7 +728,7 @@ void PairMesoCNT::sort(int *list, int size)
 {
   int i,j,temp1,temp2;
   tagint *tag = atom->tag;
-  for (int i = 1; i < size; i++) {
+  for (i = 1; i < size; i++) {
     j = i;
     temp1 = list[j-1];
     temp2 = list[j];
@@ -1373,10 +1371,6 @@ double PairMesoCNT::dxspline(double x, double y,
   double xbar = x - xlo;
   double ybar = y - ylo;
 
-  double y0 = coeff[i][j][0][0]
-          + ybar*(coeff[i][j][0][1]
-          + ybar*(coeff[i][j][0][2]
-          + ybar*(coeff[i][j][0][3])));
   double y1 = coeff[i][j][1][0]
           + ybar*(coeff[i][j][1][1]
           + ybar*(coeff[i][j][1][2]

src/comm.cpp

@@ -34,6 +34,7 @@
 #include "accelerator_kokkos.h"
 #include "memory.h"
 #include "error.h"
+#include "update.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -688,7 +689,7 @@ double Comm::get_comm_cutoff()
 
   // print warning if neighborlist cutoff overrides user cutoff
 
-  if (me == 0) {
+  if ((me == 0) && (update->setupflag == 1)) {
     if ((cutghostuser > 0.0) && (maxcommcutoff > cutghostuser)) {
       char mesg[128];
       snprintf(mesg,128,"Communication cutoff adjusted to %g",maxcommcutoff);

src/create_atoms.cpp

@@ -802,26 +802,6 @@ void CreateAtoms::add_lattice()
   if (ymin < 0.0) jlo--;
   if (zmin < 0.0) klo--;
 
-  // rough estimate of total time used for create atoms
-  // one inner loop takes about 25ns on a typical desktop CPU core in 2019
-  // maxestimate = time in hours
-
-  double estimate = 2.5e-8/3600.0;
-  estimate *= static_cast<double> (khi-klo+1);
-  estimate *= static_cast<double> (jhi-jlo+1);
-  estimate *= static_cast<double> (ihi-ilo+1);
-  estimate *= static_cast<double> (nbasis);
-
-  double maxestimate = 0.0;
-  MPI_Reduce(&estimate,&maxestimate,1,MPI_DOUBLE,MPI_MAX,0,world);
-
-  if ((comm->me == 0) && (maxestimate > 0.01)) {
-    if (screen) fprintf(screen,"WARNING: create_atoms will take "
-                        "approx. %.2f hours to complete\n",maxestimate);
-    if (logfile) fprintf(logfile,"WARNING: create_atoms will take "
-                         "approx. %.2f hours to complete\n",maxestimate);
-  }
-
   // count lattice sites on each proc
 
   nlatt_overflow = 0;

src/fix_gravity.cpp

@@ -62,7 +62,7 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
     mstyle = CONSTANT;
   }
 
-  int iarg;
+  int iarg=4;
 
   if (strcmp(arg[4],"chute") == 0) {
     if (narg < 6) error->all(FLERR,"Illegal fix gravity command");

src/fix_wall_reflect.cpp

@@ -213,7 +213,7 @@ void FixWallReflect::post_integrate()
    this method may be overwritten by a child class
 ------------------------------------------------------------------------- */
 
-void FixWallReflect::wall_particle(int m, int which, double coord)
+void FixWallReflect::wall_particle(int /* m */, int which, double coord)
 {
   int i,dim,side;
 

src/pair.cpp

@@ -1664,6 +1664,8 @@ void Pair::write_file(int narg, char **arg)
 
   Pair *epair = force->pair_match("^eam",0);
   if (epair) epair->swap_eam(eamfp,&eamfp_hold);
+  if ((comm->me == 0) && (epair))
+    error->warning(FLERR,"EAM pair style. Table will not include embedding term");
 
   // if atom style defines charge, swap in dummy q vec