simplified neighbor list copying to avoid possible same-timestep re-build issues

doc/src/Section_commands.txt

@@ -702,6 +702,8 @@ package"_Section_start.html#start_3.
 "meso"_fix_meso.html,
 "manifoldforce"_fix_manifoldforce.html,
 "meso/stationary"_fix_meso_stationary.html,
+"nve/dot"_fix_nve_dot.html,
+"nve/dotc/langevin"_fix_nve_dotc_langevin.html,
 "nve/manifold/rattle"_fix_nve_manifold_rattle.html,
 "nvk"_fix_nvk.html,
 "nvt/manifold/rattle"_fix_nvt_manifold_rattle.html,
@@ -1035,6 +1037,11 @@ package"_Section_start.html#start_3.
 "morse/soft"_pair_morse.html,
 "multi/lucy"_pair_multi_lucy.html,
 "multi/lucy/rx"_pair_multi_lucy_rx.html,
+"oxdna/coaxstk"_pair_oxdna.html,
+"oxdna/excv"_pair_oxdna.html,
+"oxdna/hbond"_pair_oxdna.html,
+"oxdna/stk"_pair_oxdna.html,
+"oxdna/xstk"_pair_oxdna.html,
 "quip"_pair_quip.html,
 "reax/c (k)"_pair_reax_c.html,
 "smd/hertz"_pair_smd_hertz.html,
@@ -1083,7 +1090,8 @@ if "LAMMPS is built with the appropriate
 package"_Section_start.html#start_3.
 
 "harmonic/shift (o)"_bond_harmonic_shift.html,
-"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html :tb(c=4,ea=c)
+"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html,
+"oxdna/fene"_bond_oxdna_fene.html :tb(c=4,ea=c)
 
 :line
 

doc/src/Section_packages.txt

@@ -84,6 +84,7 @@ Package, Description, Author(s), Doc page, Example, Library
 "PERI"_#PERI, Peridynamics models, Mike Parks (Sandia), "pair_style peri"_pair_peri.html, peri, -
 "POEMS"_#POEMS, coupled rigid body motion, Rudra Mukherjee (JPL), "fix poems"_fix_poems.html, rigid, lib/poems
 "PYTHON"_#PYTHON, embed Python code in an input script, -, "python"_python.html, python, lib/python
+"REAX"_#REAX, ReaxFF potential, Aidan Thompson (Sandia), "pair_style reax"_pair_reax.html, reax, lib/reax
 "REPLICA"_#REPLICA, multi-replica methods, -, "Section 6.6.5"_Section_howto.html#howto_5, tad, -
 "RIGID"_#RIGID, rigid bodies, -, "fix rigid"_fix_rigid.html, rigid, -
 "SHOCK"_#SHOCK, shock loading methods, -, "fix msst"_fix_msst.html, -, -
@@ -1286,27 +1287,26 @@ him directly if you have questions.
 
 USER-CGDNA package :link(USER-CGDNA),h5
 
-Contents: The CGDNA package implements coarse-grained force fields
+Contents: The CGDNA package implements coarse-grained force fields for
-for single- and double-stranded DNA. This is at the moment mainly
+single- and double-stranded DNA. This is at the moment mainly the
-the oxDNA model, developed by Doye, Louis and Ouldridge at the 
+oxDNA model, developed by Doye, Louis and Ouldridge at the University
-University of Oxford. 
+of Oxford.  The package also contains Langevin-type rigid-body
-The package also contains Langevin-type rigid-body integrators
+integrators with improved stability.
-with improved stability.
 
 See these doc pages to get started:
 
 "bond_style oxdna_fene"_bond_oxdna_fene.html
-
 "pair_style oxdna_excv"_pair_oxdna_excv.html
+"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html :ul
 
-"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
+Supporting info: /src/USER-CGDNA/README, "bond_style
-
+oxdna_fene"_bond_oxdna_fene.html, "pair_style
-Supporting info: /src/USER-CGDNA/README, "bond_style oxdna_fene"_bond_oxdna_fene.html,
+oxdna_excv"_pair_oxdna_excv.html, "fix
-"pair_style oxdna_excv"_pair_oxdna_excv.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
+nve/dotc/langevin"_fix_nve_dotc_langevin.html
 
 Author: Oliver Henrich at the University of Edinburgh, UK (o.henrich
-at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk). Contact him directly 
+at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk).  Contact him directly if
-if you have any questions. 
+you have any questions.
 
 :line
 

doc/src/compute_coord_atom.txt

@@ -12,19 +12,16 @@ compute coord/atom command :h3
 
 compute ID group-ID coord/atom cstyle args ... :pre
 
-ID, group-ID are documented in "compute"_compute.html command
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
-coord/atom = style name of this compute command
+coord/atom = style name of this compute command :l
-one cstyle must be appended :ul
+cstyle = {cutoff} or {orientorder} :l
-
-cstyle = {cutoff} or {orientorder}
-
   {cutoff} args = cutoff typeN
     cutoff = distance within which to count coordination neighbors (distance units)
