add support for loading run and min styles

doc/src/plugin.rst

@@ -17,7 +17,7 @@ Syntax
 
      *load* file = load plugin(s) from shared object in *file*
      *unload* style name = unload plugin *name* of style *style*
-         *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *compute* or *fix* or *region* or *command*
+         *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *compute* or *fix* or *region* or *command* or *run* or *min*
      *list* = print a list of currently loaded plugins
      *clear* = unload all currently loaded plugins
 

examples/plugins/CMakeLists.txt

@@ -94,7 +94,10 @@ add_library(zero2plugin MODULE zero2plugin.cpp pair_zero2.cpp bond_zero2.cpp
                                angle_zero2.cpp dihedral_zero2.cpp improper_zero2.cpp)
 target_link_libraries(zero2plugin PRIVATE lammps)
 
-set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin PROPERTIES PREFIX "" SUFFIX ".so")
+add_library(runminplugin MODULE runminplugin.cpp min_cg2.cpp verlet2.cpp)
+target_link_libraries(runminplugin PRIVATE lammps)
+
+set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin runminplugin PROPERTIES PREFIX "" SUFFIX ".so")
 
 # MacOS seems to need this
 if(CMAKE_SYSTEM_NAME STREQUAL Darwin)
@@ -105,13 +108,13 @@ elseif(CMAKE_SYSTEM_NAME STREQUAL "Windows")
   set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin
     PROPERTIES WINDOWS_EXPORT_ALL_SYMBOLS TRUE)
   if(CMAKE_CROSSCOMPILING)
-    set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin
+    set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin runminplugin
       PROPERTIES LINK_FLAGS "-Wl,--export-all-symbols")
   endif()
 else()
-  set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin PROPERTIES
+  set_target_properties(morse2plugin nve2plugin helloplugin zero2plugin runminplugin PROPERTIES
     LINK_FLAGS "-rdynamic")
 endif()
 
 add_custom_target(plugins ALL ${CMAKE_COMMAND} -E echo "Building Plugins"
-  DEPENDS morse2plugin nve2plugin helloplugin zero2plugin morse2plugin)
+  DEPENDS morse2plugin nve2plugin helloplugin zero2plugin morse2plugin runminplugin)

examples/plugins/min_cg2.cpp

