git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8434 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/USER-CG-CMM/Package.sh

@@ -15,10 +15,10 @@ for file in *.cpp *.h; do
   if (test $file = pair_cg_cmm_coul_long.h -a ! -e ../pair_lj_cut_coul_long.h) then
     continue
   fi
-  if (test $file = angle_sdk.cpp -a ! -e ../pair_angle_harmonic.cpp) then
+  if (test $file = angle_sdk.cpp -a ! -e ../angle_harmonic.cpp) then
     continue
   fi
-  if (test $file = angle_sdk.h -a ! -e ../pair_angle_harmonic.h) then
+  if (test $file = angle_sdk.h -a ! -e ../angle_harmonic.h) then
     continue
   fi
   if (test $file = pair_lj_sdk_coul_long.cpp -a ! -e ../pair_lj_cut_coul_long.cpp) then

src/USER-OMP/fix_qeq_comb_omp.cpp

@@ -19,11 +19,17 @@
 #include "mpi.h"
 #include "math.h"
 #include "fix_qeq_comb_omp.h"
+#include "fix_omp.h"
 #include "atom.h"
+#include "comm.h"
 #include "force.h"
 #include "group.h"
 #include "memory.h"
+#include "modify.h"
 #include "error.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
 #include "respa.h"
 #include "update.h"
 #include "pair_comb_omp.h"
@@ -58,13 +64,28 @@ void FixQEQCombOMP::init()
 
   ngroup = group->count(igroup);
   if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
+
+  // determine status of neighbor flag of the omp package command
+  int ifix = modify->find_fix("package_omp");
+  int use_omp = 0;
+  if (ifix >=0) {
+     FixOMP * fix = static_cast<FixOMP *>(lmp->modify->fix[ifix]);
+     if (fix->get_neighbor()) use_omp = 1;
+  }
+
+  int irequest = neighbor->request(this);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->fix = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->omp = use_omp;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixQEQCombOMP::post_force(int vflag)
 {
-  int i,iloop,loopmax;
+  int i,ii,iloop,loopmax,inum,*ilist;
   double heatpq,qmass,dtq,dtq2;
   double enegchkall,enegmaxall;
 
@@ -89,8 +110,8 @@ void FixQEQCombOMP::post_force(int vflag)
   // more loops for first-time charge equilibrium
 
   iloop = 0; 
-  if (firstflag) loopmax = 5000;
+  if (firstflag) loopmax = 500;
-  else loopmax = 2000;
+  else loopmax = 200;
 
   // charge-equilibration loop
 
@@ -99,8 +120,8 @@ void FixQEQCombOMP::post_force(int vflag)
             update->ntimestep);
 
   heatpq = 0.05;
-  qmass  = 0.000548580;
+  qmass  = 0.016;
-  dtq    = 0.0006;
+  dtq    = 0.01;
   dtq2   = 0.5*dtq*dtq/qmass;
 
   double enegchk = 0.0;
@@ -109,29 +130,38 @@ void FixQEQCombOMP::post_force(int vflag)
 
   double *q = atom->q;
   int *mask = atom->mask;
-  int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++)
+  inum = comb->list->inum;
+  ilist = comb->list->ilist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
     q1[i] = q2[i] = qf[i] = 0.0;
+  }
 
   for (iloop = 0; iloop < loopmax; iloop ++ ) {
-    for (i = 0; i < nlocal; i++)
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
       if (mask[i] & groupbit) {
         q1[i] += qf[i]*dtq2 - heatpq*q1[i];
         q[i]  += q1[i]; 
       }
+    }
+    comm->forward_comm_fix(this);
 
-    enegtot = comb->yasu_char(qf,igroup);
+    if(comb) enegtot = comb->yasu_char(qf,igroup);
     enegtot /= ngroup;
     enegchk = enegmax = 0.0;
     
-    for (i = 0; i < nlocal ; i++)
+    for (ii = 0; ii < inum ; ii++) {
+      i = ilist[ii];
       if (mask[i] & groupbit) {
         q2[i] = enegtot-qf[i];
         enegmax = MAX(enegmax,fabs(q2[i]));
         enegchk += fabs(q2[i]);
         qf[i] = q2[i];
       }
+    }
 
     MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
     enegchk = enegchkall/ngroup;
@@ -145,10 +175,12 @@ void FixQEQCombOMP::post_force(int vflag)
               "enegmax %.6g, fq deviation: %.6g\n",
               iloop,enegtot,enegmax,enegchk); 
     