-  typeN = atom type for Nth coordination count (see asterisk form below) :pre
+    typeN = atom type for Nth coordination count (see asterisk form below)
-
   {orientorder} args = orientorderID threshold
-  orientorderID = ID of a previously defined orientorder/atom compute
+    orientorderID = ID of an orientorder/atom compute
-  threshold = minimum value of the scalar product between two 'connected' atoms (see text for explanation) :pre
+    threshold = minimum value of the product of two "connected" atoms :pre
+:ule
 
 [Examples:]
 
@@ -35,21 +32,21 @@ compute 1 all coord/atom orientorder 2 0.5 :pre
 
 [Description:]
 
-This compute performs generic calculations between neighboring atoms. So far,
+This compute performs calculations between neighboring atoms to
-there are two cstyles implemented: {cutoff} and {orientorder}.
+determine a coordination value.  The specific calculation and the
-The {cutoff} cstyle calculates one or more coordination numbers
+meaning of the resulting value depend on the {cstyle} keyword used.
-for each atom in a group.
 
-A coordination number is defined as the number of neighbor atoms with
+The {cutoff} cstyle calculates one or more traditional coordination
-specified atom type(s) that are within the specified cutoff distance
+numbers for each atom.  A coordination number is defined as the number
-from the central atom.  Atoms not in the group are included in a
+of neighbor atoms with specified atom type(s) that are within the
-coordination number of atoms in the group.
+specified cutoff distance from the central atom.  Atoms not in the
+specified group are included in the coordination number tally.
 
-The {typeN} keywords allow you to specify which atom types contribute
+The {typeN} keywords allow specification of which atom types
-to each coordination number.  One coordination number is computed for
+contribute to each coordination number.  One coordination number is
-each of the {typeN} keywords listed.  If no {typeN} keywords are
+computed for each of the {typeN} keywords listed.  If no {typeN}
-listed, a single coordination number is calculated, which includes
+keywords are listed, a single coordination number is calculated, which
-atoms of all types (same as the "*" format, see below).
+includes atoms of all types (same as the "*" format, see below).
 
 The {typeN} keywords can be specified in one of two ways.  An explicit
 numeric value can be used, as in the 2nd example above.  Or a
@@ -61,16 +58,27 @@ from 1 to N.  A leading asterisk means all types from 1 to n
 (inclusive).  A middle asterisk means all types from m to n
 (inclusive).
 