@@ -0,0 +1,192 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "min_cg2.h"
+
+#include "error.h"
+#include "output.h"
+#include "timer.h"
+#include "update.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+
+// EPS_ENERGY = minimum normalization for energy tolerance
+
+static constexpr double EPS_ENERGY = 1.0e-8;
+
+/* ---------------------------------------------------------------------- */
+
+MinCG2::MinCG2(LAMMPS *lmp) : MinLineSearch(lmp) {}
+
+/* ----------------------------------------------------------------------
+   minimization via conjugate gradient iterations
+------------------------------------------------------------------------- */
+
+int MinCG2::iterate(int maxiter)
+{
+  int i,m,n,fail,ntimestep;
+  double beta,gg,dot[2],dotall[2],fdotf;
+  double *fatom,*gatom,*hatom;
+
+  // nlimit = max # of CG iterations before restarting
+  // set to ndoftotal unless too big
+
+  int nlimit = static_cast<int> (MIN(MAXSMALLINT,ndoftotal));
+
+  // initialize working vectors
+
+  for (i = 0; i < nvec; i++) h[i] = g[i] = fvec[i];
+  if (nextra_atom)
+    for (m = 0; m < nextra_atom; m++) {
+      fatom = fextra_atom[m];
+      gatom = gextra_atom[m];
+      hatom = hextra_atom[m];
+      n = extra_nlen[m];
+      for (i = 0; i < n; i++) hatom[i] = gatom[i] = fatom[i];
+    }
+  if (nextra_global)
+    for (i = 0; i < nextra_global; i++) hextra[i] = gextra[i] = fextra[i];
+
+  gg = fnorm_sqr();
+
+  for (int iter = 0; iter < maxiter; iter++) {
+
+    if (timer->check_timeout(niter))
+      return TIMEOUT;
+
+    ntimestep = ++update->ntimestep;
+    niter++;
+
+    // line minimization along direction h from current atom->x
+
+    eprevious = ecurrent;
+    fail = (this->*linemin)(ecurrent,alpha_final);
+    if (fail) return fail;
+
+    // function evaluation criterion
+
+    if (neval >= update->max_eval) return MAXEVAL;
+
+    // energy tolerance criterion
+
+    if (fabs(ecurrent-eprevious) <
+        update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
+      return ETOL;
+
+    // force tolerance criterion
+
+    dot[0] = dot[1] = 0.0;
+    for (i = 0; i < nvec; i++) {
+      dot[0] += fvec[i]*fvec[i];
+      dot[1] += fvec[i]*g[i];
+    }
+
+    if (nextra_atom)
+      for (m = 0; m < nextra_atom; m++) {
+        fatom = fextra_atom[m];
+        gatom = gextra_atom[m];
+        n = extra_nlen[m];
+        for (i = 0; i < n; i++) {
+          dot[0] += fatom[i]*fatom[i];
+          dot[1] += fatom[i]*gatom[i];
+        }
+      }
+    MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
+    if (nextra_global)
+      for (i = 0; i < nextra_global; i++) {
+        dotall[0] += fextra[i]*fextra[i];
+        dotall[1] += fextra[i]*gextra[i];
+      }
+
+    fdotf = 0.0;
+    if (update->ftol > 0.0) {
+      if (normstyle == MAX) fdotf = fnorm_max();        // max force norm
+      else if (normstyle == INF) fdotf = fnorm_inf();   // infinite force norm
+      else if (normstyle == TWO) fdotf = dotall[0];     // same as fnorm_sqr(), Euclidean force 2-norm
+      else error->all(FLERR,"Illegal min_modify command");
+      if (fdotf < update->ftol*update->ftol) return FTOL;
+    }
+
+    // update new search direction h from new f = -Grad(x) and old g
+    // this is Polak-Ribieri formulation
+    // beta = dotall[0]/gg would be Fletcher-Reeves
+    // reinitialize CG every ndof iterations by setting beta = 0.0
+
+    beta = MAX(0.0,(dotall[0] - dotall[1])/gg);
+    if ((niter+1) % nlimit == 0) beta = 0.0;
+    gg = dotall[0];
+
+    for (i = 0; i < nvec; i++) {
+      g[i] = fvec[i];
+      h[i] = g[i] + beta*h[i];
+    }
+    if (nextra_atom)
+      for (m = 0; m < nextra_atom; m++) {
+        fatom = fextra_atom[m];
+        gatom = gextra_atom[m];
+        hatom = hextra_atom[m];
+        n = extra_nlen[m];
+        for (i = 0; i < n; i++) {
+          gatom[i] = fatom[i];
+          hatom[i] = gatom[i] + beta*hatom[i];
+        }
+      }
+    if (nextra_global)
+      for (i = 0; i < nextra_global; i++) {
+        gextra[i] = fextra[i];
+        hextra[i] = gextra[i] + beta*hextra[i];
+      }
+
+    // reinitialize CG if new search direction h is not downhill
+
+    dot[0] = 0.0;
+    for (i = 0; i < nvec; i++) dot[0] += g[i]*h[i];
+    if (nextra_atom)
+      for (m = 0; m < nextra_atom; m++) {
+        gatom = gextra_atom[m];
+        hatom = hextra_atom[m];
+        n = extra_nlen[m];
+        for (i = 0; i < n; i++) dot[0] += gatom[i]*hatom[i];
+      }
+    MPI_Allreduce(dot,dotall,1,MPI_DOUBLE,MPI_SUM,world);
+    if (nextra_global)
+      for (i = 0; i < nextra_global; i++)
+        dotall[0] += gextra[i]*hextra[i];
+
+    if (dotall[0] <= 0.0) {
+      for (i = 0; i < nvec; i++) h[i] = g[i];
+      if (nextra_atom)
+        for (m = 0; m < nextra_atom; m++) {
+          gatom = gextra_atom[m];
+          hatom = hextra_atom[m];
+          n = extra_nlen[m];
+          for (i = 0; i < n; i++) hatom[i] = gatom[i];
+        }
+      if (nextra_global)
+        for (i = 0; i < nextra_global; i++) hextra[i] = gextra[i];
+    }
+
+    // output for thermo, dump, restart files
+
+    if (output->next == ntimestep) {
+      timer->stamp();
+      output->write(ntimestep);
+      timer->stamp(Timer::OUTPUT);
+    }
+  }
+
+  return MAXITER;
+}

examples/plugins/min_cg2.h

@@ -0,0 +1,29 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_MIN_CG2_H
+#define LMP_MIN_CG2_H
+
+#include "min_linesearch.h"
+
+namespace LAMMPS_NS {
+
+class MinCG2 : public MinLineSearch {
+ public:
+  MinCG2(class LAMMPS *);
+  int iterate(int) override;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif

examples/plugins/runminplugin.cpp

@@ -0,0 +1,45 @@
+
+#include "lammpsplugin.h"
+
+#include "version.h"
+
+#include <cstring>
+
+#include "min_cg2.h"
+#include "verlet2.h"
+
+using namespace LAMMPS_NS;
+
+static Min *min_cg2creator(LAMMPS *lmp)
+{
+  return new MinCG2(lmp);
+}
+
+static Integrate *verlet2creator(LAMMPS *lmp, int argc, char **argv)
+{
+  return new Verlet2(lmp, argc, argv);
+}
+
+extern "C" void lammpsplugin_init(void *lmp, void *handle, void *regfunc)
+{
+  lammpsplugin_t plugin;
+  lammpsplugin_regfunc register_plugin = (lammpsplugin_regfunc) regfunc;
+
+  plugin.version = LAMMPS_VERSION;
+
+  plugin.style = "min";
+  plugin.name = "cg2";
+  plugin.info = "CG2 minimize style v1.0";
+  plugin.author = "Axel Kohlmeyer (akohlmey@gmail.com)";
+  plugin.creator.v1 = (lammpsplugin_factory1 *) &min_cg2creator;
+  plugin.handle = handle;
+  (*register_plugin)(&plugin, lmp);
+
+  plugin.style = "run";
+  plugin.name = "verlet2";
+  plugin.info = "Verlet2 run style v1.0";
+  plugin.author = "Axel Kohlmeyer (akohlmey@gmail.com)";
+  plugin.creator.v2 = (lammpsplugin_factory2 *) &verlet2creator;
+  plugin.handle = handle;
+  (*register_plugin)(&plugin, lmp);
+}

examples/plugins/verlet2.cpp

@@ -0,0 +1,426 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "verlet2.h"
+
+#include "angle.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "bond.h"
+#include "comm.h"
+#include "dihedral.h"
+#include "domain.h"
+#include "error.h"
+#include "fix.h"
+#include "force.h"
+#include "improper.h"
+#include "kspace.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "output.h"
+#include "pair.h"
+#include "timer.h"
+#include "update.h"
+
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+Verlet2::Verlet2(LAMMPS *lmp, int narg, char **arg) :
+  Integrate(lmp, narg, arg) {}
+
+Verlet2::~Verlet2()
+{
+        // do nothing
+}
+
+/* ----------------------------------------------------------------------
+   initialization before run
+------------------------------------------------------------------------- */
+
+void Verlet2::init()
+{
+  Integrate::init();
+
+  // warn if no fixes doing time integration
+
+  bool do_time_integrate = false;
+  for (const auto &fix : modify->get_fix_list())
+    if (fix->time_integrate) do_time_integrate = true;
+
+  if (!do_time_integrate && (comm->me == 0))
+    error->warning(FLERR,"No fixes with time integration, atoms won't move" + utils::errorurl(28));
+
+  // virial_style:
+  // VIRIAL_PAIR if computed explicitly in pair via sum over pair interactions
+  // VIRIAL_FDOTR if computed implicitly in pair by
+  //   virial_fdotr_compute() via sum over ghosts
+
+  if (force->newton_pair) virial_style = VIRIAL_FDOTR;
+  else virial_style = VIRIAL_PAIR;
+
+  // setup lists of computes for global and per-atom PE and pressure
+
+  ev_setup();
+
+  // detect if fix omp is present for clearing force arrays
+
+  if (modify->get_fix_by_id("package_omp")) external_force_clear = 1;
+
+  // set flags for arrays to clear in force_clear()
+
+  torqueflag = extraflag = 0;
+  if (atom->torque_flag) torqueflag = 1;
+  if (atom->avec->forceclearflag) extraflag = 1;
+
+  // orthogonal vs triclinic simulation box
+
+  triclinic = domain->triclinic;
+}
+
+/* ----------------------------------------------------------------------
+   setup before run
+------------------------------------------------------------------------- */
+
+void Verlet2::setup(int flag)
+{
+  if (comm->me == 0 && screen) {
+    fputs("Setting up Verlet2 run ...\n",screen);
+    if (flag) {
+      utils::print(screen,"  Unit style    : {}\n"
+                        "  Current step  : {}\n"
+                        "  Time step     : {}\n",
+                 update->unit_style,update->ntimestep,update->dt);
+      timer->print_timeout(screen);
+    }
+  }
+
+  if (lmp->kokkos)
+    error->all(FLERR,"KOKKOS package requires run_style verlet2/kk");
+
+  update->setupflag = 1;
+
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+
+  atom->setup();
+  modify->setup_pre_exchange();
+  if (triclinic) domain->x2lamda(atom->nlocal);
+  domain->pbc();
+  domain->reset_box();
+  comm->setup();
+  if (neighbor->style) neighbor->setup_bins();
+  comm->exchange();
+  if (atom->sortfreq > 0) atom->sort();
+  comm->borders();
+  if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+  domain->image_check();
+  domain->box_too_small_check();
+  modify->setup_pre_neighbor();
+  neighbor->build(1);
+  modify->setup_post_neighbor();
+  neighbor->ncalls = 0;
+
+  // compute all forces
+
+  force->setup();
+  ev_set(update->ntimestep);
+  force_clear();
+  modify->setup_pre_force(vflag);
+
+  if (pair_compute_flag) force->pair->compute(eflag,vflag);
+  else if (force->pair) force->pair->compute_dummy(eflag,vflag);
+
+  if (atom->molecular != Atom::ATOMIC) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+  }
+
+  if (force->kspace) {
+    force->kspace->setup();
+    if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
+    else force->kspace->compute_dummy(eflag,vflag);
+  }
+
+  modify->setup_pre_reverse(eflag,vflag);
+  if (force->newton) comm->reverse_comm();
+
+  modify->setup(vflag);
+  output->setup(flag);
+  update->setupflag = 0;
+}
+
+/* ----------------------------------------------------------------------
+   setup without output
+   flag = 0 = just force calculation
+   flag = 1 = reneighbor and force calculation
+------------------------------------------------------------------------- */
+
+void Verlet2::setup_minimal(int flag)
+{
+  update->setupflag = 1;
+
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+
+  if (flag) {
+    modify->setup_pre_exchange();
+    if (triclinic) domain->x2lamda(atom->nlocal);
+    domain->pbc();
+    domain->reset_box();
+    comm->setup();
+    if (neighbor->style) neighbor->setup_bins();
+    comm->exchange();
+    comm->borders();
+    if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+    domain->image_check();
+    domain->box_too_small_check();
+    modify->setup_pre_neighbor();
+    neighbor->build(1);
+    modify->setup_post_neighbor();
+    neighbor->ncalls = 0;
+  }
+
+  // compute all forces
+
+  ev_set(update->ntimestep);
+  force_clear();
+  modify->setup_pre_force(vflag);
+
+  if (pair_compute_flag) force->pair->compute(eflag,vflag);
+  else if (force->pair) force->pair->compute_dummy(eflag,vflag);
+
+  if (atom->molecular != Atom::ATOMIC) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+  }
+
+  if (force->kspace) {
+    force->kspace->setup();
+    if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
+    else force->kspace->compute_dummy(eflag,vflag);
+  }
+
+  modify->setup_pre_reverse(eflag,vflag);
+  if (force->newton) comm->reverse_comm();
+
+  modify->setup(vflag);
+  update->setupflag = 0;
+}
+
+/* ----------------------------------------------------------------------
+   run for N steps
+------------------------------------------------------------------------- */
+
+void Verlet2::run(int n)
+{
+  bigint ntimestep;
+  int nflag,sortflag;
+
+  int n_post_integrate = modify->n_post_integrate;
+  int n_pre_exchange = modify->n_pre_exchange;
+  int n_pre_neighbor = modify->n_pre_neighbor;
+  int n_post_neighbor = modify->n_post_neighbor;
+  int n_pre_force = modify->n_pre_force;
+  int n_pre_reverse = modify->n_pre_reverse;
+  int n_post_force_any = modify->n_post_force_any;
+  int n_end_of_step = modify->n_end_of_step;
+
+  if (atom->sortfreq > 0) sortflag = 1;
+  else sortflag = 0;
+
+  for (int i = 0; i < n; i++) {
+    if (timer->check_timeout(i)) {
+      update->nsteps = i;
+      break;
+    }
+
+    ntimestep = ++update->ntimestep;
+    ev_set(ntimestep);
+
+    // initial time integration
+
+    timer->stamp();
+    modify->initial_integrate(vflag);
+    if (n_post_integrate) modify->post_integrate();
+    timer->stamp(Timer::MODIFY);
+
+    // regular communication vs neighbor list rebuild
+
+    nflag = neighbor->decide();
+
+    if (nflag == 0) {
+      timer->stamp();
+      comm->forward_comm();
+      timer->stamp(Timer::COMM);
+    } else {
+      if (n_pre_exchange) {
+        timer->stamp();
+        modify->pre_exchange();
+        timer->stamp(Timer::MODIFY);
+      }
+      if (triclinic) domain->x2lamda(atom->nlocal);
+      domain->pbc();
+      if (domain->box_change) {
+        domain->reset_box();
+        comm->setup();
+        if (neighbor->style) neighbor->setup_bins();
+      }
+      timer->stamp();
+      comm->exchange();
+      if (sortflag && ntimestep >= atom->nextsort) atom->sort();
+      comm->borders();
+      if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+      timer->stamp(Timer::COMM);
+      if (n_pre_neighbor) {
+        modify->pre_neighbor();
+        timer->stamp(Timer::MODIFY);
+      }
+      neighbor->build(1);
+      timer->stamp(Timer::NEIGH);
+      if (n_post_neighbor) {
+        modify->post_neighbor();
+        timer->stamp(Timer::MODIFY);
+      }
+    }
+
+    // force computations
+    // important for pair to come before bonded contributions
+    // since some bonded potentials tally pairwise energy/virial
+    // and Pair:ev_tally() needs to be called before any tallying
+
+    force_clear();
+
+    timer->stamp();
+
+    if (n_pre_force) {
+      modify->pre_force(vflag);
+      timer->stamp(Timer::MODIFY);
+    }
+
+    if (pair_compute_flag) {
+      force->pair->compute(eflag,vflag);
+      timer->stamp(Timer::PAIR);
+    }
+
+    if (atom->molecular != Atom::ATOMIC) {
+      if (force->bond) force->bond->compute(eflag,vflag);
+      if (force->angle) force->angle->compute(eflag,vflag);
+      if (force->dihedral) force->dihedral->compute(eflag,vflag);
+      if (force->improper) force->improper->compute(eflag,vflag);
+      timer->stamp(Timer::BOND);
+    }
+
+    if (kspace_compute_flag) {
+      force->kspace->compute(eflag,vflag);
+      timer->stamp(Timer::KSPACE);
+    }
+
+    if (n_pre_reverse) {
+      modify->pre_reverse(eflag,vflag);
+      timer->stamp(Timer::MODIFY);
+    }
+
+    // reverse communication of forces
+
+    if (force->newton) {
+      comm->reverse_comm();
+      timer->stamp(Timer::COMM);
+    }
+
+    // force modifications, final time integration, diagnostics
+
+    if (n_post_force_any) modify->post_force(vflag);
+    modify->final_integrate();
+    if (n_end_of_step) modify->end_of_step();
+    timer->stamp(Timer::MODIFY);
+
+    // all output
+
+    if (ntimestep == output->next) {
+      timer->stamp();
+      output->write(ntimestep);
+      timer->stamp(Timer::OUTPUT);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Verlet2::cleanup()
+{
+  modify->post_run();
+  domain->box_too_small_check();
+  update->update_time();
+}
+
+/* ----------------------------------------------------------------------
+   clear force on own & ghost atoms
+   clear other arrays as needed
+------------------------------------------------------------------------- */
+
+void Verlet2::force_clear()
+{
+  size_t nbytes;
+
+  if (external_force_clear) return;
+
+  // clear force on all particles
+  // if either newton flag is set, also include ghosts
+  // when using threads always clear all forces.
+
+  int nlocal = atom->nlocal;
+
+  if (neighbor->includegroup == 0) {
+    nbytes = sizeof(double) * nlocal;
+    if (force->newton) nbytes += sizeof(double) * atom->nghost;
+
+    if (nbytes) {
+      memset(&atom->f[0][0],0,3*nbytes);
+      if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+      if (extraflag) atom->avec->force_clear(0,nbytes);
+    }
+
+  // neighbor includegroup flag is set
+  // clear force only on initial nfirst particles
+  // if either newton flag is set, also include ghosts
+
+  } else {
+    nbytes = sizeof(double) * atom->nfirst;
+
+    if (nbytes) {
+      memset(&atom->f[0][0],0,3*nbytes);
+      if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+      if (extraflag) atom->avec->force_clear(0,nbytes);
+    }
+
+    if (force->newton) {
+      nbytes = sizeof(double) * atom->nghost;
+
+      if (nbytes) {
+        memset(&atom->f[nlocal][0],0,3*nbytes);
+        if (torqueflag) memset(&atom->torque[nlocal][0],0,3*nbytes);
+        if (extraflag) atom->avec->force_clear(nlocal,nbytes);
+      }
+    }
+  }
+}

examples/plugins/verlet2.h

@@ -0,0 +1,39 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_VERLET2_H
+#define LMP_VERLET2_H
+
+#include "integrate.h"
+
+namespace LAMMPS_NS {
+
+class Verlet2 : public Integrate {
+ public:
+  Verlet2(class LAMMPS *, int, char **);
+  ~Verlet2() override;
+  void init() override;
+  void setup(int flag) override;
+  void setup_minimal(int) override;
+  void run(int) override;
+  void force_clear() override;
+  void cleanup() override;
+
+ protected:
+  int triclinic;    // 0 if domain is orthog, 1 if triclinic
+  int torqueflag, extraflag;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif

src/PLUGIN/plugin.cpp

@@ -21,6 +21,7 @@
 #include "force.h"
 #include "input.h"
 #include "modify.h"
+#include "update.h"
 
 #include <cstring>
 #include <list>
@@ -249,6 +250,22 @@ void plugin_register(lammpsplugin_t *plugin, void *ptr)
     }
     (*command_map)[plugin->name] = (Input::CommandCreator) plugin->creator.v1;
 
+  } else if (pstyle == "run") {
+    auto integrate_map = lmp->update->integrate_map;
+    if (integrate_map->find(plugin->name) != integrate_map->end()) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR, "Overriding built-in run style {} from plugin", plugin->name);
+    }
+    (*integrate_map)[plugin->name] = (Update::IntegrateCreator) plugin->creator.v2;
+
+  } else if (pstyle == "min") {
+    auto minimize_map = lmp->update->minimize_map;
+    if (minimize_map->find(plugin->name) != minimize_map->end()) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR, "Overriding built-in run style {} from plugin", plugin->name);
+    }
+    (*minimize_map)[plugin->name] = (Update::MinimizeCreator) plugin->creator.v1;
+
   } else {
     utils::logmesg(lmp, "Loading plugins for {} styles not yet implemented\n", pstyle);
     pluginlist.pop_back();
@@ -272,7 +289,8 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp)
       (strcmp(style, "angle") != 0) && (strcmp(style, "dihedral") != 0) &&
       (strcmp(style, "improper") != 0) && (strcmp(style, "kspace") != 0) &&
       (strcmp(style, "compute") != 0) && (strcmp(style, "fix") != 0) &&