-    for (i = 0; i < nlocal; i++)
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
       if (mask[i] & groupbit)
         q1[i] += qf[i]*dtq2 - heatpq*q1[i]; 
     }
+  } 
   
   if (me == 0 && fp) {
     if (iloop == loopmax)

src/USER-OMP/neigh_full_omp.cpp

@@ -178,6 +178,7 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list)
 
     // loop over all atoms, owned and ghost
     // skip i = j
+    // no molecular test when i = ghost atom
 
     if (i < nlocal) {
       for (j = 0; j < nall; j++) {
@@ -210,8 +211,7 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list)
         delz = ztmp - x[j][2];
         rsq = delx*delx + dely*dely + delz*delz;
 
-        if (rsq <= cutneighghostsq[itype][jtype])
+        if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
-          neighptr[n++] = j;
       }
     }
 
@@ -401,6 +401,7 @@ void Neighbor::full_bin_ghost_omp(NeighList *list)
     // loop over all atoms in surrounding bins in stencil including self
     // when i is a ghost atom, must check if stencil bin is out of bounds
     // skip i = j
+    // no molecular test when i = ghost atom
 
     if (i < nlocal) {
       ibin = coord2bin(x[i]);
@@ -448,8 +449,7 @@ void Neighbor::full_bin_ghost_omp(NeighList *list)
           delz = ztmp - x[j][2];
           rsq = delx*delx + dely*dely + delz*delz;
 
-          if (rsq <= cutneighghostsq[itype][jtype])
+          if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
-            neighptr[n++] = j;
         }
       }
     }

src/USER-OMP/neigh_half_bin_omp.cpp

@@ -129,6 +129,149 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list)
   list->inum = nlocal;
 }
 