-The {orientorder} cstyle calculates the number of 'connected' atoms j 
+The {orientorder} cstyle calculates the number of "connected" neighbor
-around each atom i. The atom j is connected to i if the scalar product
+atoms J around each central atom I.  For this {cstyle}, connected is
-({Ybar_lm(i)},{Ybar_lm(j)}) is larger than {threshold}. Thus, this cstyle
+defined by the orientational order parameter calculated by the
-will work only if a "compute orientorder/atom"_compute_orientorder_atom.html
+"compute orientorder/atom"_compute_orientorder_atom.html command.
-has been previously defined. This cstyle allows one to apply the 
+This {cstyle} thus allows one to apply the ten Wolde's criterion to
-ten Wolde's criterion to identify cristal-like atoms in a system 
+identify crystal-like atoms in a system, as discussed in "ten
-(see "ten Wolde et al."_#tenWolde).
+Wolde"_#tenWolde.
 
-The value of all coordination numbers will be 0.0 for atoms not in the
+The ID of the previously specified "compute
-specified compute group.
+orientorder/atom"_compute_orientorder/atom command is specified as
+{orientorderID}.  The compute must invoke its {components} option to
+calculate components of the {Ybar_lm} vector for each atoms, as
+described in its documenation.  Note that orientorder/atom compute
+defines its own criteria for identifying neighboring atoms.  If the
+scalar product ({Ybar_lm(i)},{Ybar_lm(j)}), calculated by the
+orientorder/atom compute is larger than the specified {threshold},
+then I and J are connected, and the coordination value of I is
+incremented by one.
+
+For all {cstyle} settings, all coordination values will be 0.0 for
+atoms not in the specified compute group.
 
 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (i.e. each time a snapshot of atoms
@@ -92,21 +100,23 @@ the neighbor list.
 
 [Output info:]
 
-If single {type1} keyword is specified (or if none are specified),
+For {cstyle} cutoff, this compute can calculate a per-atom vector or
-this compute calculates a per-atom vector.  If multiple {typeN}
+array.  If single {type1} keyword is specified (or if none are
-keywords are specified, this compute calculates a per-atom array, with
+specified), this compute calculates a per-atom vector.  If multiple
-N columns.  These values can be accessed by any command that uses
+{typeN} keywords are specified, this compute calculates a per-atom
-per-atom values from a compute as input.  See "Section
+array, with N columns.
+
+For {cstyle} orientorder, this compute calculates a per-atom vector.
+
+These values can be accessed by any command that uses per-atom values
+from a compute as input.  See "Section
 6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output
 options.
 
 The per-atom vector or array values will be a number >= 0.0, as
 explained above.
 
-[Restrictions:]
+[Restrictions:] none
-The cstyle {orientorder} can only be used if a 
-"compute orientorder/atom"_compute_orientorder_atom.html command
-was previously defined. Otherwise, an error message will be issued.
 
 [Related commands:]
 
@@ -118,4 +128,5 @@ was previously defined. Otherwise, an error message will be issued.
 :line
 
 :link(tenWolde)
-[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996).
+[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel,
+J. Chem. Phys. 104, 9932 (1996).

doc/src/compute_orientorder_atom.txt

@@ -19,7 +19,7 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
   {cutoff} value = distance cutoff
   {nnn} value = number of nearest neighbors
   {degrees} values = nlvalues, l1, l2,...
-  {components} value = l  :pre
+  {components} value = ldegree  :pre
 
 :ule
 
@@ -64,21 +64,21 @@ specified distance cutoff are used.
 The optional keyword {degrees} defines the list of order parameters to
 be computed.  The first argument {nlvalues} is the number of order
 parameters. This is followed by that number of integers giving the
-degree of each order parameter. Because {Q}2 and all odd-degree
+degree of each order parameter. Because {Q}2 and all odd-degree order
-order parameters are zero for atoms in cubic crystals
+parameters are zero for atoms in cubic crystals (see
-(see "Steinhardt"_#Steinhardt), the default order parameters
+"Steinhardt"_#Steinhardt), the default order parameters are {Q}4,
-are {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. For the
+{Q}6, {Q}8, {Q}10, and {Q}12. For the FCC crystal with {nnn}=12, {Q}4
-FCC crystal with {nnn}=12, {Q}4 = sqrt(7/3)/8 = 0.19094....
+= sqrt(7/3)/8 = 0.19094....  The numerical values of all order
-The numerical values of all order parameters up to {Q}12
+parameters up to {Q}12 for a range of commonly encountered
-for a range of commonly encountered high-symmetry structures are given
+high-symmetry structures are given in Table I of "Mickel et
-in Table I of "Mickel et al."_#Mickel.
+al."_#Mickel.
 
-The optional keyword {components} will output the components of
+The optional keyword {components} will output the components of the
-the normalized complex vector {Ybar_lm} of degree {l}, which must be
+normalized complex vector {Ybar_lm} of degree {ldegree}, which must be
 explicitly included in the keyword {degrees}. This option can be used
 in conjunction with "compute coord_atom"_compute_coord_atom.html to
-calculate the ten Wolde's criterion to identify crystal-like particles
+calculate the ten Wolde's criterion to identify crystal-like
-(see "ten Wolde et al."_#tenWolde96).
+particles, as discussed in "ten Wolde"_#tenWolde.
 
 The value of {Ql} is set to zero for atoms not in the
 specified compute group, as well as for atoms that have less than
@@ -104,14 +104,16 @@ the neighbor list.
 
 [Output info:]
 
-This compute calculates a per-atom array with {nlvalues} columns, giving the
+This compute calculates a per-atom array with {nlvalues} columns,
-{Ql} values for each atom, which are real numbers on the range 0 <= {Ql} <= 1.
+giving the {Ql} values for each atom, which are real numbers on the
+range 0 <= {Ql} <= 1.
 
-If the keyword {components} is set, then the real and imaginary parts of each
+If the keyword {components} is set, then the real and imaginary parts
-component of (normalized) {Ybar_lm} will be added to the output array in the
+of each component of (normalized) {Ybar_lm} will be added to the
-following order:
+output array in the following order: Re({Ybar_-m}) Im({Ybar_-m})
-Re({Ybar_-m}) Im({Ybar_-m}) Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}).
+Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}).  This
-This way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1) columns.
+way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1)
+columns.
 
 These values can be accessed by any command that uses
 per-atom values from a compute as input.  See "Section
@@ -122,19 +124,25 @@ options.
 
 [Related commands:]
 
-"compute coord/atom"_compute_coord_atom.html, "compute centro/atom"_compute_centro_atom.html, "compute hexorder/atom"_compute_hexorder_atom.html
+"compute coord/atom"_compute_coord_atom.html, "compute
+centro/atom"_compute_centro_atom.html, "compute
+hexorder/atom"_compute_hexorder_atom.html
 
 [Default:]
 
-The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
+The option defaults are {cutoff} = pair style cutoff, {nnn} = 12,
+{degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
 
 :line
 
 :link(Steinhardt)
-[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
+[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti,
+Phys. Rev. B 28, 784 (1983).
 
 :link(Mickel)
-[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke, J. Chem. Phys. 138, 044501 (2013).
+[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke,
+J. Chem. Phys. 138, 044501 (2013).
 
-:link(tenWolde96)
+:link(tenWolde)
-[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996).
+[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel,
+J. Chem. Phys. 104, 9932 (1996).

src/KOKKOS/nbin_kokkos.cpp

@@ -74,10 +74,7 @@ void NBinKokkos<DeviceType>::bin_atoms_setup(int nall)
 
     k_bincount = DAT::tdual_int_1d("Neighbor::d_bincount",mbins);
     bincount = k_bincount.view<DeviceType>();
-    last_bin_memory = update->ntimestep;
   }
-
-  last_bin = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------
@@ -87,6 +84,8 @@ void NBinKokkos<DeviceType>::bin_atoms_setup(int nall)
 template<class DeviceType>
 void NBinKokkos<DeviceType>::bin_atoms()
 {
+  last_bin = update->ntimestep;
+
   h_resize() = 1;
 
   while(h_resize() > 0) {
@@ -116,7 +115,6 @@ void NBinKokkos<DeviceType>::bin_atoms()
       k_bins = DAT::tdual_int_2d("bins", mbins, atoms_per_bin);
       bins = k_bins.view<DeviceType>();
       c_bins = bins;
-      last_bin_memory = update->ntimestep;
     }
   }
 }

src/MC/fix_atom_swap.cpp

@@ -57,7 +57,10 @@ using namespace MathConst;
 
 FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), local_swap_jatom_list(NULL), local_swap_atom_list(NULL), random_equal(NULL), random_unequal(NULL), c_pe(NULL)
+  idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), 
+  sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), 
+  local_swap_jatom_list(NULL), local_swap_atom_list(NULL), 
+  random_equal(NULL), random_unequal(NULL), c_pe(NULL)
 {
   if (narg < 10) error->all(FLERR,"Illegal fix atom/swap command");
 

src/USER-CGDNA/README

@@ -18,20 +18,19 @@ gradient thermostats for rigid body dynamics", J. Chem. Phys. 142,
 144114 (2015).
 
 Example input and data files can be found in
-/examples/USER/cgdna/examples/duplex1/ and /duplex2/.
+/examples/USER/cgdna/examples/duplex1/ and /duplex2/.  A simple python
-A simple python setup tool which creates single straight or helical DNA 
+setup tool which creates single straight or helical DNA strands as
-strands as well as DNA duplexes and arrays of duplexes can be found in 
+well as DNA duplexes and arrays of duplexes can be found in
-/examples/USER/cgdna/util/.
+/examples/USER/cgdna/util/.  A technical report with more information
-A technical report with more information on the model, the structure 
+on the model, the structure of the input and data file, the setup tool
-of the input and data file, the setup tool and the performance of the 
+and the performance of the LAMMPS-implementation of oxDNA can be found
-LAMMPS-implementation of oxDNA can be found in
+in /doc/src/PDF/USER-CGDNA-overview.pdf.
-/doc/src/PDF/USER-CGDNA-overview.pdf.
 
 IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
-with the MOLECULE and ASPHERE packages.  These should be included 
+with the MOLECULE and ASPHERE packages.  These should be included in
-in the LAMMPS build by typing "make yes-asphere yes-molecule" prior 
+the LAMMPS build by typing "make yes-asphere yes-molecule" prior to
-to the usual compilation (see the "Including/excluding packages" 
+the usual compilation (see the "Including/excluding packages" section
-section of the LAMMPS manual).
+of the LAMMPS manual).
 
 The creator of this package is:
 
@@ -39,14 +38,14 @@ Dr Oliver Henrich
 University of Edinburgh, UK
 ohenrich@ph.ed.ac.uk
 o.henrich@epcc.ed.ac.uk
+
 --------------------------------------------------------------------------
 
-Bond styles provided by this package:
+** Bond styles provided by this package:
 
 bond_oxdna_fene.cpp:  backbone connectivity, a modified FENE potential
 
-
+** Pair styles provided by this package:
-Pair styles provided by this package:
 
 pair_oxdna_excv.cpp:  excluded volume interaction between the nucleotides
 
@@ -60,8 +59,7 @@ pair_oxdna_xstk.cpp:  cross-stacking interaction between nucleotides
 
 pair_oxdna_coaxstk.cpp:  coaxial stacking interaction between nucleotides
 
-
+** Fixes provided by this package:
-Fixes provided by this package:
 
 fix_nve_dotc_langevin.cpp:  fix for Langevin-type rigid body integrator "C"
                             in above Ref. [3] 

src/USER-DPD/nbin_ssa.cpp

@@ -18,6 +18,7 @@
 
 #include "nbin_ssa.h"
 #include "atom.h"
+#include "update.h"
 #include "group.h"
 #include "memory.h"
 #include "error.h"
@@ -59,6 +60,8 @@ void NBinSSA::bin_atoms()
   int *mask = atom->mask;
   int *ssaAIR = atom->ssaAIR;
 
+  last_bin = update->ntimestep;
+
   for (i = 0; i < mbins; i++) {
     gbinhead_ssa[i] = -1;
     binhead_ssa[i] = -1;
@@ -126,4 +129,3 @@ bigint NBinSSA::memory_usage()
   }
   return bytes;
 }
-

src/USER-INTEL/nbin_intel.cpp

@@ -86,7 +86,6 @@ void NBinIntel::bin_atoms_setup(int nall)
          nocopy(binhead:length(maxbin+1) alloc_if(1) free_if(0))
     }
     #endif
-    last_bin_memory = update->ntimestep;
   }
 
   // bins = per-atom vector
