diff --git a/src/FLD/pair_lubricate.cpp b/src/FLD/pair_lubricate.cpp index 19f2b0d8ec..ce6c00c55a 100644 --- a/src/FLD/pair_lubricate.cpp +++ b/src/FLD/pair_lubricate.cpp @@ -108,8 +108,6 @@ void PairLubricate::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - int overlaps = 0; - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -299,10 +297,6 @@ void PairLubricate::compute(int eflag, int vflag) h_sep = r - 2.0*radi; - // check for overlaps - - if (h_sep < 0.0) overlaps++; - // if less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -445,16 +439,6 @@ void PairLubricate::compute(int eflag, int vflag) } } - // to DEBUG: set print_overlaps to 1 - - int print_overlaps = 0; - if (print_overlaps) { - int overlaps_all; - MPI_Allreduce(&overlaps,&overlaps_all,1,MPI_INT,MPI_SUM,world); - if (overlaps_all && comm->me == 0) - printf("Number of overlaps = %d\n",overlaps); - } - if (vflag_fdotr) virial_fdotr_compute(); } diff --git a/src/FLD/pair_lubricateU.cpp b/src/FLD/pair_lubricateU.cpp index 146b033aa8..2918a803cc 100644 --- a/src/FLD/pair_lubricateU.cpp +++ b/src/FLD/pair_lubricateU.cpp @@ -181,21 +181,17 @@ void PairLubricateU::compute(int eflag, int vflag) void PairLubricateU::stage_one() { - int i,j,ii,inum,itype; + int i,j,ii,inum; double **x = atom->x; double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; - double *radius = atom->radius; - int *type = atom->type; int newton_pair = force->newton_pair; int *ilist; - double radi; - inum = list->inum; ilist = list->ilist; @@ -358,9 +354,6 @@ void PairLubricateU::stage_one() for (ii=0;iiv; double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; - double *radius = atom->radius; - int *type = atom->type; int newton_pair = force->newton_pair; int *ilist; - double radi; - inum = list->inum; ilist = list->ilist; @@ -553,9 +542,6 @@ void PairLubricateU::stage_two(double **x) for (ii=0;iinghost; int newton_pair = force->newton_pair; + double radi; + double vxmu2f = force->vxmu2f; - int overlaps = 0; double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh; inum = list->inum; @@ -750,10 +737,6 @@ void PairLubricateU::compute_Fh(double **x) h_sep = r - 2.0*radi; - // check for overlaps - - if (h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -823,7 +806,6 @@ void PairLubricateU::compute_RU() int newton_pair = force->newton_pair; double vxmu2f = force->vxmu2f; - int overlaps = 0; double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu; inum = list->inum; @@ -954,10 +936,6 @@ void PairLubricateU::compute_RU() h_sep = r - 2.0*radi; - // check for overlaps - - if(h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -1099,7 +1077,6 @@ void PairLubricateU::compute_RU(double **x) int newton_pair = force->newton_pair; double vxmu2f = force->vxmu2f; - int overlaps = 0; double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu; inum = list->inum; @@ -1230,10 +1207,6 @@ void PairLubricateU::compute_RU(double **x) h_sep = r - 2.0*radi; - // check for overlaps - - if(h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -1376,8 +1349,7 @@ void PairLubricateU::compute_RE() int newton_pair = force->newton_pair; double vxmu2f = force->vxmu2f; - int overlaps = 0; - double xl[3],a_sq,a_sh,a_pu; + double xl[3],a_sq,a_sh; inum = list->inum; ilist = list->ilist; @@ -1421,10 +1393,6 @@ void PairLubricateU::compute_RE() h_sep = r - 2.0*radi; - // check for overlaps - - if(h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -1445,7 +1413,6 @@ void PairLubricateU::compute_RE() if (flaglog) { a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } // Relative velocity at the point of closest approach due to Ef only @@ -1517,16 +1484,10 @@ void PairLubricateU::compute_RE() torque[j][1] -= vxmu2f*ty; torque[j][2] -= vxmu2f*tz; } - - // NOTE No a_pu term needed as they add up to zero } } } } - - int print_overlaps = 0; - if (print_overlaps && overlaps) - printf("Number of overlaps=%d\n",overlaps); } /* ---------------------------------------------------------------------- @@ -1553,8 +1514,7 @@ void PairLubricateU::compute_RE(double **x) int newton_pair = force->newton_pair; double vxmu2f = force->vxmu2f; - int overlaps = 0; - double xl[3],a_sq,a_sh,a_pu; + double xl[3],a_sq,a_sh; if (!flagHI) return; @@ -1598,10 +1558,6 @@ void PairLubricateU::compute_RE(double **x) h_sep = r - 2.0*radi; - // check for overlaps - - if(h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -1622,7 +1578,6 @@ void PairLubricateU::compute_RE(double **x) if (flaglog) { a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } // Relative velocity at the point of closest approach due to Ef only @@ -1694,16 +1649,10 @@ void PairLubricateU::compute_RE(double **x) torque[j][1] -= vxmu2f*ty; torque[j][2] -= vxmu2f*tz; } - - // NOTE No a_pu term needed as they add up to zero } } } } - - int print_overlaps = 0; - if (print_overlaps && overlaps) - printf("Number of overlaps=%d\n",overlaps); } @@ -1829,13 +1778,12 @@ void PairLubricateU::init_style() // require that atom radii are identical within each type // require monodisperse system with same radii for all types - double radi, radtype; + double radtype; for (int i = 1; i <= atom->ntypes; i++) { if (!atom->radius_consistency(i,radtype)) error->all(FLERR,"Pair lubricateU requires monodisperse particles"); if (i > 1 && radtype != rad) error->all(FLERR,"Pair lubricateU requires monodisperse particles"); - radi = radtype; } // check for fix deform, if exists it must use "remap v" @@ -2029,20 +1977,13 @@ void PairLubricateU::copy_vec_uo(int inum, double *xcg, { int i,j,ii; int *ilist; - int itype; - double radi; - double inertia; double *rmass = atom->rmass; - int *type = atom->type; ilist = list->ilist; for (ii=0;iiradius[i]; - inertia = 0.4*rmass[i]*radi*radi; for (j=0;j<3;j++) { v[i][j] = xcg[6*ii+j]; diff --git a/src/FLD/pair_lubricateU_poly.cpp b/src/FLD/pair_lubricateU_poly.cpp index 545f1e3207..f78d149522 100644 --- a/src/FLD/pair_lubricateU_poly.cpp +++ b/src/FLD/pair_lubricateU_poly.cpp @@ -154,7 +154,6 @@ void PairLubricateUPoly::iterate(double **x, int stage) int inum = list->inum; int *ilist = list->ilist; - int *type = atom->type; int newton_pair = force->newton_pair; double alpha,beta; @@ -164,7 +163,6 @@ void PairLubricateUPoly::iterate(double **x, int stage) double **v = atom->v; double **f = atom->f; double **omega = atom->omega; - double **angmom = atom->angmom; double **torque = atom->torque; // First compute R_FE*E @@ -329,7 +327,6 @@ void PairLubricateUPoly::compute_Fh(double **x) int nlocal = atom->nlocal; int nghost = atom->nghost; int newton_pair = force->newton_pair; - int overlaps = 0; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; double rsq,r,h_sep,radi,radj; @@ -511,10 +508,6 @@ void PairLubricateUPoly::compute_Fh(double **x) // Find the scalar resistances a_sq, a_sh and a_pu - // check for overlaps - - if (h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -606,7 +599,6 @@ void PairLubricateUPoly::compute_RU(double **x) int *type = atom->type; int nlocal = atom->nlocal; int nghost = atom->nghost; - int overlaps = 0; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; double rsq,r,radi,radj,h_sep; @@ -617,7 +609,6 @@ void PairLubricateUPoly::compute_RU(double **x) double **v = atom->v; double **f = atom->f; double **omega = atom->omega; - double **angmom = atom->angmom; double **torque = atom->torque; double *radius = atom->radius; @@ -760,10 +751,6 @@ void PairLubricateUPoly::compute_RU(double **x) h_sep = r - radi-radj; - // check for overlaps - - if(h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -913,7 +900,6 @@ void PairLubricateUPoly::compute_RE(double **x) int *ilist,*jlist,*numneigh,**firstneigh; int *type = atom->type; - int overlaps = 0; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; double rsq,r,h_sep,radi,radj; @@ -930,7 +916,6 @@ void PairLubricateUPoly::compute_RE(double **x) double vxmu2f = force->vxmu2f; double a_sq = 0.0; double a_sh = 0.0; - double a_pu = 0.0; if (!flagHI) return; @@ -976,10 +961,6 @@ void PairLubricateUPoly::compute_RE(double **x) h_sep = r - radi-radj; - // check for overlaps - - if (h_sep < 0.0) overlaps++; - // If less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -1015,28 +996,6 @@ void PairLubricateUPoly::compute_RE(double **x) a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0 +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2] +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13; - - a_pu = pre[1]*((0.4*beta[0][1]+0.1*beta[0][2])*beta[1][1] - +(0.128-0.132*beta[0][1]+0.332*beta[0][2] - +0.172*beta[0][3])*h_sep)*log_h_sep_beta13; - - /*//a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep)); - a_sq = beta0*beta0/beta1/beta1/h_sep - +(1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*lhsep; - a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) - +pow(beta0,4))/21.0/pow(beta1,4)*h_sep*lhsep; - a_sq *= 6.0*MY_PI*mu*radi; - - a_sh = 4.0*beta0*(2.0+beta0 - +2.0*beta0*beta0)/15.0/pow(beta1,3)*log(1.0/h_sep); - a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) - +16.0*pow(beta0,4))/375.0/pow(beta1,4)*h_sep*log(1.0/h_sep); - a_sh *= 6.0*MY_PI*mu*radi; - - a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); - a_pu += (32.0-33.0*beta0+83.0*beta0*beta0 - +43.0*pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); - a_pu *= 8.0*MY_PI*mu*pow(radi,3);*/ } else a_sq = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep)); @@ -1098,15 +1057,10 @@ void PairLubricateUPoly::compute_RE(double **x) torque[i][1] -= vxmu2f*ty; torque[i][2] -= vxmu2f*tz; - // NOTE No a_pu term needed as they add up to zero } } } } - - int print_overlaps = 0; - if (print_overlaps && overlaps) - printf("Number of overlaps=%d\n",overlaps); } /*----------------------------------------------------------------------- diff --git a/src/FLD/pair_lubricate_poly.cpp b/src/FLD/pair_lubricate_poly.cpp index d37cc4c005..5f368f4786 100644 --- a/src/FLD/pair_lubricate_poly.cpp +++ b/src/FLD/pair_lubricate_poly.cpp @@ -93,8 +93,6 @@ void PairLubricatePoly::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - int overlaps = 0; - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -285,10 +283,6 @@ void PairLubricatePoly::compute(int eflag, int vflag) h_sep = r - radi-radj; - // check for overlaps - - if (h_sep < 0.0) overlaps++; - // if less than the minimum gap use the minimum gap instead if (r < cut_inner[itype][jtype]) @@ -429,16 +423,6 @@ void PairLubricatePoly::compute(int eflag, int vflag) omega[i][2] -= 0.5*h_rate[5]; } } - - // to DEBUG: set print_overlaps to 1 - - int print_overlaps = 0; - if (print_overlaps) { - int overlaps_all; - MPI_Allreduce(&overlaps,&overlaps_all,1,MPI_INT,MPI_SUM,world); - if (overlaps_all && comm->me == 0) - printf("Number of overlaps = %d\n",overlaps); - } } /* ----------------------------------------------------------------------