+/* ----------------------------------------------------------------------
+   binned neighbor list construction with partial Newton's 3rd law
+   include neighbors of ghost atoms, but no "special neighbors" for ghosts
+   owned and ghost atoms check own bin and other bins in stencil
+   pair stored once if i,j are both owned and i < j
+   pair stored by me if i owned and j ghost (also stored by proc owning j)
+   pair stored once if i,j are both ghost and i < j
+------------------------------------------------------------------------- */
+
+void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list)
+{
+  const int nlocal = atom->nlocal;
+  const int nall = nlocal + atom->nghost;
+
+  NEIGH_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(list)
+#endif
+  NEIGH_OMP_SETUP(nall);
+
+  int i,j,k,n,itype,jtype,ibin,which;
+  int xbin,ybin,zbin,xbin2,ybin2,zbin2;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr;
+
+  // bin local & ghost atoms
+
+  bin_atoms();
+
+  // loop over each atom, storing neighbors
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int molecular = atom->molecular;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  int nstencil = list->nstencil;
+  int *stencil = list->stencil;
+  int **stencilxyz = list->stencilxyz;
+
+  // each thread works on its own page
+  int npage = tid;
+  int npnt = 0;
+
+  for (i = ifrom; i < ito; i++) {
+
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
+    if (pgsize - npnt < oneatom) {
+      npnt = 0;
+      npage += nthreads;
+      if (npage >= list->maxpage) list->add_pages(nthreads);
+    }
+
+    neighptr = &(list->pages[npage][npnt]);
+    n = 0;
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    // loop over all atoms in other bins in stencil including self
+    // when i is a ghost atom, must check if stencil bin is out of bounds
+    // only store pair if i < j
+    // stores own/own pairs only once
+    // stores own/ghost pairs with owned atom only, on both procs
+    // stores ghost/ghost pairs only once
+    // no molecular test when i = ghost atom
+
+    if (i < nlocal) {
+      ibin = coord2bin(x[i]);
+
+      for (k = 0; k < nstencil; k++) {
+        for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
+          if (j <= i) continue;
+
+          jtype = type[j];
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          
+          if (rsq <= cutneighsq[itype][jtype]) {
+            if (molecular) {
+              which = find_special(special[i],nspecial[i],tag[j]);
+              if (which == 0) neighptr[n++] = j;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = j;
+              else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = j;
+          }
+        }
+      }
+
+    } else {
+      ibin = coord2bin(x[i],xbin,ybin,zbin);
+      for (k = 0; k < nstencil; k++) {
+        xbin2 = xbin + stencilxyz[k][0];
+        ybin2 = ybin + stencilxyz[k][1];
+        zbin2 = zbin + stencilxyz[k][2];
+        if (xbin2 < 0 || xbin2 >= mbinx ||
+            ybin2 < 0 || ybin2 >= mbiny ||
+            zbin2 < 0 || zbin2 >= mbinz) continue;
+        for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
+          if (j <= i) continue;
+
+          jtype = type[j];
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+
+          if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
+        }
+      }
+    }
+
+    ilist[i] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (n > oneatom)
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+  NEIGH_OMP_CLOSE;
+  list->inum = nlocal;
+  list->gnum = nall - atom->nlocal;
+}
+
 /* ----------------------------------------------------------------------
    binned neighbor list construction with full Newton's 3rd law
    each owned atom i checks its own bin and other bins in Newton stencil

src/USER-OMP/neigh_half_nsq_omp.cpp

@@ -32,6 +32,7 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
   const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+  const int nall = atom->nlocal + atom->nghost;
 
   NEIGH_OMP_INIT;
 #if defined(_OPENMP)
@@ -43,8 +44,6 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
-  // loop over each atom, storing neighbors
-
   int **special = atom->special;
   int **nspecial = atom->nspecial;
   int *tag = atom->tag;
@@ -53,16 +52,18 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
   int *type = atom->type;
   int *mask = atom->mask;
   int *molecule = atom->molecule;
-  int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
 
+  // each thread works on its own page
   int npage = tid;
   int npnt = 0;
 
+  // loop over owned atoms, storing neighbors
+
   for (i = ifrom; i < ito; i++) {
 
 #if defined(_OPENMP)
@@ -83,6 +84,7 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
     ztmp = x[i][2];
 
     // loop over remaining atoms, owned and ghost
+    // only store pair if i < j
 
     for (j = i+1; j < nall; j++) {
       if (includegroup && !(mask[j] & bitmask)) continue;
@@ -116,6 +118,125 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
   list->inum = nlocal;
 }
 
+/* ----------------------------------------------------------------------
+   N^2 / 2 search for neighbor pairs with partial Newton's 3rd law
+   include neighbors of ghost atoms, but no "special neighbors" for ghosts
+   pair stored once if i,j are both owned and i < j
+   pair stored by me if i owned and j ghost (also stored by proc owning j)
+   pair stored once if i,j are both ghost and i < j
+------------------------------------------------------------------------- */
+
+void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list)
+{
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+  const int nall = nlocal + atom->nghost;
+
+  NEIGH_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(list)
+#endif
+  NEIGH_OMP_SETUP(nall);
+
+  int i,j,n,itype,jtype,which;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr;
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int molecular = atom->molecular;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+
+  // each thread works on its own page
+  int npage = tid;
+  int npnt = 0;
+
+  // loop over owned & ghost atoms, storing neighbors
+
+  for (i = ifrom; i < ito; i++) {
+
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
+    if (pgsize - npnt < oneatom) {
+      npnt = 0;
+      npage += nthreads;
+      if (npage >= list->maxpage) list->add_pages(nthreads);
+    }
+
+    neighptr = &(list->pages[npage][npnt]);
+    n = 0;
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    // loop over remaining atoms, owned and ghost
+    // only store pair if i < j
+    // stores own/own pairs only once
+    // stores own/ghost pairs with owned atom only, on both procs
+    // stores ghost/ghost pairs only once
+    // no molecular test when i = ghost atom
+
+    if (i < nlocal) {
+      for (j = i+1; j < nall; j++) {
+        if (includegroup && !(mask[j] & bitmask)) continue;
+        jtype = type[j];
+        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+        
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx*delx + dely*dely + delz*delz;
+        
+        if (rsq <= cutneighsq[itype][jtype]) {
+          if (molecular) {
+            which = find_special(special[i],nspecial[i],tag[j]);
+            if (which == 0) neighptr[n++] = j;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = j;
+            else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = j;
+        }
+      }
+
+    } else {
+      for (j = i+1; j < nall; j++) {
+        if (includegroup && !(mask[j] & bitmask)) continue;
+        jtype = type[j];
+        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+        
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx*delx + dely*dely + delz*delz;
+        
+        if (rsq <= cutneighsq[itype][jtype]) neighptr[n++] = j;
+      }
+    }
+    
+    ilist[i] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (n > oneatom)
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+  NEIGH_OMP_CLOSE;
+  list->inum = atom->nlocal;
+  list->gnum = nall - atom->nlocal;
+}
+
 /* ----------------------------------------------------------------------
    N^2 / 2 search for neighbor pairs with full Newton's 3rd law
    every pair stored exactly once by some processor

src/USER-OMP/pair_gran_hertz_history_omp.cpp

@@ -16,6 +16,7 @@
 #include "pair_gran_hertz_history_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "fix_rigid.h"
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -48,6 +49,13 @@ void PairGranHertzHistoryOMP::compute(int eflag, int vflag)
   computeflag = 1;
   const int shearupdate = (update->setupflag) ? 0 : 1;
 
+  // update body ptr and values for ghost atoms if using FixRigid masses
+
+  if (fix_rigid && neighbor->ago == 0) {
+    body = fix_rigid->body;
+    comm->forward_comm_pair(this);
+  }
+
 #if defined(_OPENMP)
 #pragma omp parallel default(none) shared(eflag,vflag)
 #endif
@@ -78,7 +86,7 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
-  double meff,damp,ccel,tor1,tor2,tor3;
+  double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
   double fn,fs,fs1,fs2,fs3;
   double shrmag,rsht,polyhertz;
   int *ilist,*jlist,*numneigh,**firstneigh;
@@ -171,19 +179,27 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
         wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
         wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
 
-        // normal force = Hertzian contact + normal velocity damping
+        // meff = effective mass of pair of particles
+        // if I or J part of rigid body, use body mass
+        // if I or J is frozen, meff is other particle
 
         if (rmass) {
-          meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+          mi = rmass[i];
-          if (mask[i] & freeze_group_bit) meff = rmass[j];
+          mj = rmass[j];
-          if (mask[j] & freeze_group_bit) meff = rmass[i];
         } else {
-          itype = type[i];
+          mi = mass[type[i]];
-          jtype = type[j];
+          mj = mass[type[j]];
-          meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
-          if (mask[i] & freeze_group_bit) meff = mass[jtype];
-          if (mask[j] & freeze_group_bit) meff = mass[itype];
         }
+        if (fix_rigid) {
+          if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
+          if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
+        }
+
+        meff = mi*mj / (mi+mj);
+        if (mask[i] & freeze_group_bit) meff = mj;
+        if (mask[j] & freeze_group_bit) meff = mi;
+
+        // normal force = Hertzian contact + normal velocity damping
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;

src/USER-OMP/pair_gran_hooke_history_omp.cpp

@@ -16,6 +16,7 @@
 #include "pair_gran_hooke_history_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "fix_rigid.h"
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -49,6 +50,13 @@ void PairGranHookeHistoryOMP::compute(int eflag, int vflag)
   computeflag = 1;
   const int shearupdate = (update->setupflag) ? 0 : 1;
 
+  // update body ptr and values for ghost atoms if using FixRigid masses
+
+  if (fix_rigid && neighbor->ago == 0) {
+    body = fix_rigid->body;
+    comm->forward_comm_pair(this);
+  }
+
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
   const int inum = list->inum;
@@ -84,7 +92,7 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
-  double meff,damp,ccel,tor1,tor2,tor3;
+  double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
   double fn,fs,fs1,fs2,fs3;
   double shrmag,rsht;
   int *ilist,*jlist,*numneigh,**firstneigh;
@@ -176,18 +184,27 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
         wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
         wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
 
-        // normal forces = Hookian contact + normal velocity damping
+        // meff = effective mass of pair of particles
+        // if I or J part of rigid body, use body mass
+        // if I or J is frozen, meff is other particle
 
         if (rmass) {
-          meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+          mi = rmass[i];
-          if (mask[i] & freeze_group_bit) meff = rmass[j];
+          mj = rmass[j];
-          if (mask[j] & freeze_group_bit) meff = rmass[i];
         } else {
-          jtype = type[j];
+          mi = mass[type[i]];
-          meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+          mj = mass[type[j]];
-          if (mask[i] & freeze_group_bit) meff = mass[jtype];
-          if (mask[j] & freeze_group_bit) meff = mass[itype];
         }
+        if (fix_rigid) {
+          if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
+          if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
+        }
+
+        meff = mi*mj / (mi+mj);
+        if (mask[i] & freeze_group_bit) meff = mj;
+        if (mask[j] & freeze_group_bit) meff = mi;
+
+        // normal forces = Hookian contact + normal velocity damping
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;

src/USER-OMP/pair_gran_hooke_omp.cpp

@@ -16,6 +16,7 @@
 #include "pair_gran_hooke_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "fix_rigid.h"
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -44,6 +45,13 @@ void PairGranHookeOMP::compute(int eflag, int vflag)
   const int nthreads = comm->nthreads;
   const int inum = list->inum;
 
+  // update body ptr and values for ghost atoms if using FixRigid masses
+
+  if (fix_rigid && neighbor->ago == 0) {
+    body = fix_rigid->body;
+    comm->forward_comm_pair(this);
+  }
+
 #if defined(_OPENMP)
 #pragma omp parallel default(none) shared(eflag,vflag)
 #endif
@@ -74,7 +82,7 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr)
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
-  double meff,damp,ccel,tor1,tor2,tor3;
+  double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
   double fn,fs,ft,fs1,fs2,fs3;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
@@ -150,19 +158,27 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr)
         wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
         wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
 
-        // normal forces = Hookian contact + normal velocity damping
+        // meff = effective mass of pair of particles
+        // if I or J part of rigid body, use body mass
+        // if I or J is frozen, meff is other particle
 
         if (rmass) {
-          meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+          mi = rmass[i];
-          if (mask[i] & freeze_group_bit) meff = rmass[j];
+          mj = rmass[j];
-          if (mask[j] & freeze_group_bit) meff = rmass[i];
         } else {
-          itype = type[i];
+          mi = mass[type[i]];
-          jtype = type[j];
+          mj = mass[type[j]];
-          meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
-          if (mask[i] & freeze_group_bit) meff = mass[jtype];
-          if (mask[j] & freeze_group_bit) meff = mass[itype];
         }
+        if (fix_rigid) {
+          if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
+          if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
+        }
+
+        meff = mi*mj / (mi+mj);
+        if (mask[i] & freeze_group_bit) meff = mj;
+        if (mask[j] & freeze_group_bit) meff = mi;
+
+        // normal forces = Hookian contact + normal velocity damping
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;

src/accelerator_omp.h

@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -20,9 +20,11 @@
 #ifdef LMP_USER_OMP
 
   void half_nsq_no_newton_omp(class NeighList *);
+  void half_nsq_no_newton_ghost_omp(class NeighList *);
   void half_nsq_newton_omp(class NeighList *);
 
   void half_bin_no_newton_omp(class NeighList *);
+  void half_bin_no_newton_ghost_omp(class NeighList *);
   void half_bin_newton_omp(class NeighList *);
   void half_bin_newton_tri_omp(class NeighList *);
 
@@ -57,9 +59,11 @@
 // needed for compiling Neighbor class when USER-OMP is not installed
 
   void half_nsq_no_newton_omp(class NeighList *) {}
+  void half_nsq_no_newton_ghost_omp(class NeighList *) {}
   void half_nsq_newton_omp(class NeighList *) {}
 
   void half_bin_no_newton_omp(class NeighList *) {}
+  void half_bin_no_newton_ghost_omp(class NeighList *) {}
   void half_bin_newton_omp(class NeighList *) {}
   void half_bin_newton_tri_omp(class NeighList *) {}
 

src/neighbor.cpp

@@ -945,8 +945,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
             else if (includegroup)
               error->all(FLERR,"Neighbor include group not allowed "
                          "with ghost neighbors");
-            // NOTE: missing OMP version of this one
+            else pb = &Neighbor::half_nsq_no_newton_ghost_omp;
-            else pb = &Neighbor::half_nsq_no_newton_ghost;
           } else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp;
         } else if (rq->newton == 1) {
           pb = &Neighbor::half_nsq_newton_omp;
@@ -955,8 +954,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
           else if (includegroup)
             error->all(FLERR,"Neighbor include group not allowed "
                        "with ghost neighbors");
-          // NOTE: missing OMP version of this one
+          else pb = &Neighbor::half_nsq_no_newton_ghost_omp;
-          else pb = &Neighbor::half_nsq_no_newton_ghost;
         }
       } else if (style == BIN) {
         if (rq->newton == 0) {
@@ -965,8 +963,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
             else if (includegroup)
               error->all(FLERR,"Neighbor include group not allowed "
                          "with ghost neighbors");
-            // NOTE: missing OMP version of this one
+            else pb = &Neighbor::half_bin_no_newton_ghost_omp;
-            else pb = &Neighbor::half_bin_no_newton_ghost;
           } else if (triclinic == 0) {
             pb = &Neighbor::half_bin_newton_omp;
           } else if (triclinic == 1) 
@@ -979,8 +976,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
           else if (includegroup)
             error->all(FLERR,"Neighbor include group not allowed "
                        "with ghost neighbors");
-          // NOTE: missing OMP version of this one
+          else pb = &Neighbor::half_bin_no_newton_ghost_omp;
-          else pb = &Neighbor::half_bin_no_newton_ghost;
         }
       } else if (style == MULTI) {
         if (rq->ghost == 1)

src/region.cpp

@@ -379,4 +379,8 @@ void Region::options(int narg, char **arg)
 
   if (moveflag || rotateflag) dynamic = 1;
   else dynamic = 0;
+
+  // initialize option variables in case region is used between runs
+
+  init();
 }