@@ -127,11 +126,7 @@ void NBinIntel::bin_atoms_setup(int nall)
       _fix->get_single_buffers()->set_bininfo(_atombin,_binpacked);
     else
       _fix->get_double_buffers()->set_bininfo(_atombin,_binpacked);
-
-    last_bin_memory = update->ntimestep;
   }
-
-  last_bin = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------
@@ -140,6 +135,8 @@ void NBinIntel::bin_atoms_setup(int nall)
 
 void NBinIntel::bin_atoms()
 {
+  last_bin = update->ntimestep;
+
   if (_precision_mode == FixIntel::PREC_MODE_MIXED)
     bin_atoms(_fix->get_mixed_buffers());
   else if (_precision_mode == FixIntel::PREC_MODE_SINGLE)

src/USER-INTEL/npair_intel.cpp

@@ -46,9 +46,7 @@ NPairIntel::~NPairIntel() {
   #endif
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   copy needed info from NStencil class to this build class
-------------------------------------------------------------------------- */
 
 #ifdef _LMP_INTEL_OFFLOAD
 void NPairIntel::grow_stencil()

src/compute_coord_atom.cpp

@@ -81,11 +81,11 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
     if (iorientorder < 0)
       error->all(FLERR,"Could not find compute coord/atom compute ID");
     if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
-      error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom");
+      error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom");
 
     threshold = force->numeric(FLERR,arg[5]);
     if (threshold <= -1.0 || threshold >= 1.0)
-      error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
+      error->all(FLERR,"Compute coord/atom threshold not between -1 and 1");
 
     ncol = 1;
     typelo = new int[ncol];
@@ -126,7 +126,7 @@ void ComputeCoordAtom::init()
     comm_forward = 2*(2*l+1);
     if (c_orientorder->iqlcomp < 0)
       error->all(FLERR,"Compute coord/atom requires components "
-                 "option in compute orientorder/atom be defined");
+                 "option in compute orientorder/atom");
   }
 
   if (force->pair == NULL)