-      (strcmp(style, "region") != 0) && (strcmp(style, "command") != 0)) {
+      (strcmp(style, "region") != 0) && (strcmp(style, "command") != 0) &&
+      (strcmp(style, "run") != 0) && (strcmp(style, "min") != 0)) {
     if (me == 0)
       utils::logmesg(lmp, "Ignoring unload: {} is not a supported plugin style\n", style);
     return;
@@ -395,6 +413,18 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp)
     auto command_map = lmp->input->command_map;
     auto found = command_map->find(name);
     if (found != command_map->end()) command_map->erase(name);
+
+  } else if (pstyle == "run") {
+
+    auto integrate_map = lmp->update->integrate_map;
+    auto found = integrate_map->find(name);
+    if (found != integrate_map->end()) integrate_map->erase(name);
+
+  } else if (pstyle == "min") {
+
+    auto minimize_map = lmp->update->minimize_map;
+    auto found = minimize_map->find(name);
+    if (found != minimize_map->end()) minimize_map->erase(name);
   }
 
   // if reference count is down to zero, close DSO handle.
@@ -410,7 +440,7 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp)
 
 void plugin_clear(LAMMPS *lmp)
 {
-  verbose = false;
+  verbose = true;
   while (pluginlist.size() > 0) {
     auto p = pluginlist.begin();
     plugin_unload(p->style, p->name, lmp);