From 945fc5a36a654de8a685ed4532443f6da63e0ef3 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 29 Jun 2012 18:01:14 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8434 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-CG-CMM/Package.sh | 4 +- src/USER-OMP/fix_qeq_comb_omp.cpp | 70 ++++++--- src/USER-OMP/neigh_full_omp.cpp | 8 +- src/USER-OMP/neigh_half_bin_omp.cpp | 143 +++++++++++++++++++ src/USER-OMP/neigh_half_nsq_omp.cpp | 127 +++++++++++++++- src/USER-OMP/pair_gran_hertz_history_omp.cpp | 36 +++-- src/USER-OMP/pair_gran_hooke_history_omp.cpp | 35 +++-- src/USER-OMP/pair_gran_hooke_omp.cpp | 36 +++-- src/accelerator_omp.h | 6 +- src/neighbor.cpp | 12 +- src/region.cpp | 4 + 11 files changed, 415 insertions(+), 66 deletions(-) diff --git a/src/USER-CG-CMM/Package.sh b/src/USER-CG-CMM/Package.sh index 9b6101979d..3176482366 100644 --- a/src/USER-CG-CMM/Package.sh +++ b/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 diff --git a/src/USER-OMP/fix_qeq_comb_omp.cpp b/src/USER-OMP/fix_qeq_comb_omp.cpp index db3c9fcf77..dbf073db53 100644 --- a/src/USER-OMP/fix_qeq_comb_omp.cpp +++ b/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(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; @@ -88,9 +109,9 @@ void FixQEQCombOMP::post_force(int vflag) // more loops for first-time charge equilibrium - iloop = 0; - if (firstflag) loopmax = 5000; - else loopmax = 2000; + iloop = 0; + if (firstflag) loopmax = 500; + else loopmax = 200; // charge-equilibration loop @@ -99,8 +120,8 @@ void FixQEQCombOMP::post_force(int vflag) update->ntimestep); heatpq = 0.05; - qmass = 0.000548580; - dtq = 0.0006; + qmass = 0.016; + 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]; + 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; @@ -143,13 +173,15 @@ void FixQEQCombOMP::post_force(int vflag) if (me == 0 && fp) fprintf(fp," iteration: %d, enegtot %.6g, " "enegmax %.6g, fq deviation: %.6g\n", - iloop,enegtot,enegmax,enegchk); - - for (i = 0; i < nlocal; i++) + iloop,enegtot,enegmax,enegchk); + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; if (mask[i] & groupbit) - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; - } - + q1[i] += qf[i]*dtq2 - heatpq*q1[i]; + } + } + if (me == 0 && fp) { if (iloop == loopmax) fprintf(fp,"Charges did not converge in %d iterations\n",iloop); diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp index 21643bd29e..8ba2ffaf2b 100644 --- a/src/USER-OMP/neigh_full_omp.cpp +++ b/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]) - neighptr[n++] = j; + if (rsq <= cutneighghostsq[itype][jtype]) 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]) - neighptr[n++] = j; + if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j; } } } diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp index 42205d2400..45dd79d092 100644 --- a/src/USER-OMP/neigh_half_bin_omp.cpp +++ b/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 diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp index 024ba4c212..44b777273e 100644 --- a/src/USER-OMP/neigh_half_nsq_omp.cpp +++ b/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 diff --git a/src/USER-OMP/pair_gran_hertz_history_omp.cpp b/src/USER-OMP/pair_gran_hertz_history_omp.cpp index 7c30321735..b1fc4feabf 100644 --- a/src/USER-OMP/pair_gran_hertz_history_omp.cpp +++ b/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]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + mi = rmass[i]; + mj = rmass[j]; } else { - itype = type[i]; - jtype = 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]; + mi = mass[type[i]]; + mj = mass[type[j]]; } + 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; diff --git a/src/USER-OMP/pair_gran_hooke_history_omp.cpp b/src/USER-OMP/pair_gran_hooke_history_omp.cpp index 7fcfd5d0da..1996feb05f 100644 --- a/src/USER-OMP/pair_gran_hooke_history_omp.cpp +++ b/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]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + mi = rmass[i]; + mj = rmass[j]; } else { - jtype = 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]; + mi = mass[type[i]]; + mj = mass[type[j]]; } + 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; diff --git a/src/USER-OMP/pair_gran_hooke_omp.cpp b/src/USER-OMP/pair_gran_hooke_omp.cpp index 519384105e..eff2aeb0a6 100644 --- a/src/USER-OMP/pair_gran_hooke_omp.cpp +++ b/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]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + mi = rmass[i]; + mj = rmass[j]; } else { - itype = type[i]; - jtype = 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]; + mi = mass[type[i]]; + mj = mass[type[j]]; } + 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; diff --git a/src/accelerator_omp.h b/src/accelerator_omp.h index 89db35cae2..5f224776e2 100644 --- a/src/accelerator_omp.h +++ b/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 *) {} diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 386d962694..100ba8a8a6 100644 --- a/src/neighbor.cpp +++ b/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; + else pb = &Neighbor::half_nsq_no_newton_ghost_omp; } 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; + else pb = &Neighbor::half_nsq_no_newton_ghost_omp; } } 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; + else pb = &Neighbor::half_bin_no_newton_ghost_omp; } 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; + else pb = &Neighbor::half_bin_no_newton_ghost_omp; } } else if (style == MULTI) { if (rq->ghost == 1) diff --git a/src/region.cpp b/src/region.cpp index 9ce1d9fd2a..a9d27e1a38 100644 --- a/src/region.cpp +++ b/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(); }