@@ -169,9 +169,6 @@ void ComputeCoordAtom::compute_peratom()
 
   invoked_peratom = update->ntimestep;
 
-//  printf("Number of degrees %i components degree %i",nqlist,l);
-//  printf("Particle \t %i \t Norm \t %g \n",0,norm[0][0]);
-
   // grow coordination array if necessary
 
   if (atom->nmax > nmax) {
@@ -195,8 +192,6 @@ void ComputeCoordAtom::compute_peratom()
     }
     nqlist = c_orientorder->nqlist;
     int ltmp = l;
-//    l = c_orientorder->qlcomp;
-    if (ltmp != l) error->all(FLERR,"Debug error, ltmp != l\n");
     normv = c_orientorder->array_atom;
     comm->forward_comm_compute(this);
   }
@@ -340,7 +335,6 @@ void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf)
       normv[i][j] = buf[m++];
     }
   }
-
 }
 
 /* ----------------------------------------------------------------------

src/compute_orientorder_atom.cpp

@@ -73,7 +73,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
   int iarg = 3;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"nnn") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal compute orientorder/atom command");
       if (strcmp(arg[iarg+1],"NULL") == 0) {
         nnn = 0;
       } else {
@@ -91,7 +92,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
       memory->destroy(qlist);
       memory->create(qlist,nqlist,"orientorder/atom:qlist");
       iarg += 2;
-      if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (iarg+nqlist > narg) 
+        error->all(FLERR,"Illegal compute orientorder/atom command");
       qmax = 0;
       for (int iw = 0; iw < nqlist; iw++) {
         qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
@@ -157,11 +159,12 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
 void ComputeOrientOrderAtom::init()
 {
   if (force->pair == NULL)
-    error->all(FLERR,"Compute orientorder/atom requires a pair style be defined");
+    error->all(FLERR,"Compute orientorder/atom requires a "
+               "pair style be defined");
   if (cutsq == 0.0) cutsq = force->pair->cutforce * force->pair->cutforce;
   else if (sqrt(cutsq) > force->pair->cutforce)
-    error->all(FLERR,
+    error->all(FLERR,"Compute orientorder/atom cutoff is "
-               "Compute orientorder/atom cutoff is longer than pairwise cutoff");
+               "longer than pairwise cutoff");
 
   memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r");
   memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i");
@@ -488,7 +491,8 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
    calculate scalar distance
 ------------------------------------------------------------------------- */
 
