From 747ea7246185aa5ec372d2cfa972609bb1f37566 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Wed, 12 Oct 2011 22:20:00 -0400
Subject: [PATCH] force arrays and alike are cleared in parallel when the omp
 package command is active. FixOMP does no longer need to be a friend.

---
 src/integrate.cpp | 1 +
 src/integrate.h   | 2 +-
 src/min.cpp       | 9 ++++++++-
 src/min.h         | 2 +-
 src/respa.cpp     | 8 +++++++-
 src/verlet.cpp    | 9 +++++++--
 6 files changed, 25 insertions(+), 6 deletions(-)

diff --git a/src/integrate.cpp b/src/integrate.cpp
index 80bf89f676..6de4df26c3 100644
--- a/src/integrate.cpp
+++ b/src/integrate.cpp
@@ -26,6 +26,7 @@ Integrate::Integrate(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 {
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
+  external_force_clear = 0;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/integrate.h b/src/integrate.h
index 0950c82148..3c75176a10 100644
--- a/src/integrate.h
+++ b/src/integrate.h
@@ -19,7 +19,6 @@
 namespace LAMMPS_NS {
 
 class Integrate : protected Pointers {
- friend class FixOMP;
  public:
   Integrate(class LAMMPS *, int, char **);
   virtual ~Integrate();
@@ -34,6 +33,7 @@ class Integrate : protected Pointers {
  protected:
   int eflag,vflag;                  // flags for energy/virial computation
   int virial_style;                 // compute virial explicitly or implicitly
+  int external_force_clear;         // clear forces locally or externally
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;
diff --git a/src/min.cpp b/src/min.cpp
index 01a835aa40..f35ceed957 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -57,6 +57,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
+  external_force_clear = 0;
 
   nextra_global = 0;
   fextra = NULL;
@@ -133,6 +134,10 @@ void Min::init()
 
   ev_setup();
 
+  // detect if fix omp is present and will clear force arrays for us
+  int ifix = modify->find_fix("package_omp");
+  if (ifix >= 0) external_force_clear = 1;
+
   // set flags for what arrays to clear in force_clear()
   // need to clear additionals arrays if they exist
 
@@ -520,6 +525,8 @@ double Min::energy_force(int resetflag)
 
 void Min::force_clear()
 {
+  if (external_force_clear) return;
+
   int i;
 
   // clear global force array
@@ -529,7 +536,7 @@ void Min::force_clear()
   if (force->newton) nall = atom->nlocal + atom->nghost;
   else nall = atom->nlocal;
 
-  size_t nbytes = sizeof(double) * nall * comm->nthreads;
+  size_t nbytes = sizeof(double) * nall;
 
   memset(&(atom->f[0][0]),0,3*nbytes);
   if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);
diff --git a/src/min.h b/src/min.h
index c05c58fe56..47956059d8 100644
--- a/src/min.h
+++ b/src/min.h
@@ -19,7 +19,6 @@
 namespace LAMMPS_NS {
 
 class Min : protected Pointers {
- friend class FixOMP;
  public:
   double einitial,efinal,eprevious;
   double fnorm2_init,fnorminf_init,fnorm2_final,fnorminf_final;
@@ -50,6 +49,7 @@ class Min : protected Pointers {
  protected:
   int eflag,vflag;            // flags for energy/virial computation
   int virial_style;           // compute virial explicitly or implicitly
+  int external_force_clear;   // clear forces locally or externally
 
   double dmax;                // max dist to move any atom in one step
   int linestyle;              // 0 = backtrack, 1 = quadratic
diff --git a/src/respa.cpp b/src/respa.cpp
index f247118b20..6ef711e2bb 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -277,6 +277,10 @@ void Respa::init()
 
   ev_setup();
 
+  // detect if fix omp is present and will clear force arrays for us
+  int ifix = modify->find_fix("package_omp");
+  if (ifix >= 0) external_force_clear = 1;
+
   // set flags for what arrays to clear in force_clear()
   // need to clear additionals arrays if they exist
 
@@ -615,6 +619,8 @@ void Respa::recurse(int ilevel)
 
 void Respa::force_clear(int newtonflag)
 {
+  if (external_force_clear) return;
+
   int i;
 
   // clear global force array
@@ -624,7 +630,7 @@ void Respa::force_clear(int newtonflag)
   if (newtonflag) nall = atom->nlocal + atom->nghost;
   else nall = atom->nlocal;
 
-  size_t nbytes = sizeof(double) * nall * comm->nthreads;
+  size_t nbytes = sizeof(double) * nall;
 
   memset(&(atom->f[0][0]),0,3*nbytes);
   if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 284288e88a..0023fe339b 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -62,6 +62,10 @@ void Verlet::init()
 
   ev_setup();
 
+  // detect if fix omp is present and will clear force arrays for us
+  int ifix = modify->find_fix("package_omp");
+  if (ifix >= 0) external_force_clear = 1;
+
   // set flags for what arrays to clear in force_clear()
   // need to clear additionals arrays if they exist
 
@@ -321,18 +325,19 @@ void Verlet::cleanup()
 
 void Verlet::force_clear()
 {
+  if (external_force_clear) return;
   int i;
 
   // clear force on all particles
   // if either newton flag is set, also include ghosts
   // when using threads always clear all forces.
 
-  if (neighbor->includegroup == 0 || comm->nthreads > 1) {
+  if (neighbor->includegroup == 0) {
     int nall;
     if (force->newton) nall = atom->nlocal + atom->nghost;
     else nall = atom->nlocal;
 
-    size_t nbytes = sizeof(double) * nall * comm->nthreads;
+    size_t nbytes = sizeof(double) * nall;
 
     memset(&(atom->f[0][0]),0,3*nbytes);
     if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);