-double ComputeOrientOrderAtom::dist(const double r[]) {
+double ComputeOrientOrderAtom::dist(const double r[])
+{
   return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
 }
 
@@ -497,8 +501,8 @@ double ComputeOrientOrderAtom::dist(const double r[]) {
    Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi)
 ------------------------------------------------------------------------- */
 
-double ComputeOrientOrderAtom::
+double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta) 
-polar_prefactor(int l, int m, double costheta) {
+{
   const int mabs = abs(m);
 
   double prefactor = 1.0;
@@ -517,8 +521,8 @@ polar_prefactor(int l, int m, double costheta) {
    associated legendre polynomial
 ------------------------------------------------------------------------- */
 
-double ComputeOrientOrderAtom::
+double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) 
-associated_legendre(int l, int m, double x) {
+{
   if (l < m) return 0.0;
 
   double p(1.0), pm1(0.0), pm2(0.0);

src/compute_rdf.cpp

@@ -40,7 +40,8 @@ using namespace MathConst;
 ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
   rdfpair(NULL), nrdfpair(NULL), ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL),
-  hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), duplicates(NULL)
+  hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), 
+  duplicates(NULL)
 {
   if (narg < 4 || (narg-4) % 2) error->all(FLERR,"Illegal compute rdf command");
 

src/nbin.cpp

@@ -24,7 +24,7 @@ using namespace LAMMPS_NS;
 
 NBin::NBin(LAMMPS *lmp) : Pointers(lmp)
 {
-  last_setup = last_bin = last_bin_memory = -1;
+  last_bin = -1;
   maxbin = maxatom = 0;
   binhead = NULL;
   bins = NULL;
@@ -71,7 +71,6 @@ void NBin::bin_atoms_setup(int nall)
     maxbin = mbins;
     memory->destroy(binhead);
     memory->create(binhead,maxbin,"neigh:binhead");
-    last_bin_memory = update->ntimestep;
   }
 
   // bins = per-atom vector
@@ -80,10 +79,7 @@ void NBin::bin_atoms_setup(int nall)
     maxatom = nall;
     memory->destroy(bins);
     memory->create(bins,maxatom,"neigh:bins");
-    last_bin_memory = update->ntimestep;
   }
-
-  last_bin = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------

src/nbin.h

@@ -21,9 +21,7 @@ namespace LAMMPS_NS {
 class NBin : protected Pointers {
  public:
   int istyle;                      // 1-N index into binnames
-
+  bigint last_bin;                 // last timestep atoms were binned
-  bigint last_setup,last_bin;      // timesteps for last operations performed
-  bigint last_bin_memory;
 
   int nbinx,nbiny,nbinz;           // # of global bins
   int mbins;                       // # of local bins and offset on this proc

src/nbin_standard.cpp

@@ -54,8 +54,6 @@ NBinStandard::NBinStandard(LAMMPS *lmp) : NBin(lmp) {}
 
 void NBinStandard::setup_bins(int style)
 {
-  last_setup = update->ntimestep;
-
   // bbox = size of bbox of entire domain
   // bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost
   // for triclinic:
@@ -197,6 +195,7 @@ void NBinStandard::bin_atoms()
 {
   int i,ibin;
 
+  last_bin = update->ntimestep;
   for (i = 0; i < mbins; i++) binhead[i] = -1;
 
   // bin in reverse order so linked list will be in forward order

src/neighbor.cpp

@@ -695,27 +695,38 @@ void Neighbor::init_pair()
   }
       
   // (C) rule: 
-  // for fix/compute requests, occasional or not does not matter
+  // for fix/compute requests to be copied:
   // 1st check:
+  // occasional or not does not matter
+  // Kokkos host/device flags must match
+  // SSA flag must match
+  // newton setting must match (else list has different neighbors in it)
+  // 2nd check:
   // if request = half and non-skip/copy pair half/respaouter request exists,
   // or if request = full and non-skip/copy pair full request exists,
   // or if request = gran and non-skip/copy pair gran request exists,
   // then morph to copy of the matching parent list
-  // 2nd check: only if no match to 1st check
+  // 3rd check: only if no match to 1st check
   // if request = half and non-skip/copy pair full request exists,
   // then morph to half_from_full of the matching parent list
   // for 1st or 2nd check, parent can be copy list or pair or fix
   
+  int inewton,jnewton;
+
   for (i = 0; i < nrequest; i++) {
     if (!requests[i]->fix && !requests[i]->compute) continue;
     for (j = 0; j < nrequest; j++) {
-      // Kokkos flags must match
       if (requests[i]->kokkos_device != requests[j]->kokkos_device) continue;
       if (requests[i]->kokkos_host != requests[j]->kokkos_host) continue;
 
       if (requests[i]->ssa != requests[j]->ssa) continue;
-      // newton 2 and newton 0 both are newton off
+
-      if ((requests[i]->newton & 2) != (requests[j]->newton & 2)) continue;
+      // IJ newton = 1 for newton on, 2 for newton off
+      inewton = requests[i]->newton;
+      if (inewton == 0) inewton = force->newton_pair ? 1 : 2; 
+      jnewton = requests[i]->newton;
+      if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
+      if (inewton != jnewton) continue;
 
       if (requests[i]->half && requests[j]->pair && 
           !requests[j]->skip && requests[j]->half && !requests[j]->copy)
@@ -899,10 +910,12 @@ void Neighbor::init_pair()
   }
 
   // reorder plist vector if necessary
-  // relevant for lists that copy/skip/half-full from parent
+  // relevant for lists that are derived from a parent list:
+  //   half-full,copy,skip
   // the child index must appear in plist after the parent index
   // swap two indices within plist when dependency is mis-ordered
-  // done when entire pass thru plist results in no swaps
+  // start double loop check again whenever a swap is made
+  // done when entire double loop test results in no swaps
 
   NeighList *ptr;
 
@@ -910,11 +923,11 @@ void Neighbor::init_pair()
   while (!done) {
     done = 1;
     for (i = 0; i < npair_perpetual; i++) {
+      for (k = 0; k < 3; k++) {
         ptr = NULL;
-      if (lists[plist[i]]->listfull) ptr = lists[plist[i]]->listfull;
+        if (k == 0) ptr = lists[plist[i]]->listcopy;
-      if (lists[plist[i]]->listcopy) ptr = lists[plist[i]]->listcopy;
+        if (k == 1) ptr = lists[plist[i]]->listskip;
-      // listskip check must be after listfull check
+        if (k == 2) ptr = lists[plist[i]]->listfull;
-      if (lists[plist[i]]->listskip) ptr = lists[plist[i]]->listskip;
         if (ptr == NULL) continue;
         for (m = 0; m < nrequest; m++)
           if (ptr == lists[m]) break;
@@ -927,6 +940,8 @@ void Neighbor::init_pair()
         done = 0;
         break;
       }
+      if (!done) break;
+    }
   }
 
   // debug output
@@ -1841,7 +1856,7 @@ void Neighbor::build_one(class NeighList *mylist, int preflag)
   // create stencil if hasn't been created since last setup_bins() call
 
   NStencil *ns = np->ns;
-  if (ns && ns->last_create < last_setup_bins) {
+  if (ns && ns->last_stencil < last_setup_bins) {
     ns->create_setup();
     ns->create();
   }
@@ -1874,29 +1889,22 @@ void Neighbor::set(int narg, char **arg)
 /* ----------------------------------------------------------------------
    reset timestamps in all NeignBin, NStencil, NPair classes
    so that neighbor lists will rebuild properly with timestep change
+   ditto for lastcall and last_setup_bins
 ------------------------------------------------------------------------- */
 
 void Neighbor::reset_timestep(bigint ntimestep)
 {
-  for (int i = 0; i < nbin; i++) {
+  for (int i = 0; i < nbin; i++)
-    neigh_bin[i]->last_setup = -1;
     neigh_bin[i]->last_bin = -1;
-    neigh_bin[i]->last_bin_memory = -1;
+  for (int i = 0; i < nstencil; i++)
-  }
+    neigh_stencil[i]->last_stencil = -1;
-
-  for (int i = 0; i < nstencil; i++) {
-    neigh_stencil[i]->last_create = -1;
-    neigh_stencil[i]->last_stencil_memory = -1;
-    neigh_stencil[i]->last_copy_bin = -1;
-  }
-
   for (int i = 0; i < nlist; i++) {
     if (!neigh_pair[i]) continue;
     neigh_pair[i]->last_build = -1;
-    neigh_pair[i]->last_copy_bin_setup = -1;
-    neigh_pair[i]->last_copy_bin = -1;
-    neigh_pair[i]->last_copy_stencil = -1;
   }
+
+  lastcall = -1;
+  last_setup_bins = -1;
 }
 
 /* ----------------------------------------------------------------------

src/npair.cpp

@@ -28,8 +28,6 @@ NPair::NPair(LAMMPS *lmp)
   : Pointers(lmp), nb(NULL), ns(NULL), bins(NULL), stencil(NULL)
 {
   last_build = -1;
-  last_copy_bin_setup = last_copy_bin = last_copy_stencil = -1;
-
   molecular = atom->molecular;
 }
 
@@ -76,10 +74,10 @@ void NPair::copy_neighbor_info()
 }
 
 /* ----------------------------------------------------------------------
-   copy bin geometry info from NBin class to this build class
+   copy info from NBin class to this build class
 ------------------------------------------------------------------------- */
 
-void NPair::copy_bin_setup_info()
+void NPair::copy_bin_info()
 {
   nbinx = nb->nbinx;
   nbiny = nb->nbiny;
@@ -95,20 +93,13 @@ void NPair::copy_bin_setup_info()
   bininvx = nb->bininvx;
   bininvy = nb->bininvy;
   bininvz = nb->bininvz;
-}
 
-/* ----------------------------------------------------------------------
-   copy per-atom and per-bin vectors from NBin class to this build class
-------------------------------------------------------------------------- */
-
-void NPair::copy_bin_info()
-{
   bins = nb->bins;
   binhead = nb->binhead;
 }
 
 /* ----------------------------------------------------------------------
-   copy needed info from NStencil class to this build class
+   copy info from NStencil class to this build class
 ------------------------------------------------------------------------- */
 
 void NPair::copy_stencil_info()
@@ -122,23 +113,15 @@ void NPair::copy_stencil_info()
 }
 
 /* ----------------------------------------------------------------------
-   copy needed info from NStencil class to this build class
+   copy info from NBin and NStencil classes to this build class
 ------------------------------------------------------------------------- */
 
 void NPair::build_setup()
 {
-  if (nb && last_copy_bin_setup < nb->last_setup) {
+  if (nb) copy_bin_info();
-    copy_bin_setup_info();
+  if (ns) copy_stencil_info();
-    last_copy_bin_setup = update->ntimestep;
+
-  }
+  // set here, since build_setup() always called before build()
-  if (nb && ((last_copy_bin < nb->last_bin_memory) || (bins != nb->bins))) {
-    copy_bin_info();
-    last_copy_bin = update->ntimestep;
-  }
-  if (ns && ((last_copy_stencil < ns->last_create) || (stencil != ns->stencil))) {
-    copy_stencil_info();
-    last_copy_stencil = update->ntimestep;
-  }
 
   last_build = update->ntimestep;
 }

src/npair.h

@@ -23,11 +23,7 @@ class NPair : protected Pointers {
   int istyle;                   // 1-N index into pairnames
   class NBin *nb;               // ptr to NBin instance I depend on
   class NStencil *ns;           // ptr to NStencil instance I depend on
-
   bigint last_build;            // last timestep build performed
-  bigint last_copy_bin_setup;   // last timestep I invoked copy_bin_setup_info()
-  bigint last_copy_bin;         // last step I invoked copy_bin_info()
-  bigint last_copy_stencil;     // last step I invoked copy_bin_stencil_info()
 
   NPair(class LAMMPS *);
   virtual ~NPair() {}
@@ -93,7 +89,6 @@ class NPair : protected Pointers {
 
   // methods for all NPair variants
 
-  void copy_bin_setup_info();
   virtual void copy_bin_info();
   virtual void copy_stencil_info();
 

src/nstencil.cpp

@@ -56,11 +56,9 @@ enum{NSQ,BIN,MULTI};     // also in Neighbor
 
 NStencil::NStencil(LAMMPS *lmp) : Pointers(lmp)
 {
-  last_create = last_stencil_memory = -1;
+  last_stencil = -1;
-  last_copy_bin = -1;
 
   xyzflag = 0;
-
   maxstencil = maxstencil_multi = 0;
   stencil = NULL;
   stencilxyz = NULL;
@@ -126,10 +124,8 @@ void NStencil::copy_bin_info()
 
 void NStencil::create_setup()
 {
-  if (nb && last_copy_bin < nb->last_setup) {
+  if (nb) copy_bin_info();
-    copy_bin_info();
+  last_stencil = update->ntimestep;
-    last_copy_bin = update->ntimestep;
-  }
 
   // sx,sy,sz = max range of stencil in each dim
   // smax = max possible size of entire 3d stencil
@@ -157,7 +153,6 @@ void NStencil::create_setup()
         memory->destroy(stencilxyz);
         memory->create(stencilxyz,maxstencil,3,"neighstencil:stencilxyz");
       }
-      last_stencil_memory = update->ntimestep;
     }
 
   } else {
@@ -172,7 +167,6 @@ void NStencil::create_setup()
         stencil_multi[i] = NULL;
         distsq_multi[i] = NULL;
       }
-      last_stencil_memory = update->ntimestep;
     }
     if (smax > maxstencil_multi) {
       maxstencil_multi = smax;
@@ -183,12 +177,9 @@ void NStencil::create_setup()
                        "neighstencil:stencil_multi");
         memory->create(distsq_multi[i],maxstencil_multi,
                        "neighstencil:distsq_multi");
-        last_stencil_memory = update->ntimestep;
       }
     }
   }
-
-  last_create = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------

src/nstencil.h

@@ -22,10 +22,7 @@ class NStencil : protected Pointers {
  public:
   int istyle;                      // 1-N index into binnames
   class NBin *nb;                  // ptr to NBin instance I depend on
-
+  bigint last_stencil;             // last timestep stencil was created
-  bigint last_create;              // timesteps for last operations performed
-  bigint last_stencil_memory;
-  bigint last_copy_bin;
 
   int nstencil;                    // # of bins in stencil
   int *stencil;                    // list of bin offsets